su2clebsch.cpp

   1 #include <iostream>
   2
   3 #include "clebsch.h"
   4
   5 #include "su2clebsch.hpp"
   6
   7 using namespace std;
   8
   9 namespace su2clebsch {
  10   int ipow(int base, int exp)
  11   {
  12     int result = 1;
  13     while (exp)
  14       {
  15         if (exp & 1)
  16           result *= base;
  17         exp >>= 1;
  18         base *= base;
  19       }
  20
  21     return result;
  22   }
  23
  24   prestore::~prestore ()
  25   {
  26     delete[] cgcdata;
  27   }
  28
  29   prestore::prestore(unsigned int _max2j)
  30   {
  31     max2j = _max2j;
  32     cgcdata = new double[ipow(max2j+1,3)*ipow(2*max2j+1,3)];
  33
  34     for (int factor1_2j=0; factor1_2j<=max2j; factor1_2j++)
  35       for (int factor1_mp=0; factor1_mp<=2*max2j; factor1_mp++)
  36         for (int factor2_2j=0; factor2_2j<=max2j; factor2_2j++)
  37           for (int factor2_mp=0; factor2_mp<=2*max2j; factor2_mp++)
  38             for (int irrep_2j=0; irrep_2j<=max2j; irrep_2j++)
  39               for (int irrep_mp=0; irrep_mp<=2*max2j; irrep_mp++)
  40                 {
  41                   cgcdata[ factor1_2j*ipow(max2j+1,2)*ipow(2*max2j+1,3) +
  42                            factor1_mp*ipow(max2j+1,2)*ipow(2*max2j+1,2) +
  43                            factor2_2j*(max2j+1)*ipow(2*max2j+1,2) +
  44                            factor2_mp*(max2j+1)*(2*max2j+1) +
  45                            irrep_2j*(2*max2j+1) +
  46                            irrep_mp ]
  47                     =
  48                     cgc( factor1_2j, factor1_mp-max2j, factor2_2j, factor2_mp-max2j,
  49                          irrep_2j, irrep_mp-max2j );
  50                 }
  51   }
  52
  53   bool stateValid(int state_2j, int state_2m)
  54   {
  55     if ( abs(state_2m) > state_2j || (state_2j + state_2m) % 2 != 0 )
  56       return false;
  57     else
  58       return true;
  59   }
  60
  61   double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
  62               int irrep_2j, int irrep_2m, const prestore& Prestore)
  63   {
  64     if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) )
  65       return 0.0;
  66
  67     if ( factor1_2j > Prestore.max2j || factor2_2j > Prestore.max2j || irrep_2j > Prestore.max2j )
  68       return cgc(factor1_2j, factor1_2m, factor2_2j, factor2_2m, irrep_2j, irrep_2m);
  69
  70     return Prestore.cgcdata[ factor1_2j*ipow(Prestore.max2j+1,2)*ipow(2*Prestore.max2j+1,3) +
  71                              (factor1_2m+Prestore.max2j) * ipow(Prestore.max2j+1,2)*ipow(2*Prestore.max2j+1,2) +
  72                              factor2_2j*(Prestore.max2j+1)*ipow(2*Prestore.max2j+1,2) +
  73                              (factor2_2m+Prestore.max2j) * (Prestore.max2j+1)*(2*Prestore.max2j+1) +
  74                              irrep_2j*(2*Prestore.max2j+1) +
  75                              (irrep_2m+Prestore.max2j) ];
  76   }
  77
  78   double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
  79               int irrep_2j, int irrep_2m)
  80   {
  81     if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) )
  82       return 0.0;
  83
  84     const int factor1_n = (factor1_2j + factor1_2m) / 2;
  85     const int factor2_n = (factor2_2j + factor2_2m) / 2;
  86     const int irrep_n =   (irrep_2j + irrep_2m) / 2;
  87
  88     clebsch::weight factor1_S(2);
  89     factor1_S(1) = factor1_2j;
  90     clebsch::weight factor2_S(2);
  91     factor2_S(1) = factor2_2j;
  92     clebsch::weight irrep_S(2);
  93     irrep_S(1) = irrep_2j;
  94
  95     clebsch::decomposition decomp(factor1_S, factor2_S);
  96
  97     bool b = true;
  98     for (int i = 0; i < decomp.size(); ++i) {
  99       if (decomp(i) == irrep_S) {
 100         b = false;
 101         break;
 102       }
 103     }
 104     if (b) {
 105       return 0.0;
 106     }
 107
 108     const clebsch::coefficients C(irrep_S, factor1_S, factor2_S);
 109
 110     return C( factor1_n, factor2_n, 0, irrep_n );
 111   }
 112
 113 };
 114