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 m1, int m2, int M) {
14 
15  if (m1+m2!=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-m1)*factorial(l1+m1)*factorial(l2-m2)
18  *factorial(l2+m2));
19  double F3=0.0;
20  int max=0;
21  max = std::max(max,l2+m2);
22  max = std::max(max,l1-m1);
23  max = std::max(max,l1+l2-L);
24  // max = std::max(max,l2-L-m1);
25  // max = std::max(max,l1+m2-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-m1-k;
36  if (T>=0) term *= factorial(T); else skipTerm=true;
37  }
38  if (!skipTerm)
39  {
40  int T=l2+m2-k;
41  if (T>=0) term *= factorial(T); else skipTerm=true;
42  }
43  if (!skipTerm)
44  {
45  int T=L-l2+m1+k;
46  if (T>=0) term *= factorial(T); else skipTerm=true;
47  }
48  if (!skipTerm)
49  {
50  int T=L-l1-m2+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 }