Added Wigner symbols source files.

src/fortranLinkage.h

@@ -0,0 +1,10 @@
+#ifndef FORTRANLINKAGE_H
+#define FORTRANLINKAGE_H
+
+extern "C"
+{
+ extern void drc3jj_wrap(double,double,double,double,double*,double*,double*,int,int*);
+ extern void drc6j_wrap(double,double,double,double,double*,double*,double*,int,int*);
+}
+
+#endif // FORTRANLINKAGE_H

src/testWigner.cpp

@@ -0,0 +1,220 @@
+#include "wignerSymbols.h"
+
+#include <iostream>
+#include <iomanip>
+#include <armadillo>
+#include <time.h>
+
+#define EPS 1.0e-10
+
+using namespace std;
+using namespace arma;
+
+// Orthogonality relations and sum rules for the 3j symbols [Brink & Satchler, pp. 136+139].
+double firstOrthoRelation3j(double l1, double l2, double l3, double m3);
+double sumOverM23j(double l1, double l2, double l3, double m3);
+
+// Orthogonality relations for the 6j symbols [Brink & Satchler pp. 142-3].
+double firstSumOverL3(double l1, double l2, double l6);
+double seconSumOverL3(double l1, double l2, double l6);
+
+// Timing computation of Wigner symbols. 
+double timingWignerSymbols(double l2, double l3, double m1, double m2, double m3);
+double timingWignerSymbolsF(double l2, double l3, double m1, double m2, double m3);
+
+int main(int argc, char* argv[])
+{
+ // We time the computation for different sizes.
+ int N = 500;
+ double l2(500.0),m1(-10.0),m2(60.0),m3(-50.0);
+ mat times(N,2);
+ mat timesF(N,2);
+
+ for (int i=0;i<N;i++)
+ {
+ int l3 = i+std::fabs(m3);
+ times(i,0) = (int)std::ceil(l2+l3-std::max(std::fabs(l2-l3),std::fabs(m1)))+1;
+ timesF(i,0) = times(i,0);
+ times(i,1) = timingWignerSymbols(l2,l3,m1,m2,m3);
+ timesF(i,1) = timingWignerSymbolsF(l2,l3,m1,m2,m3);
+ }
+
+ times.save("times.dat", raw_ascii);
+ timesF.save("timesF.dat",raw_ascii);
+
+ // Generate values for l1, l2 and derive the rest of the allowed values.
+ double lMax = atof(argv[1]);
+ colvec lVec = linspace(0,lMax,lMax+1);
+
+ // Run over l1
+ for (int i=0;i<lMax+1;i++)
+ {
+ double l1 = lVec(i);
+
+ // Run over l2
+ for (int j=0;j<lMax+1;j++)
+ {
+
+ double l2 = lVec(j);
+ double l3min = std::fabs(l1-l2);
+ double l3max = l1+l2;
+
+ // Run over l3
+ for (int k=0;k<(l3max-l3min+1);k++)
+ {
+ double l3 = l3min+k;
+
+ double test = firstSumOverL3(l1,l2,l3); 
+
+ if (!test)
+ cout << "The first sum over L3 test failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl6 = " << l3 << endl;
+
+ // Run over m3
+ mat precisionWigner(2*l3+1,2);
+ for (int l=0;l<(2*l3+1);l++)
+ {
+ double m3 = -l3+l;
+ test = firstOrthoRelation3j(l1,l2,l3,m3);
+
+ precisionWigner(l,0) = m3;
+ precisionWigner(l,1) = test;
+
+ if (!test)
+ cout << "The first orthogonality test failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl3 = " << l3 << "\n\tm3 = " << m3 << endl;
+
+ //if (l3!=0.0) test = sumOverM23j(l1,l2,l3,m3);
+
+ if (!test)
+ cout << "The sum over M2 test failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl3 = " << l3 << "\n\tm3 = " << m3 << endl;
+ }
+ precisionWigner.save("precision.dat", raw_ascii);
+ }
+
+ double l6max = std::min(2.*l1,2.*l2);
+
+ // Run over l6
+ for (int k=0;k<l6max+1;k++)
+ {
+ bool test = seconSumOverL3(l1,l2,k);
+
+ if (!test)
+ cout << "The second sum over L3 test has failed for \n\tl1 = " << l1 << "\n\tl2 = " << l2 << "\n\tl6 = " << k << endl;
+ }
+ }
+ }
+
+ return 0;
+}
+
+// Sum_{m1,m2} (2*l3+1)*Wigner3j(l1,l2,l3,m1,m2,m3)*3j(l1,l2,l3',m1,m2,m3') = delta(l3,l3')*delta(m3,m3')
+double firstOrthoRelation3j(double l1, double l2, double l3, double m3)
+{
+ // We determine all possible values of m1 and m2.
+ int m1Size = ((int)(2.0*l1+1.0));
+ int m2Size = ((int)(2.0*l2+1.0)); 
+
+ double sum = 0.0;
+ for (int i=0;i<m1Size;i++)
+ {
+ for (int j=0;j<m2Size;j++)
+ {
+ sum += (2.0*l3+1.0)*varPhase::c_wigner3j(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3)
+ *varPhase::c_wigner3j(l1,l2,l3,((double)(-l1+i)),((double)(-l2+j)),m3);
+ }
+ }
+
+ double diff = std::fabs(1.0-sum);
+
+ return diff;
+}
+
+// sum_{m1,m2} m2*(2*l3+1)c_wigner3j(l1,l2,l3,m1,m2,-m3)^2 = m3*(l3*(l3+1)+l2*(l2+1)-l1*(l1+1))/(2*l3*(l3+1))
+double sumOverM23j(double l1, double l2, double l3, double m3)
+{
+ double sum = 0.0;
+ // We sum over all allowable values of m1+m2+m3=0.
+ for (int i=0;i<(2*l1+1);i++)
+ {
+ double m1 = -l1+i;
+
+ for (int j=0;j<(2*l2+1);j++)
+ {
+ double m2 = -l2+j;
+
+ if (std::fabs(m1+m2-m3)<1.0e-6)
+ sum += m2*(2.0*l3+1.0)*pow(varPhase::c_wigner3j(l1,l2,l3,m1,m2,-m3),2.0);
+ }
+ }
+
+ double value = m3*(l3*(l3+1.)+l2*(l2+1.)-l1*(l1+1.))/(2.*l3*(l3+1.));
+
+ return std::fabs(sum-value);
+}
+
+// sum_{l3} = (-1)^(2*l3)*(2*l3+1)*Wigner6j(l1,l2,l3,l1,l2,l6) = 1
+double firstSumOverL3(double l1, double l2, double l6)
+{
+ // We sum over the value of l3. 
+ double l3min = std::fabs(l2-l1);
+ double l3max = l1+l2;
+
+ double sum = 0.0;
+ for (int i=0;i<l3max-l3min+1;i++)
+ {
+ double l3 = l3min+i;
+ sum += pow(-1.,2.0*l3)*(2.0*l3+1.0)*varPhase::c_wigner6j(l1,l2,l3,l1,l2,l6);
+ }
+
+ double diff = std::fabs(1.0-sum);
+
+ return diff;
+}
+
+// sum_{l3} = (-1)^(l1+l2+l3)*(2*l3+1)*Wigner6j(l1,l2,l3,l2,l1,l6) = delta(l6,0)*sqrt((2*l1+1)*(2*l2+1))
+double seconSumOverL3(double l1, double l2, double l6)
+{
+ // We sum over the value of l3. 
+ double l3min = std::fabs(l1-l2);
+ double l3max = l1+l2;
+
+ double sum = 0.0;
+ for (int i=0;i<l3max-l3min+1;i++)
+ {
+ double l3 = l3min+i;
+ sum += pow(-1.,l1+l2+l3)*(2.0*l3+1.0)*varPhase::c_wigner6j(l1,l2,l3,l2,l1,l6);
+ }
+
+ double value = (l6==0.0 ? sqrt((2.*l1+1.)*(2.*l2+1.)) : 0.0) ;
+
+ double diff = std::fabs(value-sum);
+
+ return std::fabs(value-sum);
+
+}
+
+double timingWignerSymbols(double l2, double l3, double m1, double m2, double m3)
+{
+ // We start the timer. 
+ clock_t start = clock();
+
+ // We call the subroutine.
+ colvec test = varPhase::c_wigner3j(l2,l3,m1,m2,m3);
+
+ clock_t end = clock();
+
+ return ((double)end-start)/CLOCKS_PER_SEC;
+}
+
+double timingWignerSymbolsF(double l2, double l3, double m1, double m2, double m3)
+{
+ // We start the timer. 
+ clock_t start = clock();
+
+ // We call the subroutine.
+ double test = varPhase::wigner3j(l2+l3,l2,l3,m1,m2,m3);
+
+ clock_t end = clock();
+
+ return ((double)end-start)/CLOCKS_PER_SEC;
+}
+