qpms/wigner.c

   1 #include "quaternions.h"
   2 #include "tiny_inlines.h"
   3 #include "kahansum.h"
   4 #define WIGNER_ATOL (1e-15)
   5
   6 complex double qpms_wignerD_elem(const qpms_quat_t R,
   7     const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m) {
   8   // TODO do some optimisations... The combinatoric coeffs could be precomputed.
   9   QPMS_ENSURE(abs(m) <= l || abs(mp) <= l, "Got invalid values l = %d, m = %d", (int)l, (int)m);
  10   const double mRa = cabs(R.a), mRb = cabs(R.b);
  11   if (mRa < WIGNER_ATOL) {
  12     if (m != -mp) return 0;
  13     else return min1pow(l+m) * cpow(R.b, 2*m);
  14   } else if (mRb < WIGNER_ATOL) {
  15     if (m != mp) return 0;
  16     else return cpow(R.a, 2*m);
  17   } else if (mRa >= mRb) {
  18     complex double prefac = exp(.5*(lgamma(l+m+1) + lgamma(l-m+1)
  19         - lgamma(l+mp+1) - lgamma(l-mp+1)))
  20       * pow(mRa, 2*l-2*m)
  21       * cpow(R.a, mp+m)
  22       * cpow(R.b, -mp+m);
  23     double sum, sumc; kahaninit( &sum, &sumc);
  24     // FIXME this is probably quite stupid way to calculate the sum
  25     for(int rho=MAX(0, -m+mp); rho <= MIN(l-m, l+mp); ++rho)
  26       kahanadd(&sum, &sumc,
  27           exp(
  28             lgamma(l+mp+1) - lgamma(rho+1) - lgamma(l+mp-rho+1)
  29             + lgamma(l-mp+1) - lgamma(l-rho-m+1) - lgamma(-mp+rho+m+1)
  30           )
  31           * pow(-mRb*mRb/(mRa*mRa), rho)
  32       );
  33     return prefac * sum;
  34   } else { // (mRa < mRb)
  35     complex double prefac = min1pow(l-m) * exp(.5*(lgamma(l+m+1) + lgamma(l-m+1)
  36         - lgamma(l+mp+1) - lgamma(l-mp+1)))
  37       * pow(mRb, 2*l-2*m)
  38       * cpow(R.a, mp+m)
  39       * cpow(R.b, -mp+m);
  40     double sum, sumc; kahaninit( &sum, &sumc);
  41     // FIXME this is probably quite stupid way to calculate the sum
  42     for(int rho=MAX(0, -m-mp); rho <= MIN(l-mp, l-m); ++rho)
  43       kahanadd(&sum, &sumc,
  44           exp(
  45             lgamma(l+mp+1) - lgamma(l-m-rho+1) - lgamma(mp+m+rho+1)
  46             + lgamma(l-mp+1) - lgamma(rho+1) - lgamma(l-mp-rho+1)
  47           )
  48           * pow(-mRa*mRa/(mRb*mRb), rho)
  49       );
  50     return prefac * sum;
  51   }
  52 }
  53
  54 complex double qpms_vswf_irot_elem_from_irot3(const qpms_irot3_t q,
  55     const qpms_l_t l, const qpms_m_t mp, const qpms_m_t m, bool pseudo) {
  56 #ifndef NDEBUG
  57   qpms_irot3_checkdet(q);
  58 #endif
  59   complex double res = qpms_wignerD_elem(q.rot, l, mp, m);
  60   if (q.det == -1) {
  61     res *= min1pow(l);
  62     if (pseudo)
  63       res *= -1;
  64   }
  65   return res;
  66 }