67 lines
2.3 KiB
C
67 lines
2.3 KiB
C
#include "quaternions.h"
|
|
#include "tiny_inlines.h"
|
|
#include "kahansum.h"
|
|
#define WIGNER_ATOL (1e-15)
|
|
|
|
complex double qpms_wignerD_elem(const qpms_quat_t R,
|
|
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m) {
|
|
// TODO do some optimisations... The combinatoric coeffs could be precomputed.
|
|
QPMS_ENSURE(abs(m) <= l || abs(mp) <= l, "Got invalid values l = %d, m = %d", (int)l, (int)m);
|
|
const double mRa = cabs(R.a), mRb = cabs(R.b);
|
|
if (mRa < WIGNER_ATOL) {
|
|
if (m != -mp) return 0;
|
|
else return min1pow(l+m) * cpow(R.b, 2*m);
|
|
} else if (mRb < WIGNER_ATOL) {
|
|
if (m != mp) return 0;
|
|
else return cpow(R.a, 2*m);
|
|
} else if (mRa >= mRb) {
|
|
complex double prefac = exp(.5*(lgamma(l+m+1) + lgamma(l-m+1)
|
|
- lgamma(l+mp+1) - lgamma(l-mp+1)))
|
|
* pow(mRa, 2*l-2*m)
|
|
* cpow(R.a, mp+m)
|
|
* cpow(R.b, -mp+m);
|
|
double sum, sumc; kahaninit( &sum, &sumc);
|
|
// FIXME this is probably quite stupid way to calculate the sum
|
|
for(int rho=MAX(0, -m+mp); rho <= MIN(l-m, l+mp); ++rho)
|
|
kahanadd(&sum, &sumc,
|
|
exp(
|
|
lgamma(l+mp+1) - lgamma(rho+1) - lgamma(l+mp-rho+1)
|
|
+ lgamma(l-mp+1) - lgamma(l-rho-m+1) - lgamma(-mp+rho+m+1)
|
|
)
|
|
* pow(-mRb*mRb/(mRa*mRa), rho)
|
|
);
|
|
return prefac * sum;
|
|
} else { // (mRa < mRb)
|
|
complex double prefac = min1pow(l-m) * exp(.5*(lgamma(l+m+1) + lgamma(l-m+1)
|
|
- lgamma(l+mp+1) - lgamma(l-mp+1)))
|
|
* pow(mRb, 2*l-2*m)
|
|
* cpow(R.a, mp+m)
|
|
* cpow(R.b, -mp+m);
|
|
double sum, sumc; kahaninit( &sum, &sumc);
|
|
// FIXME this is probably quite stupid way to calculate the sum
|
|
for(int rho=MAX(0, -m-mp); rho <= MIN(l-mp, l-m); ++rho)
|
|
kahanadd(&sum, &sumc,
|
|
exp(
|
|
lgamma(l+mp+1) - lgamma(l-m-rho+1) - lgamma(mp+m+rho+1)
|
|
+ lgamma(l-mp+1) - lgamma(rho+1) - lgamma(l-mp-rho+1)
|
|
)
|
|
* pow(-mRa*mRa/(mRb*mRb), rho)
|
|
);
|
|
return prefac * sum;
|
|
}
|
|
}
|
|
|
|
complex double qpms_vswf_irot_elem_from_irot3(const qpms_irot3_t q,
|
|
const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, bool pseudo) {
|
|
#ifndef NDEBUG
|
|
qpms_irot3_checkdet(q);
|
|
#endif
|
|
complex double res = qpms_wignerD_elem(q.rot, l, mp, m);
|
|
if (q.det == -1) {
|
|
res *= min1pow(l);
|
|
if (pseudo)
|
|
res *= -1;
|
|
}
|
|
return res;
|
|
}
|