Angular Namespace Reference

Detailed Description

Angular provides functions and classes for calculating and storing angular factors (3,6,9-J symbols etc.)

Provides functions to:

• Calculate wigner 3,6,9-J symbols + Clebsh-Gordon coefs etc..
• Also provied look-up tables, which are faster than calculating symbols on-the-fly
• Wrapper functions to calculate wigner 3,6,9-J symbols.
• Uses GSL: https://www.gnu.org/software/gsl/doc/html/specfunc.html#coupling-coefficients

The general equations are defined:

\begin{align} \frac{1}{r_{12}} &= \sum_{kq} \frac{r_<^k}{r_>^{k+1}}(-1)^q C^k_{-q}(\hat{n}_1)C^k_{q}(\hat{n}_2)\\ C^k_{q} &\equiv \sqrt{\frac{4\pi}{2k+1}} Y_{kq}(\hat{n}), \end{align}

and

\begin{align} C^k_{ab} &\equiv \langle{\kappa_a}||C^k||{\kappa_b}\rangle \equiv (-1)^{j_a+1/2} \widetilde C^k_{ab} ,\\ &= (-1)^{j_a+1/2}\sqrt{[j_a][j_b]} \begin{pmatrix} {j_a} & {j_b} & {k} \\ {-1/2} & {1/2} &{0} \end{pmatrix} \pi(l_a+l_b+k). \end{align}

Classes
class	CkTable
	Lookup table for Ck_ab reduced matrix elements and the special 3j symbol (j_a, j_b, k; -1/2, 1/2, 0). More...
class	SixJTable
	Lookup table for Wigner 6j symbols. More...

Functions
constexpr int	jindex_to_twoj (int jindex)
	Converts jindex to 2j; jindex = 0,1,2,... maps to 2j = 1,3,5,...
constexpr int	twoj_to_jindex (int twoj)
	Converts 2j to jindex; 2j = 1,3,5,... maps to jindex = 0,1,2,...
constexpr int	kappa_to_jindex (int ka)
	Converts kappa to jindex; e.g. {-1, 1, -2, 2} -> {0, 0, 1, 1}.
template<class A >
int	twojk (const A &a)
	Returns 2*k if a is an integer, or 2*j if a is a DiracSpinor.
template<class A , class B , class C , class D , class E , class F >
double	SixJ (const A &a, const B &b, const C &c, const D &d, const E &e, const F &f)
	"Magic" 6j symbol accepting integers or DiracSpinors.
constexpr int	l_k (int ka)
	returns l given kappa
constexpr int	twoj_k (int ka)
	returns 2j given kappa
constexpr int	parity_k (int ka)
	returns parity [(-1)^l] given kappa
constexpr int	parity_l (int l)
	returns parity [(-1)^l] given l
constexpr int	l_tilde_k (int ka)
	"Complimentary" l := (2j-l) = l+/-1, for j=l+/-1/2
constexpr int	kappa_twojl (int twoj, int l)
	returns kappa, given 2j and l
constexpr int	kappa_twojpi (const int twoj, const int pi)
	returns kappa, given 2*j and parity (+/-1),
constexpr bool	zeroQ (double x, double eps=1.0e-10)
	Checks if a double is within eps of 0.0 [eps=1.0e-10 by default].
constexpr std::uint64_t	kappa_to_kindex (int ka)
	returns kappa_index, given kappa
constexpr int	kindex_to_kappa (std::uint64_t i)
	Returns kappa, given kappa_index.
constexpr int	kindex_to_twoj (std::uint64_t i)
	returns 2j, given kappa_index
constexpr int	kindex_to_l (std::uint64_t i)
	returns l, given kappa_index
constexpr std::uint64_t	l_to_max_kindex (int l)
	Maximum kappa index for a given orbital angular momentum l.
