-
Notifications
You must be signed in to change notification settings - Fork 15
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
First commit of working code. Needs cleanup.
- Loading branch information
1 parent
b235ba3
commit b262e89
Showing
7 changed files
with
1,788 additions
and
0 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Original file line number | Diff line number | Diff line change |
---|---|---|
@@ -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; | ||
} | ||
|
Oops, something went wrong.