-
Notifications
You must be signed in to change notification settings - Fork 15
/
wignerSymbols-fortran.cpp
86 lines (70 loc) · 2.5 KB
/
wignerSymbols-fortran.cpp
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
/*******************************************************-/
* This source code is subject to the terms of the GNU -/
* Lesser Public License. If a copy of the LGPL was not -/
* distributed with this file, you can obtain one at -/
* https://www.gnu.org/licenses/lgpl.html. -/
********************************************************/
#include "../include/wignerSymbols/wignerSymbols-fortran.h"
namespace WignerSymbols {
/*! Computes the Wigner-3j symbol for given l1,l2,l3,m1,m2,m3. We
* explicitly enforce the selection rules. */
double wigner3j_f(double l1, double l2, double l3,
double m1, double m2, double m3)
{
// We enforce the selection rules.
bool select(true);
select = (
m1+m2+m3==0.0
&& std::floor(l1+l2+l3)==(l1+l2+l3)
&& l3 >= std::fabs(l1-l2)
&& l3 <= l1+l2
&& std::fabs(m1) <= l1
&& std::fabs(m2) <= l2
&& std::fabs(m3) <= l3
);
if (!select) return 0.0;
// We compute the size of the resulting array.
int size = (int)std::ceil(l2+l3-std::max(std::fabs(l2-l3),std::fabs(m1)))+1;
// We prepare the output values.
double l1min, l1max;
double thrcof [size];
int ierr;
// External function call.
drc3jj_wrap(l2,l3,m2,m3,&l1min,&l1max,thrcof,size,&ierr);
// We fetch and return the value with the proper l1 value.
int index = (int)l1-l1min;
return thrcof[index];
}
double wigner6j_f(double l1, double l2, double l3,
double l4, double l5, double l6)
{
// We enforce the selection rules.
bool select(true);
// Triangle relations for the four tryads
select = (
std::fabs(l1-l2) <= l3 && l3 <= l1+l2
&& std::fabs(l1-l5) <= l6 && l6 <= l1+l5
&& std::fabs(l4-l2) <= l6 && l6 <= l4+l2
&& std::fabs(l4-l5) <= l3 && l3 <= l4+l5
);
// Sum rule of the tryads
select = (
std::floor(l1+l2+l3)==(l1+l2+l3)
&& std::floor(l1+l5+l6)==(l1+l5+l6)
&& std::floor(l4+l2+l6)==(l4+l2+l6)
&& std::floor(l4+l5+l3)==(l4+l5+l3)
);
if (!select) return 0.0;
// We compute the size of the resulting array.
int size = (int)std::ceil(std::min(l2+l3,l5+l6)-std::max(std::fabs(l2-l3),std::fabs(l5-l6)))+1;
// We prepare the output values
double l1min, l1max;
double sixcof [size];
int ierr;
// External function call
drc6j_wrap(l2,l3,l4,l5,l6,&l1min,&l1max,sixcof,size,&ierr);
// We fetch and return the coefficient with the proper l1 value.
int index = (int)l1-l1min;
return sixcof[index];
}
}