ARTS 2.5.4 (git: 31ce4f0e)
wigner_functions.cc
Go to the documentation of this file.
1/* Copyright (C) 2012
2 * Richard Larsson <ric.larsson@gmail.com>
3 *
4 * This program is free software; you can redistribute it and/or modify it
5 * under the terms of the GNU General Public License as published by the
6 * Free Software Foundation; either version 2, or (at your option) any
7 * later version.
8 *
9 * This program is distributed in the hope that it will be useful,
10 * but WITHOUT ANY WARRANTY; without even the implied warranty of
11 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 * GNU General Public License for more details.
13 *
14 * You should have received a copy of the GNU General Public License
15 * along with this program; if not, write to the Free Software
16 * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307,
17 * USA. */
18
27#include "wigner_functions.h"
28
29#include <sys/errno.h>
30
31#include <algorithm>
32
33#include "arts_omp.h"
34#include "arts_conversions.h"
35#include "debug.h"
36
37#if DO_FAST_WIGNER
38#define WIGNER3 fw3jja6
39#define WIGNER6 fw6jja
40#else
41#define WIGNER3 wig3jj
42#define WIGNER6 wig6jj
43#endif
44
46 const Rational j2,
47 const Rational j3,
48 const Rational m1,
49 const Rational m2,
50 const Rational m3) {
51 errno = 0;
52
53 const int a = (2 * j1).toInt(), b = (2 * j2).toInt(), c = (2 * j3).toInt(),
54 d = (2 * m1).toInt(), e = (2 * m2).toInt(), f = (2 * m3).toInt();
55 double g;
56 const int j = std::max({std::abs(a),
57 std::abs(b),
58 std::abs(c),
59 std::abs(d),
60 std::abs(e),
61 std::abs(f)}) *
62 3 / 2 +
63 1;
64
65 wig_thread_temp_init(j);
66 g = WIGNER3(a, b, c, d, e, f);
67 wig_temp_free();
68
69 if (errno == EDOM) {
70 errno = 0;
71 ARTS_USER_ERROR("Bad state, perhaps you need to call Wigner3Init?")
72 }
73
74 return Numeric(g);
75}
76
78 const Rational j2,
79 const Rational j3,
80 const Rational l1,
81 const Rational l2,
82 const Rational l3) {
83 errno = 0;
84
85 const int a = (2 * j1).toInt(), b = (2 * j2).toInt(), c = (2 * j3).toInt(),
86 d = (2 * l1).toInt(), e = (2 * l2).toInt(), f = (2 * l3).toInt();
87 double g;
88 const int j = std::max({std::abs(a),
89 std::abs(b),
90 std::abs(c),
91 std::abs(d),
92 std::abs(e),
93 std::abs(f)});
94
95 wig_thread_temp_init(j);
96 g = WIGNER6(a, b, c, d, e, f);
97 wig_temp_free();
98
99 if (errno == EDOM) {
100 errno = 0;
101 ARTS_USER_ERROR("Bad state, perhaps you need to call Wigner6Init?")
102 }
103
104 return Numeric(g);
105}
106
107std::pair<Rational, Rational> wigner_limits(std::pair<Rational, Rational> a,
108 std::pair<Rational, Rational> b) {
109 const bool invalid = a.first.isUndefined() or b.first.isUndefined();
110 if (invalid) {
112 }
113
114 std::pair<Rational, Rational> out{max(a.first, b.first),
115 min(a.second, b.second)};
116 if (out.first > out.second) out = {RATIONAL_UNDEFINED, RATIONAL_UNDEFINED};
117 return out;
118}
119
120Index make_wigner_ready(int largest, [[maybe_unused]] int fastest, int size) {
121 if (size == 3) {
122#if DO_FAST_WIGNER
123 fastwigxj_load(FAST_WIGNER_PATH_3J, 3, NULL);
124#ifdef _OPENMP
125 fastwigxj_thread_dyn_init(3, fastest);
126#else
127 fastwigxj_dyn_init(3, fastest);
128#endif
129#endif
130 wig_table_init(largest, 3);
131
132 return largest;
133 }
134
135 if (size == 6) {
136#if DO_FAST_WIGNER
137 fastwigxj_load(FAST_WIGNER_PATH_3J, 3, NULL);
138 fastwigxj_load(FAST_WIGNER_PATH_6J, 6, NULL);
139#ifdef _OPENMP
140 fastwigxj_thread_dyn_init(3, fastest);
141 fastwigxj_thread_dyn_init(6, fastest);
142#else
143 fastwigxj_dyn_init(3, fastest);
144 fastwigxj_dyn_init(6, fastest);
145#endif
146#endif
147 wig_table_init(largest * 2, 6);
148
149 return largest;
150 }
151
152 return 0;
153}
154
155bool is_wigner_ready(int j) {
156 extern int wigxjpf_max_prime_decomp;
157 return not(j > wigxjpf_max_prime_decomp);
158}
159
161 const int test = J.toInt(6) / 2 + 1; // nb. J can be half-valued
162 return is_wigner_ready(test);
163}
164
166 const int test = J.toInt(4) + 1; // nb. J can be half-valued
167 return is_wigner_ready(test);
168}
169
171 auto CJM = [](Index J, Index M) {
172 Numeric CJM = 1.;
173 for (Index I = 1; I <= J; I++) CJM *= (1. - .5 / static_cast<Numeric>(I));
174 for (Index K = 1; K <= M; K++)
175 CJM *= static_cast<Numeric>(J + 1 - K) / static_cast<Numeric>(J + K);
176 return CJM;
177 };
178
179 Numeric GCM = 0.;
180 if (J1 < 0 or J2 < 0 or J < 0) return GCM;
181 Index JM = J1 + J2;
182 Index J0 = std::abs(J1 - J2);
183 if (J > JM or J < J0) return GCM;
184 Index JS = std::max(J1, J2);
185 Index JI = std::min(J1, J2);
186 Index MA = std::abs(M);
187 if (MA > JI) return GCM;
188 Index UN = 1 - 2 * (JS % 2);
189 Index QM = M + M;
190 Numeric CG0 = 0.;
191 GCM = static_cast<Numeric>(UN) *
192 std::sqrt(CJM(JI, MA) / CJM(JS, MA) * CJM(J0, 0) /
193 static_cast<Numeric>(JS + JS + 1));
194 Index AJ0 = J0;
195 Index AJM = JM + 1;
196 Index AJ02 = AJ0 * AJ0;
197 Index AJM2 = AJM * AJM;
198 Numeric ACG0 = 0.;
199 for (Index I = J0 + 1; I <= J; I++) {
200 Index AI = I;
201 Index AI2 = AI * AI;
202 Numeric ACG = std::sqrt((AJM2 - AI2) * (AI2 - AJ02));
203 Numeric CG1 =
204 (static_cast<Numeric>(QM) * static_cast<Numeric>(I + I - 1) * GCM -
205 ACG0 * CG0) /
206 ACG;
207 CG0 = GCM;
208 GCM = CG1;
209 ACG0 = ACG;
210 }
211 return GCM;
212}
213
215 Numeric SIXJ, TERM;
216 if (std::abs(A - C) > F or std::abs(B - D) > F or (A + C) < F or (B + D < F))
217 goto x1000;
218 switch (C - D + 2) {
219 case 2:
220 goto x2;
221 case 3:
222 goto x3;
223 case 4:
224 goto x1000;
225 default:
226 goto x1;
227 }
228
229x1:
230 switch (A - B + 2) {
231 case 2:
232 goto x11;
233 case 3:
234 goto x12;
235 default:
236 goto x10;
237 }
238
239x10:
240 TERM =
241 (static_cast<Numeric>(F + B + D + 1) * static_cast<Numeric>(F + B + D) *
242 static_cast<Numeric>(B + D - F) * static_cast<Numeric>(B + D - (1 + F)));
243 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B) *
244 static_cast<Numeric>(2 * B - 1) * static_cast<Numeric>(D) *
245 static_cast<Numeric>(2 * D - 1) * static_cast<Numeric>(2 * D + 1));
246 SIXJ = static_cast<Numeric>(pow_negative_one(A + C + F)) * std::sqrt(TERM);
247 return SIXJ;
248
249x11:
250 TERM = (static_cast<Numeric>(F + B + D + 1) *
251 static_cast<Numeric>(F + B - D + 1) *
252 static_cast<Numeric>(F + D - B) * static_cast<Numeric>(B + D - F));
253 TERM /= (4. * static_cast<Numeric>(B) * static_cast<Numeric>(2 * B + 1) *
254 static_cast<Numeric>(B + 1) * static_cast<Numeric>(D) *
255 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(2 * D - 1));
256 SIXJ =
257 static_cast<Numeric>(pow_negative_one(A + C - F - 1)) * std::sqrt(TERM);
258 return SIXJ;
259
260x12:
261 TERM =
262 (static_cast<Numeric>(F + D - B) * static_cast<Numeric>(F + D - B - 1) *
263 static_cast<Numeric>(F + B - D + 2) *
264 static_cast<Numeric>(F + B - D + 1));
265 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B + 1) *
266 static_cast<Numeric>(2 * B + 3) * static_cast<Numeric>(2 * D - 1) *
267 static_cast<Numeric>(D) * static_cast<Numeric>(2 * D + 1));
268 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
269 return SIXJ;
270
271x2:
272 switch (A - B + 2) {
273 case 2:
274 goto x21;
275 case 3:
276 goto x22;
277 default:
278 goto x20;
279 }
280
281x20:
282 TERM =
283 (static_cast<Numeric>(F + B + D + 1) * static_cast<Numeric>(D + B - F) *
284 static_cast<Numeric>(F + B - D) * static_cast<Numeric>(F + D - B + 1));
285 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B) *
286 static_cast<Numeric>(2 * B - 1) * static_cast<Numeric>(D) *
287 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(D + 1));
288 SIXJ =
289 static_cast<Numeric>(pow_negative_one(A + C - F - 1)) * std::sqrt(TERM);
290 return SIXJ;
291
292x21:
293 TERM = (static_cast<Numeric>(B) * static_cast<Numeric>(B + 1) +
294 static_cast<Numeric>(D) * static_cast<Numeric>(D + 1) -
295 static_cast<Numeric>(F) * static_cast<Numeric>(F + 1));
296 TERM /=
297 std::sqrt(4. * static_cast<Numeric>(B) * static_cast<Numeric>(B + 1) *
298 static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(D) *
299 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(D + 1));
300 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F - 1)) * TERM;
301 return SIXJ;
302
303x22:
304 TERM =
305 (static_cast<Numeric>(F + D + B + 2) *
306 static_cast<Numeric>(F + B - D + 1) *
307 static_cast<Numeric>(B + D - F + 1) * static_cast<Numeric>(F + D - B));
308 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B + 1) *
309 static_cast<Numeric>(2 * B + 3) * static_cast<Numeric>(D) *
310 static_cast<Numeric>(D + 1) * static_cast<Numeric>(2 * D + 1));
311 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
312 return SIXJ;
313
314x3:
315 switch (A - B + 2) {
316 case 2:
317 goto x31;
318 case 3:
319 goto x32;
320 default:
321 goto x30;
322 }
323
324x30:
325 TERM =
326 (static_cast<Numeric>(F + B - D) * static_cast<Numeric>(F + B - D - 1) *
327 static_cast<Numeric>(F + D - B + 2) *
328 static_cast<Numeric>(F + D - B + 1));
329 TERM /= (4. * static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(B) *
330 static_cast<Numeric>(2 * B - 1) * static_cast<Numeric>(D + 1) *
331 static_cast<Numeric>(2 * D + 1) * static_cast<Numeric>(2 * D + 3));
332 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
333 return SIXJ;
334
335x31:
336 TERM =
337 (static_cast<Numeric>(F + D + B + 2) *
338 static_cast<Numeric>(B + D - F + 1) *
339 static_cast<Numeric>(F + D - B + 1) * static_cast<Numeric>(F + B - D));
340 TERM /= (4. * static_cast<Numeric>(B) * static_cast<Numeric>(2 * B + 1) *
341 static_cast<Numeric>(B + 1) * static_cast<Numeric>(2 * D + 1) *
342 static_cast<Numeric>(D + 1) * static_cast<Numeric>(2 * D + 3));
343 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
344 return SIXJ;
345
346x32:
347 TERM = (static_cast<Numeric>(F + D + B + 3) *
348 static_cast<Numeric>(F + B + D + 2) *
349 static_cast<Numeric>(B + D - F + 2) *
350 static_cast<Numeric>(B + D - F + 1));
351 TERM /= (4. * static_cast<Numeric>(2 * B + 3) * static_cast<Numeric>(B + 1) *
352 static_cast<Numeric>(2 * B + 1) * static_cast<Numeric>(2 * D + 3) *
353 static_cast<Numeric>(D + 1) * static_cast<Numeric>(2 * D + 1));
354 SIXJ = static_cast<Numeric>(pow_negative_one(A + C - F)) * std::sqrt(TERM);
355 return SIXJ;
356
357x1000:
358 SIXJ = 0.;
359 return SIXJ;
360}
Common ARTS conversions.
Header file for helper functions for OpenMP.
Helper macros for debugging.
#define ARTS_USER_ERROR(...)
Definition: debug.h:150
#define abs(x)
#define min(a, b)
#define max(a, b)
constexpr Index pow_negative_one(Index x) noexcept
Computes std::pow(-1, x) without std::pow.
Definition: math_funcs.h:158
NUMERIC Numeric
The type to use for all floating point numbers.
Definition: matpack.h:33
INDEX Index
The type to use for all integer numbers and indices.
Definition: matpack.h:39
constexpr Numeric e
Elementary charge convenience name [C].
#define RATIONAL_UNDEFINED
Definition: rational.h:357
Numeric sqrt(const Rational r)
Square root.
Definition: rational.h:705
#define M
Definition: rng.cc:165
Implements rational numbers to work with other ARTS types.
Definition: rational.h:43
constexpr int toInt(int n=1) const noexcept
Converts the value to int by n-scaled division in Index form.
Definition: rational.h:122
#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.