New version for issue #2.

joeydumont · Jul 19, 2014 · a8cc337 · a8cc337
1 parent 4c61b93
commit a8cc337
Show file tree

Hide file tree

Showing 3 changed files with 22 additions and 11 deletions.
diff --git a/src/wignerSymbols-cpp.cpp b/src/wignerSymbols-cpp.cpp
@@ -13,7 +13,7 @@ std::vector<double> wigner3j(double l2, double l3,
  double m1, double m2, double m3)
 {
  // We compute the numeric limits of double precision. 
- double huge = std::numeric_limits<double>::max();
+ double huge = sqrt(std::numeric_limits<double>::max()/20.0);
  double srhuge = sqrt(huge);
  double tiny = std::numeric_limits<double>::min();
  double srtiny = sqrt(tiny);
@@ -103,8 +103,9 @@ std::vector<double> wigner3j(double l2, double l3,
  std::cout << "We renormalized the forward recursion." << std::endl;
  for (std::vector<double>::iterator it = thrcof.begin(); it != thrcof.begin()+i; ++it)
  {
- if (std::fabs(*it) < srtiny) *it = 0;
- else *it /= srhuge;
+ //if (std::fabs(*it) < srtiny) *it = 0;
+ //else 
+ *it /= srhuge;
  }
  }
 
@@ -160,8 +161,9 @@ std::vector<double> wigner3j(double l2, double l3,
  std::cout << "We renormalized the backward recursion." << std::endl;
  for (std::vector<double>::iterator it = thrcof.begin()+j; it != thrcof.end(); ++it)
  {
- if (std::fabs(*it) < srtiny) *it = 0;
- else *it /= srhuge;
+ //if (std::fabs(*it) < srtiny) *it = 0;
+ //else 
+ *it /= srhuge;
  }
  }
 
@@ -186,12 +188,16 @@ std::vector<double> wigner3j(double l2, double l3,
  {
  sum += (2.0*(l1min+k)+1.0)*thrcof[k]*thrcof[k];
  }
+ std::cout << sum << std::endl;
 
- double c1 = pow(-1.0,l2-l3-m1)*sgn(thrcof[size-1])/sqrt(sum);
-
+ std::cout << "(-1)^(l2-l3-m1): " << pow(-1.0,l2-l3-m1) << " sgn:" << sgn(thrcof[size-1]) << std::endl;
+ double c1 = pow(-1.0,l2-l3-m1)*sgn(thrcof[size-1]);
+ std::cout << "c1: " << c1 << std::endl;
  for (std::vector<double>::iterator it = thrcof.begin(); it != thrcof.end(); ++it)
  {
- *it *= c1;
+ std::cout << *it << ", " << c1 << ", ";
+ *it *= c1/sqrt(sum);
+ std::cout << *it << std::endl;
  }
  return thrcof;
 }

diff --git a/tests/gh-issue-1.cpp b/tests/gh-issue-1.cpp
@@ -37,7 +37,6 @@ int main (int argc, char* argv[])
  double test31 = WignerSymbols::wigner3j(772, 874, 1231, 442, 756, -1198);
  double test32 = WignerSymbols::wigner3j(874, 1231, 772, 756, -1198, 442);
  double test33 = WignerSymbols::wigner3j(1231, 772, 874, -1198, 442, 756);
-
  double test34 = WignerSymbols::wigner3j_f(772, 874, 1231, 442, 756, -1198);
  double test35 = WignerSymbols::wigner3j_f(874, 1231, 772, 756, -1198, 442);
  double test36 = WignerSymbols::wigner3j_f(1231, 772, 874, -1198, 442, 756);
@@ -59,6 +58,12 @@ int main (int argc, char* argv[])
  std::cout << "FOR impl.: " << test24 << ", " << test25 << ", " << test26 << std::endl;
  std::cout << "C++ impl.: " << test31 << ", " << test32 << ", " << test33 << std::endl;
  std::cout << "FOR impl.: " << test34 << ", " << test35 << ", " << test36 << std::endl;
+
+ std::vector<double> test4 = WignerSymbols::wigner3j(992, 1243, 196, -901, 705);
+ double test41 = WignerSymbols::wigner3j(529, 992, 1243, 196, -901, 705);
+ double test42 = WignerSymbols::wigner3j_f(529, 992, 1243, 196, -901, 705);
+ std::cout << "C++ impl.: " << test41 << std::endl;
+ std::cout << "FOR impl.: " << test42 << std::endl;
  std::cout << std::endl;
 
  return 0;

diff --git a/tests/testWigner.cpp b/tests/testWigner.cpp
@@ -134,8 +134,8 @@ double firstOrthoRelation3j(double l1, double l2, double l3, double m3)
  {
  for (int j=0;j<m2Size;j++)
  {
- sum += (2.0*l3+1.0)*WignerSymbols::wigner3j(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3)
- *WignerSymbols::wigner3j(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3);
+ sum += (2.0*l3+1.0)*WignerSymbols::wigner3j_f(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3)
+ *WignerSymbols::wigner3j_f(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3);
  }
  }