constexpr std::uint64_t	states_below_n (int n)
	Returns number of possible states below given n (i.e., n' < n)
constexpr std::uint64_t	nk_to_index (int n, int k)
	Returns nk_index, given {n, kappa}.
std::pair< int, int >	index_to_nk (std::uint64_t nk_index)
	Returns {n, kappa} given nk_index.
std::pair< int, std::uint64_t >	nkindex_to_n_kindex (std::uint64_t nk_index)
	Returns {n, kappa_index} given nk_index.
int	nkindex_to_kappa (std::uint64_t nk_index)
	Returns kappa, given nk_index.
int	nkindex_to_twoj (std::uint64_t nk_index)
	Returns 2*j, given nk_index.
int	nkindex_to_l (std::uint64_t nk_index)
	Returns l, given nk_index.
constexpr bool	evenQ (int a)
	Returns true if a is even - for integer values.
constexpr bool	evenQ_2 (int two_a)
	Returns true if a (an int) is even, given 2*a (true if two_a/2 is even)
constexpr int	neg1pow (int a)
	Evaluates (-1)^{a} (for integer a)
constexpr int	neg1pow_2 (int two_a)
	Evaluates (-1)^{two_a/2} (for integer a; two_a is even)
constexpr int	parity (int la, int lb, int k)
	Parity rule. Returns 1 only if la+lb+k is even.
constexpr int	triangle (int j1, int j2, int J)
	Returns 1 if triangle rule is satisfied. nb: works with j OR twoj!
constexpr int	sumsToZero (int m1, int m2, int m3)
	Checks if three ints sum to zero, returns 1 if they do.
double	threej_2 (int two_j1, int two_j2, int two_j3, int two_m1, int two_m2, int two_m3)
	Wigner 3j symbol. Inputs are 2*j (or 2*l) as integers.
double	special_threej_2 (int two_j1, int two_j2, int two_k)
	Special 3j symbol: (j1 j2 k; -1/2 1/2 0). Inputs are 2*j as integers.
double	cg_2 (int two_j1, int two_m1, int two_j2, int two_m2, int two_J, int two_M)
	Clebsch-Gordon coefficient <j1 m1, j2 m2 | J M>. Takes 2*j as integers.
bool	sixj_zeroQ (int a, int b, int c, int d, int e, int f)
	Returns true if the 6j symbol is zero by triangle/parity rules. Inputs are 2*j.
bool	sixjTriads (std::optional< int > a, std::optional< int > b, std::optional< int > c, std::optional< int > d, std::optional< int > e, std::optional< int > f)
	Checks 6j triangle conditions with optional arguments (wildcards). Inputs are 2*j.
double	sixj_2 (int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6)
	Wigner 6j symbol {j1 j2 j3 | j4 j5 j6}. Inputs are 2*j as integers.
double	ninej_2 (int two_j1, int two_j2, int two_j3, int two_j4, int two_j5, int two_j6, int two_j7, int two_j8, int two_j9)
	Wigner 9j symbol. Inputs are 2*j as integers.
double	Ck_kk (int k, int ka, int kb)
	Reduced relativistic angular ME: <ka||C^k||kb>. Takes kappa values.
double	Ck_kk_mmq (int k, int ka, int kb, int twoma, int twomb, int twoq)
	Full (non-reduced) C^k_q matrix element. Takes kappa and 2*m values.
bool	Ck_kk_SR (int k, int ka, int kb)
	Selection rule check for C^k: returns true if <ka||C^k||kb> may be non-zero.
double	tildeCk_kk (int k, int ka, int kb)
	Tilde variant of the C^k reduced ME: (-1)^{ja+1/2} * <ka||C^k||kb>
std::pair< int, int >	kminmax_Ck (int ka, int kb)
	Min and max non-zero multipole rank k for C^k, given kappa_a and kappa_b.
double	S_kk (int ka, int kb)
	Reduced spin-1/2 angular ME: <ka||S||kb>. Takes kappa values.

Function Documentation

jindex_to_twoj()

constexpr int Angular::jindex_to_twoj ( int  jindex )

constexpr

Converts jindex to 2j; jindex = 0,1,2,... maps to 2j = 1,3,5,...

twoj_to_jindex()

constexpr int Angular::twoj_to_jindex ( int  twoj )

constexpr

Converts 2j to jindex; 2j = 1,3,5,... maps to jindex = 0,1,2,...

kappa_to_jindex()

constexpr int Angular::kappa_to_jindex ( int  ka )

constexpr

Converts kappa to jindex; e.g. {-1, 1, -2, 2} -> {0, 0, 1, 1}.

twojk()

template<class A >

int Angular::twojk

(

const A &

a

)

Returns 2*k if a is an integer, or 2*j if a is a DiracSpinor.

SixJ()

template<class A , class B , class C , class D , class E , class F >

double Angular::SixJ	(	const A &	a,
		const B &	b,
		const C &	c,
		const D &	d,
		const E &	e,
		const F &	f
	)

"Magic" 6j symbol accepting integers or DiracSpinors.

Pass an integer for a rank \(k\), or a DiracSpinor for a half-integer \(j\). Example: SixJ(Fa, Fb, k, Fc, Fd, k) evaluates \(\{j_a,j_b,k;j_c,j_d,k\}\).

Warning

Do not pre-multiply by 2, and do not pass \(j\) or \(2j\) directly pass the DiracSpinor itself.

l_k()

constexpr int Angular::l_k ( int  ka )

constexpr

returns l given kappa

twoj_k()

constexpr int Angular::twoj_k ( int  ka )

constexpr

returns 2j given kappa

parity_k()

constexpr int Angular::parity_k ( int  ka )

constexpr

returns parity [(-1)^l] given kappa

parity_l()

constexpr int Angular::parity_l ( int  l )

constexpr

returns parity [(-1)^l] given l

l_tilde_k()

constexpr int Angular::l_tilde_k ( int  ka )

constexpr

"Complimentary" l := (2j-l) = l+/-1, for j=l+/-1/2

kappa_twojl()

constexpr int Angular::kappa_twojl ( int  twoj, int  l  )

constexpr

returns kappa, given 2j and l

kappa_twojpi()

constexpr int Angular::kappa_twojpi ( const int  twoj, const int  pi  )

constexpr

returns kappa, given 2*j and parity (+/-1),

zeroQ()

constexpr bool Angular::zeroQ ( double  x, double  eps = 1.0e-10  )

inlineconstexpr

Checks if a double is within eps of 0.0 [eps=1.0e-10 by default].

kappa_to_kindex()

constexpr std::uint64_t Angular::kappa_to_kindex ( int  ka )

constexpr

returns kappa_index, given kappa

Kappa Index: For easy array access, define 1-to-1 index for each kappa:

kappa	-1	1	-2	2	-3	3	-4	4	...
index	0	1	2	3	4	5	6	7	...

\[ i_k(\kappa)= \begin{cases} -2\kappa-2, & \kappa<0, \\ 2\kappa-1, & \kappa>0 . \end{cases} \]

\[ \kappa(i_k)= \begin{cases} -\left(\frac{i_k}{2}+1\right), & i_k\ \text{even},\\ \frac{i_k+1}{2}, & i_k\ \text{odd}. \end{cases} \]

kindex_to_kappa()

constexpr int Angular::kindex_to_kappa ( std::uint64_t  i )

constexpr

Returns kappa, given kappa_index.

see Angular::kappa_to_kindex()

kindex_to_twoj()

constexpr int Angular::kindex_to_twoj ( std::uint64_t  i )

constexpr

returns 2j, given kappa_index

see Angular::kappa_to_kindex()

kindex_to_l()

constexpr int Angular::kindex_to_l ( std::uint64_t  i )

constexpr

returns l, given kappa_index

see Angular::kappa_to_kindex()

l_to_max_kindex()

constexpr std::uint64_t Angular::l_to_max_kindex ( int  l )

constexpr

Maximum kappa index for a given orbital angular momentum l.

Parameters

l	Orbital angular momentum quantum number

Returns

Kappa index of the state with maximum (larger) j

See also

Angular::kappa_to_kindex()

states_below_n()

constexpr std::uint64_t Angular::states_below_n ( int  n )

constexpr

Returns number of possible states below given n (i.e., n' < n)

This could be made more "compact" with a maximum l. In practice, we go to large n, but never very large l.

Warning

• n must be >1 (no non-standard states)

Cannot overflow; converts to std::uint64_t.

nk_to_index()

constexpr std::uint64_t Angular::nk_to_index ( int  n, int  k  )

constexpr

Returns nk_index, given {n, kappa}.

For convenient array access we define a one-to-one mapping between the pair \((n,\kappa)\) and a non-negative index:

\[ i_{nk}(n,\kappa) = (n-1)^2 + i_k(\kappa), \]

where \(i_k(\kappa)\) is the kappa index defined by Angular::kappa_to_kindex()

Equivalently,

\[ i_{nk}(n,\kappa) = n^2 - 2n + 1 + i_k(\kappa). \]

Since

\[ n^2 - 2n + 1 = (n-1)^2 = \texttt{states\_below\_n}(n), \]

this counts the number of states with principal quantum number \(n' < n\), plus the index within the \(n\) shell.

Note

This indexing is only valid for physical bound-state quantum numbers, n >= 1, and cannot in general be used for arbitrary basis states.

Warning

To safely convert the returned index (uint64_t) to int, one must have n <= 46340.

To safely convert the returned index to Coulomb::nkIndex (uint16_t), one must have n <= 256.

index_to_nk()

std::pair< int, int > Angular::index_to_nk ( std::uint64_t  nk_index )

inline

Returns {n, kappa} given nk_index.

Inverse of Angular::nk_to_index().

nkindex_to_n_kindex()

std::pair< int, std::uint64_t > Angular::nkindex_to_n_kindex ( std::uint64_t  nk_index )

inline

Returns {n, kappa_index} given nk_index.

Inverse of Angular::nk_to_index().

nkindex_to_kappa()

int Angular::nkindex_to_kappa ( std::uint64_t  nk_index )

inline

Returns kappa, given nk_index.

See Angular::nk_to_index().

nkindex_to_twoj()

int Angular::nkindex_to_twoj ( std::uint64_t  nk_index )

inline

Returns 2*j, given nk_index.

See Angular::nk_to_index().

nkindex_to_l()

int Angular::nkindex_to_l ( std::uint64_t  nk_index )

inline

Returns l, given nk_index.

See Angular::nk_to_index().

evenQ()

constexpr bool Angular::evenQ ( int  a )

constexpr

Returns true if a is even - for integer values.

evenQ_2()

constexpr bool Angular::evenQ_2 ( int  two_a )

constexpr

Returns true if a (an int) is even, given 2*a (true if two_a/2 is even)

neg1pow()

constexpr int Angular::neg1pow ( int  a )

constexpr

Evaluates (-1)^{a} (for integer a)

neg1pow_2()

constexpr int Angular::neg1pow_2 ( int  two_a )

constexpr

Evaluates (-1)^{two_a/2} (for integer a; two_a is even)

parity()

constexpr int Angular::parity ( int  la, int  lb, int  k  )

constexpr

Parity rule. Returns 1 only if la+lb+k is even.

triangle()

constexpr int Angular::triangle ( int  j1, int  j2, int  J  )

constexpr

Returns 1 if triangle rule is satisfied. nb: works with j OR twoj!

sumsToZero()

constexpr int Angular::sumsToZero ( int  m1, int  m2, int  m3  )

constexpr

Checks if three ints sum to zero, returns 1 if they do.

threej_2()

double Angular::threej_2 ( int  two_j1, int  two_j2, int  two_j3, int  two_m1, int  two_m2, int  two_m3  )

inline

Wigner 3j symbol. Inputs are 2*j (or 2*l) as integers.

Computes the Wigner 3j symbol:

\[ \begin{pmatrix} j_1 & j_2 & j_3 \\ m_1 & m_2 & m_3 \end{pmatrix} \]

Returns zero if the triangle rule or \( m_1+m_2+m_3=0 \) is not satisfied. Wraps gsl_sf_coupling_3j.

Parameters

two_j1	2*j1
two_j2	2*j2
two_j3	2*j3
two_m1	2*m1
two_m2	2*m2
two_m3	2*m3

Returns

Value of the 3j symbol; 0 if selection rules are not met.

special_threej_2()

double Angular::special_threej_2 ( int  two_j1, int  two_j2, int  two_k  )

inline

Special 3j symbol: (j1 j2 k; -1/2 1/2 0). Inputs are 2*j as integers.

Evaluates the commonly occurring 3j symbol:

\[ \begin{pmatrix} j_1 & j_2 & k \\ -\tfrac{1}{2} & \tfrac{1}{2} & 0 \end{pmatrix} \]

Handles the \( k=0 \) case analytically; otherwise wraps gsl_sf_coupling_3j. Returns zero if the triangle rule is not satisfied.

Parameters

two_j1	2*j1
two_j2	2*j2
two_k	2*k

Returns

Value of the 3j symbol; 0 if selection rules are not met.

cg_2()

double Angular::cg_2 ( int  two_j1, int  two_m1, int  two_j2, int  two_m2, int  two_J, int  two_M  )

inline

Clebsch-Gordon coefficient <j1 m1, j2 m2 | J M>. Takes 2*j as integers.

Computes the Clebsch-Gordon coefficient:

\[ \langle j_1 m_1,\, j_2 m_2 \,|\, J M \rangle = (-1)^{j_1-j_2+M}\sqrt{2J+1} \begin{pmatrix} j_1 & j_2 & J \\ m_1 & m_2 & -M \end{pmatrix} \]

Returns zero if the triangle rule or \( m_1+m_2=M \) is not satisfied.

Parameters

two_j1	2*j1
two_m1	2*m1
two_j2	2*j2
two_m2	2*m2
two_J	2*J
two_M	2*M

Returns

Value of the CG coefficient; 0 if selection rules are not met.

sixj_zeroQ()

bool Angular::sixj_zeroQ ( int  a, int  b, int  c, int  d, int  e, int  f  )

inline

Returns true if the 6j symbol is zero by triangle/parity rules. Inputs are 2*j.

Checks all four triangle triads and integer-sum (parity) conditions for:

\[ \sixj{a}{b}{c}{d}{e}{f} \]

The triads (abc), (cde), (bdf), (efa) must each satisfy the triangle rule and sum to an even integer (inputs are 2*j).

Note

Some symbols that pass these checks are still numerically zero. GSL may return a very small non-zero value in such cases.

Warning

Only valid for inputs that are 2*(integer or half-integer).

sixjTriads()

bool Angular::sixjTriads ( std::optional< int >  a, std::optional< int >  b, std::optional< int >  c, std::optional< int >  d, std::optional< int >  e, std::optional< int >  f  )

inline

Checks 6j triangle conditions with optional arguments (wildcards). Inputs are 2*j.

Checks only the triads formed by the provided (non-null) arguments. Missing (empty) optionals act as wildcards. For example, sixjTriads(a,b,c,{},{},{}) checks whether any 6j symbol \( \{ a\, b\, c \;|\; *\, *\, * \} \) could be non-zero.

The four triads of \( \{a\,b\,c\;|\;d\,e\,f\} \) are: (abc), (cde), (aef), (bdf).

Parameters

a	2*j1 (optional)
b	2*j2 (optional)
c	2*j3 (optional)
d	2*j4 (optional)
e	2*j5 (optional)
f	2*j6 (optional)

Returns

false if any complete triad fails the triangle rule; true otherwise.

sixj_2()

double Angular::sixj_2 ( int  two_j1, int  two_j2, int  two_j3, int  two_j4, int  two_j5, int  two_j6  )

inline

Wigner 6j symbol {j1 j2 j3 | j4 j5 j6}. Inputs are 2*j as integers.

Computes the Wigner 6j symbol:

\[ \sixj{j_1}{j_2}{j_3}{j_4}{j_5}{j_6} \]

Returns zero if any triangle or parity condition is violated (via sixj_zeroQ). Wraps gsl_sf_coupling_6j.

Parameters

two_j1	2*j1
two_j2	2*j2
two_j3	2*j3
two_j4	2*j4
two_j5	2*j5
two_j6	2*j6

Returns

Value of the 6j symbol; 0 if selection rules are not met.

ninej_2()

double Angular::ninej_2 ( int  two_j1, int  two_j2, int  two_j3, int  two_j4, int  two_j5, int  two_j6, int  two_j7, int  two_j8, int  two_j9  )

inline

Wigner 9j symbol. Inputs are 2*j as integers.

Computes the Wigner 9j symbol:

\[ \begin{Bmatrix} j_1 & j_2 & j_3 \\ j_4 & j_5 & j_6 \\ j_7 & j_8 & j_9 \end{Bmatrix} \]

Wraps gsl_sf_coupling_9j.

Parameters

two_j1	2*j1
two_j2	2*j2
two_j3	2*j3
two_j4	2*j4
two_j5	2*j5
two_j6	2*j6
two_j7	2*j7
two_j8	2*j8
two_j9	2*j9

Returns

Value of the 9j symbol.

Ck_kk()

double Angular::Ck_kk ( int  k, int  ka, int  kb  )

inline

Reduced relativistic angular ME: <ka||C^k||kb>. Takes kappa values.

Computes the reduced matrix element of the spherical tensor \( C^k \) between relativistic states labelled by Dirac quantum numbers:

\[ \langle \kappa_a \| C^k \| \kappa_b \rangle = (-1)^{j_a+\tfrac{1}{2}} \sqrt{[j_a][j_b]} \begin{pmatrix} j_a & j_b & k \\ -\tfrac{1}{2} & \tfrac{1}{2} & 0 \end{pmatrix} \pi(l_a+l_b+k), \]

where \( [j] \equiv 2j+1 \) and \( \pi(l_a+l_b+k) \) is 1 if \( l_a+l_b+k \) is even, 0 otherwise.

Parameters

k	Multipole rank.
ka	Dirac quantum number \( \kappa_a \) (bra).
kb	Dirac quantum number \( \kappa_b \) (ket).

Returns

\( \langle \kappa_a \| C^k \| \kappa_b \rangle \); 0 if selection rules are not met.

Ck_kk_mmq()

double Angular::Ck_kk_mmq ( int  k, int  ka, int  kb, int  twoma, int  twomb, int  twoq  )

inline

Full (non-reduced) C^k_q matrix element. Takes kappa and 2*m values.

Computes the full matrix element including magnetic quantum numbers via the Wigner-Eckart theorem:

\[ \langle \kappa_a m_a | C^k_q | \kappa_b m_b \rangle = (-1)^{j_a-m_a} \begin{pmatrix} j_a & k & j_b \\ -m_a & q & m_b \end{pmatrix} \langle \kappa_a \| C^k \| \kappa_b \rangle \]

Parameters

k	Multipole rank.
ka	Dirac quantum number \( \kappa_a \) (bra).
kb	Dirac quantum number \( \kappa_b \) (ket).
twoma	2*m_a
twomb	2*m_b
twoq	2*q

Returns

\( \langle \kappa_a m_a | C^k_q | \kappa_b m_b \rangle \)

Ck_kk_SR()

bool Angular::Ck_kk_SR ( int  k, int  ka, int  kb  )

inline

Selection rule check for C^k: returns true if <ka||C^k||kb> may be non-zero.

Checks parity and triangle conditions only; does not compute the value. Equivalent to Ck_kk(k,ka,kb)!=0 but avoids the full computation.

Parameters

k	Multipole rank.
ka	Dirac quantum number \( \kappa_a \).
kb	Dirac quantum number \( \kappa_b \).

Returns

false if \( \langle \kappa_a \| C^k \| \kappa_b \rangle = 0 \) by selection rules; true otherwise.

tildeCk_kk()

double Angular::tildeCk_kk ( int  k, int  ka, int  kb  )

inline

Tilde variant of the C^k reduced ME: (-1)^{ja+1/2} * <ka||C^k||kb>

Computes \( \widetilde{C}^k_{ab} \equiv (-1)^{j_a+1/2} \langle \kappa_a \| C^k \| \kappa_b \rangle \). This removes the \( (-1)^{j_a+1/2} \) phase from Ck_kk, leaving the purely geometric 3j piece times \( \sqrt{[j_a][j_b]} \).

Parameters

k	Multipole rank.
ka	Dirac quantum number \( \kappa_a \).
kb	Dirac quantum number \( \kappa_b \).

Returns

\( \widetilde{C}^k_{ab} \)

kminmax_Ck()

std::pair< int, int > Angular::kminmax_Ck ( int  ka, int  kb  )

inline

Min and max non-zero multipole rank k for C^k, given kappa_a and kappa_b.

Returns \( (k_\text{min}, k_\text{max}) \) such that \( \langle \kappa_a \| C^k \| \kappa_b \rangle \neq 0 \) only for \( k_\text{min} \leq k \leq k_\text{max} \). Steps are in increments of 2, enforced by parity.

Derived from:

\[ k_\text{min} = |j_a - j_b|, \quad k_\text{max} = j_a + j_b, \]

then adjusted so that \( l_a + l_b + k \) is even.

Parameters

ka	Dirac quantum number \( \kappa_a \).
kb	Dirac quantum number \( \kappa_b \).

Returns

\( \{k_\text{min},\, k_\text{max}\} \)

S_kk()

double Angular::S_kk ( int  ka, int  kb  )

inline

Reduced spin-1/2 angular ME: <ka||S||kb>. Takes kappa values.

Computes the reduced matrix element of the spin operator \( \mathbf{S} \) (for spin-1/2) between relativistic states:

\[ \langle \kappa_a \| S \| \kappa_b \rangle = \delta_{l_a l_b} \, (-1)^{j_a+l_a+3/2} \sqrt{[j_a][j_b]\tfrac{3}{2}} \sixj{j_a}{1}{j_b}{\tfrac{1}{2}}{l_a}{\tfrac{1}{2}} \]

The 6j symbol is evaluated analytically:

\[ \sixj{j_a}{1}{j_b}{\tfrac{1}{2}}{l_a}{\tfrac{1}{2}} = \frac{(-1)^{\kappa_a+\kappa_b}}{2} \sqrt{\left|\frac{(\kappa_a-1)^2 - \kappa_b^2}{3\kappa_a(1+2\kappa_a)}\right|} \]

This special-case formula is ~20% faster than the general 6j routine. Returns 0 if \( l_a \neq l_b \) or the triangle rule is violated.

Parameters

ka	Dirac quantum number \( \kappa_a \).
kb	Dirac quantum number \( \kappa_b \).

Returns

\( \langle \kappa_a \| S \| \kappa_b \rangle \); 0 if selection rules are not met.