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
  11   double cgc (int factor1_2j, int factor1_2m, int factor2_2j, int factor2_2m,
  12               int irrep_2j, int irrep_2m)
  13   {
  14     const int factor1_n = (factor1_2j + factor1_2m) / 2;
  15     const int factor2_n = (factor2_2j + factor2_2m) / 2;
  16     const int irrep_n =   (irrep_2j + irrep_2m) / 2;
  17
  18     clebsch::weight factor1_S(2);
  19     factor1_S(1) = factor1_2j;
  20     if ( ! (0 <= factor1_n && factor1_n < factor1_S.dimension()) ) {
  21       //std::cerr << "clebsch warning: 0 <= factor1_state && factor1_state < factor1_dimension" << std::endl;
  22       return 0.0;
  23     }
  24
  25     clebsch::weight factor2_S(2);
  26     factor2_S(1) = factor2_2j;
  27     if ( ! (0 <= factor2_n && factor2_n < factor2_S.dimension()) ) {
  28       //std::cerr << "clebsch warning: 0 <= factor2_state && factor2_state < factor2_dimension" << std::endl;
  29       return 0.0;
  30     }
  31
  32     clebsch::weight irrep_S(2);
  33     irrep_S(1) = irrep_2j;
  34     if ( ! (0 <= irrep_n && irrep_n < irrep_S.dimension()) ) {
  35       //std::cerr << "clebsch warning: 0 <= irrep_state && irrep_state < irrep_dimension" << std::endl;
  36       return 0.0;
  37     }
  38
  39     clebsch::decomposition decomp(factor1_S, factor2_S);
  40
  41     bool b = true;
  42     for (int i = 0; i < decomp.size(); ++i) {
  43       if (decomp(i) == irrep_S) {
  44         b = false;
  45         break;
  46       }
  47     }
  48     if (b) {
  49       cerr << "WARNING: 2j=" << irrep_2j
  50            << " does not occour in decomposition!" << endl;
  51       return 0.0;
  52     }
  53
  54     const clebsch::coefficients C(irrep_S, factor1_S, factor2_S);
  55
  56     return C( factor1_n, factor2_n, 0, irrep_n );
  57   }
  58
  59 };
  60