CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

ClebschGordanCoefficientSet.cc
Go to the documentation of this file.
2 #include <iostream>
3 #include <cmath>
4 #include <stdexcept>
5 inline double factorial (int N) {
6  double retVal(1.0);
7  for (int i=2;i<=N;i++) retVal*=i;
8  return retVal;
9 }
10 
11 
12 namespace Genfun {
13  double ClebschGordanCoefficientSet::calcCoefficient(int l1, int l2, int L, int xm1, int xm2, int M) {
14 
15  if (xm1+xm2!=M) return 0;
16  double F1=sqrt((2*L+1)*factorial(L+l1-l2)*factorial(L-l1+l2)*factorial(l1+l2-L)/factorial(l1+l2+L+1));
17  double F2=sqrt(factorial(L+M)*factorial(L-M)*factorial(l1-xm1)*factorial(l1+xm1)*factorial(l2-xm2)
18  *factorial(l2+xm2));
19  double F3=0.0;
20  int max=0;
21  max = std::max(max,l2+xm2);
22  max = std::max(max,l1-xm1);
23  max = std::max(max,l1+l2-L);
24  // max = std::max(max,l2-L-xm1);
25  // max = std::max(max,l1+xm2-L);
26  for (int k=0;k<=max;k++) {
27  double term = factorial(k);
28  bool skipTerm=false;
29  {
30  int T=l1 + l2 -L -k;
31  if (T>=0) term *= factorial(T); else skipTerm=true;
32  }
33  if (!skipTerm)
34  {
35  int T=l1-xm1-k;
36  if (T>=0) term *= factorial(T); else skipTerm=true;
37  }
38  if (!skipTerm)
39  {
40  int T=l2+xm2-k;
41  if (T>=0) term *= factorial(T); else skipTerm=true;
42  }
43  if (!skipTerm)
44  {
45  int T=L-l2+xm1+k;
46  if (T>=0) term *= factorial(T); else skipTerm=true;
47  }
48  if (!skipTerm)
49  {
50  int T=L-l1-xm2+k;
51  if (T>=0) term *= factorial(T); else skipTerm=true;
52  }
53  if (!skipTerm) F3+= ((k%2) ? -1:1)/term;
54  }
55 
56 
57  return F1*F2*F3;
58  }
59 }