//
// spatial_orientation_tensors.h
//
// @copyright Deepansh J. Srivastava, 2019-2020.
// Created by Deepansh J. Srivastava, Aug 26, 2019.
// Contact email = srivastava.89@osu.edu
//
#include "mrsimulator.h"
/**
* The scaled spatial orientation tensors (sSOT) from the first-order
* perturbation expansion of the nuclear shielding Hamiltonian,
* in the principal axis system (PAS), includes the
* zeroth and second-rank irreducible tensors which follows,
* @f[ \left.
* \varsigma_{0,0}^{(\sigma)} = \omega_0\delta_\text{iso}
* \right\} \text{Rank-0},
* @f]
* @f[ \left.
* \begin{aligned}
* \varsigma_{2,0}^{(\sigma)} &= -\omega_0\zeta_\sigma, \\
* \varsigma_{2,\pm1}^{(\sigma)} &= 0, \\
* \varsigma_{2,\pm2}^{(\sigma)} &= \frac{1}{\sqrt{6}}\omega_0
* \eta_\sigma \zeta_\sigma,
* \end{aligned}
* \right\} \text{Rank-2},
* @f]
* where @f$\sigma_\text{iso}@f$ is the isotropic nuclear shielding, and,
* @f$\zeta_\sigma@f$, @f$\eta_\sigma@f$ are the shielding anisotropy and
* asymmetry parameters from the symmetric second-rank irreducible nuclear
* shielding tensor defined using Haeberlen convention. Here,
* @f$\omega_0 = -\gamma_I B_0@f$ is the Larmor frequency where, @f$\gamma_I@f$
* and @f$B_0@f$ are the gyromagnetic ratio of the nucleus and the magnetic
* flux density of the external magnetic field, respectively.
*
* For non-zero Euler angles, @f$\Theta = [\alpha, \beta, \gamma]@f$,
* Wigner rotation of @f$\varsigma_{2,n}^{(\sigma)}@f$ is performed
* following,
* @f[ \mathcal{R'}_{2,n}^{(\sigma)}(\Theta) =
* \sum_{m = -2}^2 D^2_{m, n}(\Theta)
* \varsigma_{2,n}^{(\sigma)},
* @f]
* where @f$\mathcal{R'}_{2,n}^{(\sigma)}(\Theta)@f$ are the tensors in the
* frame defined by the Euler angles @f$\Theta@f$.
*
* @note
* - The method accepts a frequency physical quantity, that is,
* @f$\omega_0\sigma_\text{iso}/2\pi@f$ and
* @f$\omega_0\zeta_\sigma/2\pi@f$, as the isotropic
* nuclear shielding and nuclear shielding anisotropy, respectively.
* - When @f$\Theta = [0,0,0]@f$,
* @f$\mathcal{R'}_{2,n}^{(\sigma)}(\Theta) = \varsigma_{2,n}^{(\sigma)}@f$
* where @f$ n \in [-2,2]@f$.
* - @f$\mathcal{R'}_{0,0}^{(\sigma)}(\Theta) = \varsigma_{0,0}^{(\sigma)} ~~~
* \forall ~ \Theta@f$.
* - The method returns @f$\mathcal{R'}_{0,0}^{(\sigma)}(\Theta)/2\pi@f$ and
* @f$\mathcal{R'}_{2,n}^{(\sigma)}(\Theta)/2\pi@f$, that is, in **units of
* frequency**.
*
* @param R_0 A pointer to an array of length 1 where the zeroth-rank
* irreducible tensor,
* @f$\mathcal{R'}_{0,0}^{(\sigma)}(\Theta)/2\pi@f$, will be stored.
*
* @param R_2 A pointer to a complex array of length 5 where the second-rank
* irreducible tensor,
* @f$\mathcal{R'}_{2,n}^{(\sigma)}(\Theta)/2\pi@f$, will be stored ordered
* according to
* @f$\left[\mathcal{R'}_{2,n}^{(\sigma)}(\Theta)/2\pi\right]_{n=-2}^2@f$.
*
* @param omega_0_delta_iso_in_Hz The quantity,
* @f$\omega_0\sigma_\text{iso}/2\pi@f$, given in Hz.
*
* @param omega_0_zeta_sigma_in_Hz The quantity,
* @f$\omega_0\zeta_\sigma/2\pi@f$, given in Hz.
*
* @param eta The nuclear shielding asymmetry, @f$\eta_\sigma \in [0, 1]@f$.
*
* @param Theta A pointer to an array of length 3 where Euler angles,
* ordered as @f$[\alpha, \beta, \gamma]@f$, are stored in radians.
*/
static inline void sSOT_1st_order_nuclear_shielding_tensor_components(
double *restrict R_0, void *restrict R_2,
const double omega_0_delta_iso_in_Hz, const double omega_0_zeta_sigma_in_Hz,
const double eta, const double *Theta) {
// contribution from the zeroth rank
*R_0 = omega_0_delta_iso_in_Hz;
// contribution from the shielding symmetric second rank
vm_double_zeros(10, (double *)R_2);
double *R_2_ = (double *)R_2;
double temp = 0.4082482905 * (omega_0_zeta_sigma_in_Hz * eta);
R_2_[0] = temp; // R2-2 real
R_2_[4] = -omega_0_zeta_sigma_in_Hz; // R2 0 real
R_2_[8] = temp; // R2 2 real
// wigner rotations
if (Theta[0] == 0.0 && Theta[1] == 0.0 && Theta[2] == 0.0) {
return;
}
single_wigner_rotation(2, Theta, R_2, R_2);
}
/**
* The scaled spatial orientation tensors (sSOT) from the first-order
* perturbation expansion of the electric quadrupole Hamiltonian, in the
* principal axis system (PAS), includes the second-rank irreducible tensor
* which follows,
* @f[ \left.
* \begin{aligned}
* \varsigma_{2,0}^{(q)} &= \frac{1}{\sqrt{6}} \omega_q, \\
* \varsigma_{2,\pm1}^{(q)} &= 0, \\
* \varsigma_{2,\pm2}^{(q)} &= -\frac{1}{6} \eta_q \omega_q,
* \end{aligned}
* \right\} \text{Rank-2},
* @f]
* where @f$\omega_q = \frac{6\pi C_q}{2I(2I-1)}@f$ is the quadrupole splitting
* frequency, and @f$\eta_q@f$ is the quadrupole asymmetry parameter. Here,
* @f$I@f$ is the spin quantum number of the quadrupole nucleus, and @f$C_q@f$
* is the quadrupole coupling constant.
*
* As before, for non-zero Euler angles, @f$\Theta = [\alpha,\beta,\gamma]@f$,
* a Wigner rotation of @f$\varsigma_{2,n}^{(q)}@f$ is performed following,
* @f[ \mathcal{R'}_{2,n}^{(q)}(\Theta) = \sum_{m = -2}^2 D^2_{m, n}(\Theta)
* \varsigma_{2,n}^{(q)}.
* @f]
* where @f$\mathcal{R'}_{2,n}^{(q)}(\Theta)@f$ are the tensors in the
* frame defined by the Euler angles @f$\Theta@f$.
*
* @note
* - When @f$\Theta = [0,0,0]@f$,
* @f$\mathcal{R'}_{2,n}^{(q)}(\Theta) = \varsigma_{2,n}^{(q)}@f$
* where @f$ n \in [-2,2]@f$.
* - The method returns @f$\mathcal{R'}_{2,0}^{(q)}(\Theta)/2\pi@f$, that is,
* in **units of frequency**.
*
* @param R_2 A pointer to a complex array of length 5 where the the
* second-rank irreducible tensor,
* @f$\mathcal{R'}_{2,n}^{(q)}(\Theta)/2\pi@f$, will be stored ordered
* according to
* @f$\left[\mathcal{R'}_{2,n}^{(q)}(\Theta)/2\pi\right]_{n=-2}^2@f$.
*
* @param spin The spin quantum number, @f$I@f$.
*
* @param Cq_in_Hz The quadrupole coupling constant, @f$C_q@f$, in Hz.
*
* @param eta The quadrupole asymmetry parameter, @f$\eta_q \in [0, 1]@f$.
*
* @param Theta A pointer to an array of length 3 where Euler angles,
* ordered as @f$[\alpha, \beta, \gamma]@f$, are stored in radians.
*/
static inline void sSOT_1st_order_electric_quadrupole_tensor_components(
void *restrict R_2, const double spin, const double Cq_in_Hz,
const double eta, const double *Theta) {
/* vq is the Quadrupole coupling constant given as vq = 3*Cq/(2I(2I-1)),
where `I` is the spin quantum number. */
double vq = 3.0 * Cq_in_Hz;
double denominator = 2.0 * spin * (2.0 * spin - 1.0);
vq /= denominator;
// contribution from the symmetric second rank quadrupole tensor
vm_double_zeros(10, (double *)R_2);
double *R_2_ = (double *)R_2;
double temp = -0.1666666667 * (vq * eta);
R_2_[0] = temp; // R2-2 real
R_2_[4] = 0.4082482905 * vq; // R2 0 real
R_2_[8] = temp; // R2 2 real
// wigner rotations
if (Theta[0] == 0.0 && Theta[1] == 0.0 && Theta[2] == 0.0) {
return;
}
single_wigner_rotation(2, Theta, R_2, R_2);
}
/**
* The scaled spatial orientation tensors (sSOT) from the second-order
* perturbation expansion of the electric quadrupole Hamiltonian, in the
* principal axis system (PAS), includes the zeroth, second and fourth-rank
* irreducible tensors which follows,
* @f[\left.
* \varsigma_{0,0}^{(qq)} = \frac{\omega_q^2}{\omega_0} \frac{1}{6\sqrt{5}}
* \left(\frac{\eta_q^2}{3} + 1 \right)
* \right\} \text{Rank-0},
* @f]
* @f[ \left.
* \begin{aligned}
* \varsigma_{2,0}^{(qq)} &= \frac{\omega_q^2}{\omega_0}
* \frac{\sqrt{2}}{6\sqrt{7}}
* \left(\frac{\eta_q^2}{3} - 1 \right), \\
* \varsigma_{2,\pm1}^{(qq)} &= 0, \\
* \varsigma_{2,\pm2}^{(qq)} &= -\frac{\omega_q^2}{\omega_0}
* \frac{1}{3\sqrt{21}} \eta_q,
* \end{aligned}
* \right\} \text{Rank-2},
* @f]
* @f[ \left.
* \begin{aligned}
* \varsigma_{4,0}^{(qq)} &= \frac{\omega_q^2}{\omega_0}
* \frac{1}{\sqrt{70}} \left(\frac{\eta_q^2}{18} + 1 \right), \\
* \varsigma_{4,\pm1}^{(qq)} &= 0, \\
* \varsigma_{4,\pm2}^{(qq)} &= -\frac{\omega_q^2}{\omega_0}
* \frac{1}{6\sqrt{7}} \eta_q, \\
* \varsigma_{4,\pm3}^{(qq)} &= 0, \\
* \varsigma_{4,\pm4}^{(qq)} &= \frac{\omega_q^2}{\omega_0} \frac{1}{36}
* \eta_q^2,
* \end{aligned}
* \right\} \text{Rank-4},
* @f]
* where @f$\omega_q = \frac{6\pi C_q}{2I(2I-1)}@f$ is the quadrupole splitting
* frequency, @f$\omega_0@f$ is the Larmor angular frequency, and @f$\eta_q@f$
* is the quadrupole asymmetry parameter. Here, @f$I@f$ is the spin quantum
* number, and @f$C_q@f$ is the quadrupole coupling constant.
*
* For non-zero Euler angles, @f$\Theta = [\alpha,\beta,\gamma]@f$, Wigner
* rotation of @f$\varsigma_{2,n}^{(qq)}@f$ and
* @f$\varsigma_{4,n}^{(qq)}@f$ are performed following,
* @f[
* \mathcal{R'}_{2,n}^{(qq)}(\Theta) &= \sum_{m = -2}^2 D^2_{m, n}(\Theta)
* \varsigma_{2,n}^{(qq)}, \\
* \mathcal{R'}_{4,n}^{(qq)}(\Theta) &= \sum_{m = -4}^4 D^4_{m, n}(\Theta)
* \varsigma_{4,n}^{(qq)}.
* @f]
* where @f$\mathcal{R'}_{2,n}^{(qq)}(\Theta)@f$ and
* @f$\mathcal{R'}_{4,n}^{(qq)}(\Theta)@f$ are the tensors in the
* frame defined by the Euler angles @f$\Theta@f$.
*
* @note
* - When @f$\Theta = [0,0,0]@f$,
* @f$\mathcal{R'}_{2,n}^{(qq)}(\Theta) = \varsigma_{2,n}^{(qq)}@f$
* where @f$ n \in [-2,2]@f$.
* - When @f$\Theta = [0,0,0]@f$,
* @f$\mathcal{R'}_{4,n}^{(qq)}(\Theta) = \varsigma_{4,n}^{(qq)}@f$
* where @f$ n \in [-4,4]@f$.
* - @f$\mathcal{R'}_{0,0}^{(qq)}(\Theta) = \varsigma_{0,0}^{(qq)} ~~~
* \forall ~ \Theta@f$.
* - The method returns @f$\mathcal{R'}_{0,0}^{(qq)}(\Theta)/2\pi@f$,
* @f$\mathcal{R'}_{2,n}^{(qq)}(\Theta)/2\pi@f$, and
* @f$\mathcal{R'}_{4,n}^{(qq)}(\Theta)/2\pi@f$, that is, in **units of
* frequency**.
*
* @param R_0 A pointer to an array of length 1 where the zeroth-rank
* irreducible tensor, @f$\mathcal{R'}_{0,0}^{(qq)}(\Theta)/2\pi@f$, will
* be stored.
*
* @param R_2 A pointer to a complex array of length 5 where the second-rank
* irreducible tensor,
* @f$\mathcal{R'}_{2,n}^{(qq)}(\Theta)/2\pi@f$, will be stored ordered
* according to
* @f$\left[\mathcal{R'}_{2,n}^{(qq)}(\Theta)/2\pi\right]_{n=-2}^2@f$.
*
* @param R_4 A pointer to a complex array of length 9 where the fourth-rank
* irreducible tensor,
* @f$\mathcal{R'}_{4,n}^{(qq)}(\Theta)/2\pi@f$, will be stored ordered
* according to
* @f$\left[\mathcal{R'}_{4,n}^{(qq)}(\Theta)/2\pi\right]_{n=-4}^4@f$.
*
* @param spin The spin quantum number, @f$I@f$.
*
* @param v0_in_Hz The Larmor frequency, @f$\omega_0/2\pi@f$, in Hz.
*
* @param Cq_in_Hz The quadrupole coupling constant, @f$C_q@f$, in Hz.
*
* @param eta The quadrupole asymmetry parameter, @f$\eta_q \in [0, 1]@f$.
*
* @param Theta A pointer to an array of length 3 where Euler angles,
* ordered as @f$[\alpha, \beta, \gamma]@f$, are stored in radians.
*/
static inline void sSOT_2nd_order_electric_quadrupole_tensor_components(
double *restrict R_0, void *restrict R_2, void *restrict R_4,
const double spin, const double v0_in_Hz, const double Cq_in_Hz,
const double eta, const double *Theta) {
/* vq is the Quadrupole coupling constant given as vq = 3*Cq/(2I(2I-1)),
where `I` is the spin quantum number. */
double vq = 3.0 * Cq_in_Hz;
double denominator = 2.0 * spin * (2.0 * spin - 1.0);
vq /= denominator;
double scale = vq * vq / v0_in_Hz;
double eta2 = eta * eta;
// contribution from the zeroth rank
*R_0 = (eta2 * 0.33333333333 + 1.0) * 0.07453559925 * scale;
// contribution from the second rank
vm_double_zeros(10, (double *)R_2);
double *R_2_ = (double *)R_2;
double temp = -eta * 0.07273929675 * scale;
double temp2 = 0.08908708064 * scale * (eta2 * 0.33333333333 - 1.0);
R_2_[0] = temp; // R2-2 real
R_2_[4] = temp2; // R2 0 real
R_2_[8] = temp; // R2 2 real
// contribution from the fourth rank
vm_double_zeros(18, (double *)R_4);
double *R_4_ = (double *)R_4;
temp = eta2 * 0.02777777778 * scale;
temp2 = -0.06299407883 * eta * scale;
double temp4 = 0.1195228609 * scale * (eta2 * 0.05555555556 + 1.0);
R_4_[0] = temp; // R4-4 real
R_4_[4] = temp2; // R4-2 real
R_4_[8] = temp4; // R4 0 real
R_4_[12] = temp2; // R4 2 real
R_4_[16] = temp; // R4 4 real
// wigner rotations
if (Theta[0] == 0.0 && Theta[1] == 0.0 && Theta[2] == 0.0) {
return;
}
single_wigner_rotation(2, Theta, R_2, R_2);
single_wigner_rotation(4, Theta, R_4, R_4);
}
