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,6)];
  33
  34     for (int factor1_2j=0; factor1_2j<=max2j; factor1_2j++)
  35       for (int factor1_mp=0; factor1_mp<=max2j; factor1_mp++)
  36         for (int factor2_2j=0; factor2_2j<=max2j; factor2_2j++)
  37           for (int factor2_mp=0; factor2_mp<=max2j; factor2_mp++)
  38             for (int irrep_2j=0; irrep_2j<=max2j; irrep_2j++)
  39               for (int irrep_mp=0; irrep_mp<=max2j; irrep_mp++)
  40                 {
  41                   cgcdata[ factor1_2j*ipow(max2j+1,5) + factor1_mp*ipow(max2j+1,4) +
  42                            factor2_2j*ipow(max2j+1,3) + factor2_mp*ipow(max2j+1,2) +
  43                            irrep_2j*(max2j+1)         + irrep_mp ]
  44                     =
  45                     cgc( factor1_2j, 2*factor1_mp-max2j, factor2_2j, 2*factor2_mp-max2j,
  46                          irrep_2j, 2*irrep_mp-max2j );
  47                 }
  48   }
  49
  50   bool stateValid(int state_2j, int state_2m)
  51   {
  52     if ( abs(state_2m) > state_2j || (state_2j + state_2m) % 2 != 0 )
  53       return false;
  54     else
  55       return true;
  56   }
  57
  58   double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
  59               int irrep_2j, int irrep_2m, const prestore& Prestore)
  60   {
  61     if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) )
  62       return 0.0;
  63
  64     if ( factor1_2j > Prestore.max2j || factor2_2j > Prestore.max2j || irrep_2j > Prestore.max2j )
  65       return cgc(factor1_2j, factor1_2m, factor2_2j, factor2_2m, irrep_2j, irrep_2m);
  66
  67     return Prestore.cgcdata[ factor1_2j*ipow(Prestore.max2j+1,5) +
  68                              (factor1_2m+Prestore.max2j)/2 *ipow(Prestore.max2j+1,4) +
  69                              factor2_2j*ipow(Prestore.max2j+1,3) +
  70                              (factor2_2m+Prestore.max2j)/2 *ipow(Prestore.max2j+1,2) +
  71                              irrep_2j*(Prestore.max2j+1) +
  72                              (irrep_2m+Prestore.max2j)/2 ];
  73   }
  74
  75   double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
  76               int irrep_2j, int irrep_2m)
  77   {
  78     if ( ! ( stateValid(factor1_2j, factor1_2m) && stateValid(factor2_2j, factor2_2m) && stateValid(irrep_2j, irrep_2m) ) )
  79       return 0.0;
  80
  81     const int factor1_n = (factor1_2j + factor1_2m) / 2;
  82     const int factor2_n = (factor2_2j + factor2_2m) / 2;
  83     const int irrep_n =   (irrep_2j + irrep_2m) / 2;
  84
  85     clebsch::weight factor1_S(2);
  86     factor1_S(1) = factor1_2j;
  87     clebsch::weight factor2_S(2);
  88     factor2_S(1) = factor2_2j;
  89     clebsch::weight irrep_S(2);
  90     irrep_S(1) = irrep_2j;
  91
  92     clebsch::decomposition decomp(factor1_S, factor2_S);
  93
  94     bool b = true;
  95     for (int i = 0; i < decomp.size(); ++i) {
  96       if (decomp(i) == irrep_S) {
  97         b = false;
  98         break;
  99       }
 100     }
 101     if (b) {
 102       return 0.0;
 103     }
 104
 105     const clebsch::coefficients C(irrep_S, factor1_S, factor2_S);
 106
 107     return C( factor1_n, factor2_n, 0, irrep_n );
 108   }
 109
 110 };
 111