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