ARTS 2.5.11 (git: 725533f0)
wigner_functions.cc
Go to the documentation of this file.
1
9#include "wigner_functions.h"
10
11#include <sys/errno.h>
12
13#include <algorithm>
14
15#include "arts_omp.h"
16#include "arts_conversions.h"
17#include "debug.h"
18
19#if DO_FAST_WIGNER
20#define WIGNER3 fw3jja6
21#define WIGNER6 fw6jja
22#else
23#define WIGNER3 wig3jj
24#define WIGNER6 wig6jj
25#endif
26
27Numeric wigner3j(const Rational j1,
28 const Rational j2,
29 const Rational j3,
30 const Rational m1,
31 const Rational m2,
32 const Rational m3) {
33 errno = 0;
34
35 const int a = (2 * j1).toInt(), b = (2 * j2).toInt(), c = (2 * j3).toInt(),
36 d = (2 * m1).toInt(), e = (2 * m2).toInt(), f = (2 * m3).toInt();
37 double g;
38 const int j = std::max({std::abs(a),
39 std::abs(b),
40 std::abs(c),
41 std::abs(d),
42 std::abs(e),
43 std::abs(f)}) *
44 3 / 2 +
45 1;
46
47 wig_thread_temp_init(j);
48 g = WIGNER3(a, b, c, d, e, f);
49 wig_temp_free();
50
51 if (errno == EDOM) {
52 errno = 0;
53 ARTS_USER_ERROR("Bad state, perhaps you need to call Wigner3Init?")
54 }
55
56 return Numeric(g);
57}
58
59Numeric wigner6j(const Rational j1,
60 const Rational j2,
61 const Rational j3,
62 const Rational l1,
63 const Rational l2,
64 const Rational l3) {
65 errno = 0;
66
67 const int a = (2 * j1).toInt(), b = (2 * j2).toInt(), c = (2 * j3).toInt(),
68 d = (2 * l1).toInt(), e = (2 * l2).toInt(), f = (2 * l3).toInt();
69 double g;
70 const int j = std::max({std::abs(a),
71 std::abs(b),
72 std::abs(c),
73 std::abs(d),
74 std::abs(e),
75 std::abs(f)});
76
77 wig_thread_temp_init(j);
78 g = WIGNER6(a, b, c, d, e, f);
79 wig_temp_free();
80
81 if (errno == EDOM) {
82 errno = 0;
83 ARTS_USER_ERROR("Bad state, perhaps you need to call Wigner6Init?")
84 }
85
86 return Numeric(g);
87}
88
89std::pair<Rational, Rational> wigner_limits(std::pair<Rational, Rational> a,
90 std::pair<Rational, Rational> b) {
91 const bool invalid = a.first.isUndefined() or b.first.isUndefined();
92 if (invalid) {
93 return {RATIONAL_UNDEFINED, RATIONAL_UNDEFINED};
94 }
95
96 std::pair<Rational, Rational> out{max(a.first, b.first),
97 min(a.second, b.second)};
98 if (out.first > out.second) out = {RATIONAL_UNDEFINED, RATIONAL_UNDEFINED};
99 return out;
100}
101
102Index make_wigner_ready(int largest, [[maybe_unused]] int fastest, int size) {
103 if (size == 3) {
104#if DO_FAST_WIGNER
105 fastwigxj_load(FAST_WIGNER_PATH_3J, 3, NULL);
106#ifdef _OPENMP
107 fastwigxj_thread_dyn_init(3, fastest);
108#else
109 fastwigxj_dyn_init(3, fastest);
110#endif
111#endif
112 wig_table_init(largest, 3);
113
114 return largest;
115 }
116
117 if (size == 6) {
118#if DO_FAST_WIGNER
119 fastwigxj_load(FAST_WIGNER_PATH_3J, 3, NULL);
120 fastwigxj_load(FAST_WIGNER_PATH_6J, 6, NULL);
121#ifdef _OPENMP
122 fastwigxj_thread_dyn_init(3, fastest);
123 fastwigxj_thread_dyn_init(6, fastest);
124#else
125 fastwigxj_dyn_init(3, fastest);
126 fastwigxj_dyn_init(6, fastest);
127#endif
128#endif
129 wig_table_init(largest * 2, 6);
130
131 return largest;
132 }
133
134 return 0;
135}
136
137bool is_wigner_ready(int j) {
138 extern int wigxjpf_max_prime_decomp;
139 return not(j > wigxjpf_max_prime_decomp);
140}
141
142bool is_wigner3_ready(const Rational& J) {
143 const int test = J.toInt(6) / 2 + 1; // nb. J can be half-valued
144 return is_wigner_ready(test);
145}
146
147bool is_wigner6_ready(const Rational& J) {
148 const int test = J.toInt(4) + 1; // nb. J can be half-valued
149 return is_wigner_ready(test);
150}
151
152Numeric dwigner3j(Index M, Index J1, Index J2, Index J) {
153 auto CJM = [](Index j, Index m) {
154 Numeric cjm = 1.;
155 for (Index I = 1; I <= j; I++) cjm *= (1. - .5 / static_cast<Numeric>(I));
156 for (Index K = 1; K <= m; K++)
157 cjm *= static_cast<Numeric>(j + 1 - K) / static_cast<Numeric>(j + K);
158 return cjm;
159 };
160
161 Numeric GCM = 0.;
162 if (J1 < 0 or J2 < 0 or J < 0) return GCM;
163 Index JM = J1 + J2;
164 Index J0 = std::abs(J1 - J2);
165 if (J > JM or J < J0) return GCM;
166 Index JS = std::max(J1, J2);
167 Index JI = std::min(J1, J2);
168 Index MA = std::abs(M);
169 if (MA > JI) return GCM;
170 Index UN = 1 - 2 * (JS % 2);
171 Index QM = M + M;
172 Numeric CG0 = 0.;
173 GCM = static_cast<Numeric>(UN) *
174 std::sqrt(CJM(JI, MA) / CJM(JS, MA) * CJM(J0, 0) /
175 static_cast<Numeric>(JS + JS + 1));
176 Index AJ0 = J0;
177 Index AJM = JM + 1;
178 Index AJ02 = AJ0 * AJ0;
179 Index AJM2 = AJM * AJM;
180 Numeric ACG0 = 0.;
181 for (Index I = J0 + 1; I <= J; I++) {
182 Index AI = I;
183 Index AI2 = AI * AI;
184 Numeric ACG = std::sqrt((AJM2 - AI2) * (AI2 - AJ02));
185 Numeric CG1 =
186 (static_cast<Numeric>(QM) * static_cast<Numeric>(I + I - 1) * GCM -
187 ACG0 * CG0) /
188 ACG;
189 CG0 = GCM;
190 GCM = CG1;
191 ACG0 = ACG;
192 }
193 return GCM;
194}
195
196Numeric dwigner6j(Index A, Index B, Index C, Index D, Index F) {
197 Numeric SIXJ, TERM;
198 if (std::abs(A - C) > F or std::abs(B - D) > F or (A + C) < F or (B + D < F))
199 goto x1000;
200 switch (C - D + 2) {
201 case 2:
202 goto x2;
203 case 3:
204 goto x3;
205 case 4:
206 goto x1000;
207 default:
208 goto x1;
209 }
210
211x1:
212 switch (A - B + 2) {
213 case 2:
214 goto x11;
215 case 3:
216 goto x12;
217 default:
218 goto x10;
219 }
220
221x10:
222 TERM =
223 (static_cast<Numeric>(F + B + D + 1) * static_cast<Numeric>(F + B + D) *
224 static_cast<Numeric>(B + D - F) * static_cast<Numeric>(B + D - (1 + F)));
225 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B) *
226 static_cast<Numeric>(2 * B - 1) * static_cast<Numeric>(D) *
227 static_cast<Numeric>(2 * D - 1) * static_cast<Numeric>(2 * D + 1));
228 SIXJ = static_cast<Numeric>(pow_negative_one(A + C + F)) * std::sqrt(TERM);
229 return SIXJ;
230
231x11:
232 TERM = (static_cast<Numeric>(F + B + D + 1) *
233 static_cast<Numeric>(F + B - D + 1) *
234 static_cast<Numeric>(F + D - B) * static_cast<Numeric>(B + D - F));
235 TERM /= (4. * static_cast<Numeric>(B) * static_cast<Numeric>(2 * B + 1) *
236 static_cast<Numeric>(B + 1) * static_cast<Numeric>(D) *
237 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(2 * D - 1));
238 SIXJ =
239 static_cast<Numeric>(pow_negative_one(A + C - F - 1)) * std::sqrt(TERM);
240 return SIXJ;
241
242x12:
243 TERM =
244 (static_cast<Numeric>(F + D - B) * static_cast<Numeric>(F + D - B - 1) *
245 static_cast<Numeric>(F + B - D + 2) *
246 static_cast<Numeric>(F + B - D + 1));
247 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B + 1) *
248 static_cast<Numeric>(2 * B + 3) * static_cast<Numeric>(2 * D - 1) *
249 static_cast<Numeric>(D) * static_cast<Numeric>(2 * D + 1));
250 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
251 return SIXJ;
252
253x2:
254 switch (A - B + 2) {
255 case 2:
256 goto x21;
257 case 3:
258 goto x22;
259 default:
260 goto x20;
261 }
262
263x20:
264 TERM =
265 (static_cast<Numeric>(F + B + D + 1) * static_cast<Numeric>(D + B - F) *
266 static_cast<Numeric>(F + B - D) * static_cast<Numeric>(F + D - B + 1));
267 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B) *
268 static_cast<Numeric>(2 * B - 1) * static_cast<Numeric>(D) *
269 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(D + 1));
270 SIXJ =
271 static_cast<Numeric>(pow_negative_one(A + C - F - 1)) * std::sqrt(TERM);
272 return SIXJ;
273
274x21:
275 TERM = (static_cast<Numeric>(B) * static_cast<Numeric>(B + 1) +
276 static_cast<Numeric>(D) * static_cast<Numeric>(D + 1) -
277 static_cast<Numeric>(F) * static_cast<Numeric>(F + 1));
278 TERM /=
279 std::sqrt(4. * static_cast<Numeric>(B) * static_cast<Numeric>(B + 1) *
280 static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(D) *
281 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(D + 1));
282 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F - 1)) * TERM;
283 return SIXJ;
284
285x22:
286 TERM =
287 (static_cast<Numeric>(F + D + B + 2) *
288 static_cast<Numeric>(F + B - D + 1) *
289 static_cast<Numeric>(B + D - F + 1) * static_cast<Numeric>(F + D - B));
290 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B + 1) *
291 static_cast<Numeric>(2 * B + 3) * static_cast<Numeric>(D) *
292 static_cast<Numeric>(D + 1) * static_cast<Numeric>(2 * D + 1));
293 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
294 return SIXJ;
295
296x3:
297 switch (A - B + 2) {
298 case 2:
299 goto x31;
300 case 3:
301 goto x32;
302 default:
303 goto x30;
304 }
305
306x30:
307 TERM =
308 (static_cast<Numeric>(F + B - D) * static_cast<Numeric>(F + B - D - 1) *
309 static_cast<Numeric>(F + D - B + 2) *
310 static_cast<Numeric>(F + D - B + 1));
311 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B) *
312 static_cast<Numeric>(2 * B - 1) * static_cast<Numeric>(D + 1) *
313 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(2 * D + 3));
314 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
315 return SIXJ;
316
317x31:
318 TERM =
319 (static_cast<Numeric>(F + D + B + 2) *
320 static_cast<Numeric>(B + D - F + 1) *
321 static_cast<Numeric>(F + D - B + 1) * static_cast<Numeric>(F + B - D));
322 TERM /= (4. * static_cast<Numeric>(B) * static_cast<Numeric>(2 * B + 1) *
323 static_cast<Numeric>(B + 1) * static_cast<Numeric>(2 * D + 1) *
324 static_cast<Numeric>(D + 1) * static_cast<Numeric>(2 * D + 3));
325 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
326 return SIXJ;
327
328x32:
329 TERM = (static_cast<Numeric>(F + D + B + 3) *
330 static_cast<Numeric>(F + B + D + 2) *
331 static_cast<Numeric>(B + D - F + 2) *
332 static_cast<Numeric>(B + D - F + 1));
333 TERM /= (4. * static_cast<Numeric>(2 * B + 3) * static_cast<Numeric>(B + 1) *
334 static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(2 * D + 3) *
335 static_cast<Numeric>(D + 1) * static_cast<Numeric>(2 * D + 1));
336 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
337 return SIXJ;
338
339x1000:
340 SIXJ = 0.;
341 return SIXJ;
342}
base max(const Array< base > &x)
Max function.
Definition array.h:128
base min(const Array< base > &x)
Min function.
Definition array.h:144
Common ARTS conversions.
Header file for helper functions for OpenMP.
Helper macros for debugging.
#define ARTS_USER_ERROR(...)
Definition debug.h:153
constexpr Index pow_negative_one(Index x) noexcept
Computes std::pow(-1, x) without std::pow.
Definition math_funcs.h:143
#define d
#define a
#define c
#define b
bool is_wigner_ready(int j)
Tells if the function can deal with the input integer.
#define WIGNER6
Numeric dwigner6j(Index A, Index B, Index C, Index D, Index F)
Computes the wigner 6J symbol with floating point precision.
Numeric wigner6j(const Rational j1, const Rational j2, const Rational j3, const Rational l1, const Rational l2, const Rational l3)
Wigner 6J symbol.
bool is_wigner3_ready(const Rational &J)
Tells if the function is ready for Wigner 3J calculations.
Index make_wigner_ready(int largest, int fastest, int size)
Ready Wigner.
Numeric dwigner3j(Index M, Index J1, Index J2, Index J)
Computes the wigner 3J symbol with floating point precision.
bool is_wigner6_ready(const Rational &J)
Tells if the function is ready for Wigner 6J calculations.
Numeric wigner3j(const Rational j1, const Rational j2, const Rational j3, const Rational m1, const Rational m2, const Rational m3)
Wigner 3J symbol.
#define WIGNER3
std::pair< Rational, Rational > wigner_limits(std::pair< Rational, Rational > a, std::pair< Rational, Rational > b)
Wigner symbol interactions.