| Geant4 Cross Reference |
1 /* 1
2 * calculate coupling coefficients of angular
3 *
4 * Author:
5 * Kawano, T <kawano@mailaps.or
6 *
7 * Modified by David Brown <dbrown@bnl.gov>
8 * No longer must precompute the logarith
9 * Also renamed things to make more Pytho
10 * Finally, fixed a bunch of bugs & confu
11 *
12 * Functions:
13 *
14 * Note that arguments of those functions mus
15 *
16 * wigner_3j(j1,j2,j3,j4,j5,j6)
17 * Wigner's 3J symbol (similar to Clebsh-
18 * = / j1 j2 j3 \
19 * \ j4 j5 j6 /
20 *
21 * wigner_6j(j1,j2,j3,j4,j5,j6)
22 * Wigner's 6J symbol (similar to Racah)
23 * = { j1 j2 j3 }
24 * { j4 j5 j6 }
25 *
26 * wigner_9j(j1,j2,j3,j4,j5,j6,j7,j8,j9)
27 * Wigner's 9J symbol
28 * / j1 j2 j3 \
29 * = | j4 j5 j6 |
30 * \ j7 j8 j9 /
31 *
32 * racah(j1, j2, l2, l1, j3, l3)
33 * = W(j1, j2, l2, l1 ; j3, l3)
34 * = (-1)^(j1+j2+l1+l2) * { j1 j2
35 * { l1 l2
36 *
37 * clebsh_gordan(j1,j2,m1,m2,j3)
38 * Clebsh-Gordan coefficient
39 * = <j1,j2,m1,m2|j3,m1+m2>
40 * = (-)^(j1-j2+m1+m2) * std::sqr
41 *
42 *
43 * z_coefficient(l1,j1,l2,j2,S,L)
44 * Biedenharn's Z-coefficient coefficient
45 * = Z(l1 j1 l2 j2 | S L )
46 *
47 * reduced_matrix_element(L,S,J,l0,j0,l1,j1)
48 * Reduced Matrix Element for Tensor Oper
49 * = < l1j1 || T(YL,sigma_S)J ||
50 *
51 * References:
52 * A. R. Edmonds, Angular Momentum in Quantum
53 * E. Condon, and G. Shortley, The Theory of
54 */
55
56 #include <stdlib.h>
57 #define _USE_MATH_DEFINES
58 #include <cmath>
59
60 #include "nf_specialFunctions.h"
61
62 #if defined __cplusplus
63 #include <cmath>
64 #include "G4Exp.hh"
65 namespace GIDI {
66 using namespace GIDI;
67 #endif
68
69 static const int MAX_FACTORIAL = 200;
70 /*static const double ARRAY_OVER = 1.0e+300;
71 static const double nf_amc_log_fact[] = {0.0,
72
73 static int parity( int x );
74 static int max3( int a, int b, int c );
75 static int max4( int a, int b, int c, int d );
76 static int min3( int a, int b, int c );
77 static double w6j0( int, int * );
78 static double w6j1( int * );
79 static double cg1( int, int, int );
80 static double cg2( int, int, int, int, int, in
81 static double cg3( int, int, int, int, int, in
82 /*static double triangle( int, int, int );*/
83 /*
84 ==============================================
85 */
86 double nf_amc_log_factorial( int n ) {
87 /*
88 * returns ln( n! ).
89 */
90 if( n > MAX_FACTORIAL ) return( INFINITY )
91 if( n < 0 ) return( INFINITY );
92 return nf_amc_log_fact[n];
93 }
94 /*
95 ==============================================
96 */
97 double nf_amc_factorial( int n ) {
98 /*
99 * returns n! for pre-computed table. INF
100 */
101 return G4Exp( nf_amc_log_factorial( n ) );
102 }
103 /*
104 ==============================================
105 */
106 double nf_amc_wigner_3j( int j1, int j2, int j
107 /*
108 * Wigner's 3J symbol (similar to Clebsh-
109 * = / j1 j2 j3 \
110 * \ j4 j5 j6 /
111 */
112 double cg;
113
114 if( ( j4 + j5 + j6 ) != 0 ) return( 0.0 );
115 if( ( cg = nf_amc_clebsh_gordan( j1, j2, j
116 if( cg == INFINITY ) return( cg );
117 return( ( ( ( j1 - j2 - j6 ) % 4 == 0 ) ?
118 }
119 /*
120 ==============================================
121 */
122 double nf_amc_wigner_6j( int j1, int j2, int j
123 /*
124 * Wigner's 6J symbol (similar to Racah)
125 * = { j1 j2 j3 }
126 * { j4 j5 j6 }
127 */
128 int i, x[6];
129
130 x[0] = j1; x[1] = j2; x[2] = j3; x[3] = j4
131 for( i = 0; i < 6; i++ ) if ( x[i] == 0 )
132
133 return( w6j1( x ) );
134 }
135 /*
136 ==============================================
137 */
138 static double w6j0( int i, int *x ) {
139
140 switch( i ){
141 case 0: if ( ( x[1] != x[2] ) || ( x[4]
142 x[5] = x[3]; x[0] = x[1]; x[3]
143 case 1: if ( ( x[0] != x[2] ) || ( x[3]
144 x[5] = x[4];
145 case 2: if ( ( x[0] != x[1] ) || ( x[3]
146 break;
147 //TK fix bug and add comment on 17-05
148 //This is the case of 6.3.2 of
149 case 3: if ( ( x[1] != x[5] ) || ( x[2]
150 x[5] = x[0]; x[0] = x[4]; x[3]
151 case 4: if ( ( x[0] != x[5] ) || ( x[2]
152 x[5] = x[1];
153 case 5: if ( ( x[0] != x[4] ) || ( x[1]
154 x[5] = x[2];
155 }
156
157 if( ( x[5] > ( x[0] + x[3] ) ) || ( x[5] <
158 if( x[0] > MAX_FACTORIAL || x[3] > MAX_FAC
159 return( INFINITY );
160 }
161
162 return( 1.0 / std::sqrt( (double) ( ( x[0]
163 }
164 /*
165 ==============================================
166 */
167 static double w6j1( int *x ) {
168
169 double w6j, w;
170 int i, k, k1, k2, n, l1, l2, l3, l4, n1, n
171 static int a[3][4] = { { 0, 0, 3, 3},
172 { 1, 4, 1, 4},
173 { 2, 5, 5, 2} };
174
175 w6j = 0.0;
176
177 for ( k = 0; k < 4; k++ ){
178 x1 = x[ ( a[0][k] ) ];
179 x2 = x[ ( a[1][k] ) ];
180 x3 = x[ ( a[2][k] ) ];
181
182 n = ( x1 + x2 + x3 ) / 2;
183 if( n > MAX_FACTORIAL ) {
184 return( INFINITY ); }
185 else if( n < 0 ) {
186 return( 0.0 );
187 }
188
189 if ( ( n1 = n - x3 ) < 0 ) return( 0.0
190 if ( ( n2 = n - x2 ) < 0 ) return( 0.0
191 if ( ( n3 = n - x1 ) < 0 ) return( 0.0
192
193 y[k] = n + 2;
194 w6j += nf_amc_log_fact[n1] + nf_amc_lo
195 }
196
197 n1 = ( x[0] + x[1] + x[3] + x[4] ) / 2;
198 n2 = ( x[0] + x[2] + x[3] + x[5] ) / 2;
199 n3 = ( x[1] + x[2] + x[4] + x[5] ) / 2;
200
201 k1 = max4( y[0], y[1], y[2], y[3] ) - 1;
202 k2 = min3( n1, n2, n3 ) + 1;
203
204 l1 = k1 - y[0] + 1; m1 = n1 - k1 + 1;
205 l2 = k1 - y[1] + 1; m2 = n2 - k1 + 1;
206 l3 = k1 - y[2] + 1; m3 = n3 - k1 + 1;
207 l4 = k1 - y[3] + 1;
208
209 w6j = w = G4Exp( 0.5 * w6j + nf_amc_log_fa
210 - nf_amc_log_fact
211 if( w6j == INFINITY ) return( INFINITY );
212
213 if( k1 != k2 ){
214 k = k2 - k1;
215 m1 -= k-1; m2 -= k-1; m3 -= k-1;
216 l1 += k ; l2 += k ; l3 += k ; l4 +=
217
218 for ( i = 0; i < k; i++ )
219 w6j = w - w6j * ( ( k2 - i ) * ( m
220 / ( ( l1 - i ) * ( l2 - i
221 }
222 return( w6j );
223 }
224 /*
225 ==============================================
226 */
227 double nf_amc_wigner_9j( int j1, int j2, int j
228 /*
229 * Wigner's 9J symbol
230 * / j1 j2 j3 \
231 * = | j4 j5 j6 |
232 * \ j7 j8 j9 /
233 *
234 */
235 int i, i0, i1;
236 double rac;
237
238 i0 = max3( std::abs( j1 - j9 ), std::abs(
239 i1 = min3( ( j1 + j9 ), ( j2 + j6
240
241 rac = 0.0;
242 for ( i = i0; i <= i1; i += 2 ){
243 rac += nf_amc_racah( j1, j4, j9, j8, j
244 * nf_amc_racah( j2, j5, i, j4, j
245 * nf_amc_racah( j9, i, j3, j2, j
246 if( rac == INFINITY ) return( INFINITY
247 }
248
249 return( ( ( (int)( ( j1 + j3 + j5 + j8 ) /
250 }
251 /*
252 ==============================================
253 */
254 double nf_amc_racah( int j1, int j2, int l2, i
255 /*
256 * Racah coefficient definition in Edmond
257 * W(j1, j2, l2, l1 ; j3, l3) = (-1)^(j1+
258 *
259 * The call signature of W(...) appears j
260 *
261 * This convention is exactly that used b
262 */
263
264 double sig;
265
266 sig = ( ( ( j1 + j2 + l1 + l2 ) % 4 == 0 )
267 return sig * nf_amc_wigner_6j( j1, j2, j3,
268 }
269
270 /*
271 ==============================================
272 */
273 /*
274 static double triangle( int a, int b, int c )
275
276 int j1, j2, j3, j4;
277
278 if ( ( j1 = ( a + b - c ) / 2 ) < 0 ) retu
279 if ( ( j2 = ( a - b + c ) / 2 ) < 0 ) retu
280 if ( ( j3 = ( -a + b + c ) / 2 ) < 0 ) retu
281 j4 = ( a + b + c ) / 2 + 1;
282
283 return( std::exp( 0.5 * ( nf_amc_log_fact[j
284 }
285 */
286 /*
287 ==============================================
288 */
289 double nf_amc_clebsh_gordan( int j1, int j2, i
290 /*
291 * Clebsh-Gordan coefficient
292 * = <j1,j2,m1,m2|j3,m1+m2>
293 * = (-)^(j1-j2+m1+m2) * std::sqrt(2*
294 *
295 *
296 * Note: Last value m3 is preset to m1+m2. A
297 */
298
299 int m3, x1, x2, x3, y1, y2, y3;
300 double cg = 0.0;
301
302 if ( j1 < 0 || j2 < 0 || j3 < 0) return( 0
303 if ( j1 + j2 + j3 > 2 * MAX_FACTORIAL ) re
304
305 m3 = m1 + m2;
306
307 if ( ( x1 = ( j1 + m1 ) / 2 + 1 ) <= 0 ) r
308 if ( ( x2 = ( j2 + m2 ) / 2 + 1 ) <= 0 ) r
309 if ( ( x3 = ( j3 - m3 ) / 2 + 1 ) <= 0 ) r
310
311 if ( ( y1 = x1 - m1 ) <= 0 ) return( 0.0 )
312 if ( ( y2 = x2 - m2 ) <= 0 ) return( 0.0 )
313 if ( ( y3 = x3 + m3 ) <= 0 ) return( 0.0 )
314
315 if ( j3 == 0 ){
316 if ( j1 == j2 ) cg = ( 1.0 / std::sqrt
317 }
318 else if ( (j1 == 0 || j2 == 0 ) ){
319 if ( ( j1 + j2 ) == j3 ) cg = 1.0;
320 }
321 else {
322 if( m3 == 0 && std::abs( m1 ) <= 1 ){
323 if( m1 == 0 ) cg = cg1( x1, x2, x3
324 else cg = cg2( x1 + y1 -
325 }
326 else if ( m2 == 0 && std::abs( m1 ) <=
327 cg = cg2( x1 - y2 +
328 }
329 else if ( m1 == 0 && std::abs( m3 ) <=
330 cg = cg2( x1, x1 - 1
331 }
332 else cg = cg3( x1, x2, x3
333 }
334
335 return( cg );
336 }
337 /*
338 ==============================================
339 */
340 static double cg1( int x1, int x2, int x3 ) {
341
342 int p1, p2, p3, p4, q1, q2, q3, q4;
343 double a;
344
345 p1 = x1 + x2 + x3 - 1; if ( ( p1 % 2 ) !=
346 p2 = x1 + x2 - x3;
347 p3 =-x1 + x2 + x3;
348 p4 = x1 - x2 + x3;
349 if ( p2 <= 0 || p3 <= 0 || p4 <= 0 ) retur
350 if ( p1 >= MAX_FACTORIAL ) return( INFINIT
351
352 q1 = ( p1 + 1 ) / 2 - 1; p1--;
353 q2 = ( p2 + 1 ) / 2 - 1; p2--;
354 q3 = ( p3 + 1 ) / 2 - 1; p3--;
355 q4 = ( p4 + 1 ) / 2 - 1; p4--;
356
357 a = nf_amc_log_fact[q1]-( nf_amc_log_fact[
358 + 0.5 * ( nf_amc_log_fact[ 2 * x3 - 1
359 + nf_amc_log_fact[p2] + nf_amc_log_fac
360
361 return( ( ( ( q1 + x1 - x2 ) % 2 == 0 ) ?
362 }
363 /*
364 ==============================================
365 */
366 static double cg2( int k, int q0, int z1, int
367
368 int q1, q2, q3, q4, p1, p2, p3, p4;
369 double a;
370
371 p1 = z1 + q0 + 2;
372 p2 = z1 - q0 + 1;
373 p3 = z2 + q0 + 1;
374 p4 = -z2 + q0 + 1;
375 if ( p2 <= 0 || p3 <= 0 || p4 <= 0) return
376 if ( p1 >= MAX_FACTORIAL ) return( INFINIT
377
378 q1 = ( p1 + 1 ) / 2 - 1; p1--;
379 q2 = ( p2 + 1 ) / 2 - 1; p2--;
380 q3 = ( p3 + 1 ) / 2 - 1; p3--;
381 q4 = ( p4 + 1 ) / 2 - 1; p4--;
382
383 a = nf_amc_log_fact[q1] - ( nf_amc_log_fac
384 + 0.5 * ( nf_amc_log_fact[ w3 + 1 ] -
385 + nf_amc_log_fact[ w1 ] -
386 + nf_amc_log_fact[ w2 ] -
387 + nf_amc_log_fact[ p2 ] +
388
389 return( ( ( ( q4 + k + ( mm > 0 ) * ( p1 +
390 }
391 /*
392 ==============================================
393 */
394 static double cg3( int x1, int x2, int x3, int
395
396 int nx, i, k1, k2, q1, q2, q3, q4, p1, p2,
397 double a, cg;
398
399 nx = x1 + x2 + x3 - 1;
400 if ( ( z1 = nx - x1 - y1 ) < 0 ) return( 0
401 if ( ( z2 = nx - x2 - y2 ) < 0 ) return( 0
402 if ( ( z3 = nx - x3 - y3 ) < 0 ) return( 0
403
404 k1 = x2 - y3;
405 k2 = y1 - x3;
406
407 q1 = max3( k1, k2, 0 );
408 q2 = min3( y1, x2, z3 + 1 ) - 1;
409 q3 = q1 - k1;
410 q4 = q1 - k2;
411
412 p1 = y1 - q1 - 1;
413 p2 = x2 - q1 - 1;
414 p3 = z3 - q1;
415
416 a = cg = G4Exp( 0.5 * ( nf_amc_log_fact[ x
417 + nf_amc_log_fact[ z1 ] +
418 + nf_amc_log_fact[ x1 - 1 ] +
419 + nf_amc_log_fact[ y1 - 1 ] +
420 - nf_amc_log_fact[ p1 ] -
421 - nf_amc_log_fact[ q1 ] -
422 if( cg == INFINITY ) return( INFINITY );
423
424 if ( q1 != q2 ){
425 q3 = q2 - k1;
426 q4 = q2 - k2;
427 p1 = y1 - q2;
428 p2 = x2 - q2;
429 p3 = z3 - q2 + 1;
430 for( i = 0; i < ( q2 - q1 ); i++ )
431 cg = a - cg * ( ( p1 + i ) * ( p2
432 }
433 return( cg );
434 }
435 /*
436 ==============================================
437 */
438 double nf_amc_z_coefficient( int l1, int j1, i
439 /*
440 * Biedenharn's Z-coefficient coefficient
441 * = Z(l1 j1 l2 j2 | S L )
442 */
443 double z, clebsh_gordan = nf_amc_clebsh_go
444
445 if( ( clebsh_gordan == INFINITY ) || ( rac
446 z = ( ( ( -l1 + l2 + ll ) % 8 == 0 ) ? 1.0
447 * std::sqrt( l1 + 1.0 ) * std::sqrt( l
448
449 return( z );
450 }
451 /*
452 ==============================================
453 */
454 double nf_amc_zbar_coefficient( int l1, int j1
455 /*
456 * Lane & Thomas's Zbar-coefficient coeff
457 * = Zbar(l1 j1 l2 j2 | S L )
458 * = (-i)^( -l1 + l2 + ll ) * Z(l1 j
459 *
460 * Lane & Thomas Rev. Mod. Phys. 30, 257-
461 * Note, Lane & Thomas define this becaus
462 * Froehner uses Lane & Thomas convention
463 */
464 double zbar, clebsh_gordan = nf_amc_clebsh
465
466 if( ( clebsh_gordan == INFINITY ) || ( rac
467 zbar = std::sqrt( l1 + 1.0 ) * std::sqrt(
468
469 return( zbar );
470 }
471 /*
472 ==============================================
473 */
474 double nf_amc_reduced_matrix_element( int lt,
475 /*
476 * Reduced Matrix Element for Tensor Oper
477 * = < l1j1 || T(YL,sigma_S)J || l0j0
478 *
479 * M.B.Johnson, L.W.Owen, G.R.Satchler
480 * Phys. Rev. 142, 748 (1966)
481 * Note: definition differs from JOS by the f
482 */
483 int llt;
484 double x1, x2, x3, reduced_mat, clebsh_gor
485
486 if ( parity( lt ) != parity( l0 ) * parity
487 if ( std::abs( l0 - l1 ) > lt || ( l0 + l1
488 if ( std::abs( ( j0 - j1 ) / 2 ) > jt || (
489
490 llt = 2 * lt;
491 jt *= 2;
492 st *= 2;
493
494 if( ( clebsh_gordan = nf_amc_clebsh_gordan
495
496 reduced_mat = 1.0 / std::sqrt( 4 * M_PI )
497 * std::sqrt( ( j0 + 1.0 ) * (
498 * parity( ( j1 - j0 ) / 2 ) *
499
500 if( st == 2 ){
501 x1 = ( l0 - j0 / 2.0 ) * ( j0 + 1.0 );
502 x2 = ( l1 - j1 / 2.0 ) * ( j1 + 1.0 );
503 if ( jt == llt ){
504 x3 = ( lt == 0 ) ? 0 : ( x1 - x2
505 }
506 else if ( jt == ( llt - st ) ){
507 x3 = ( lt == 0 ) ? 0 : -( lt + x1
508 }
509 else if ( jt == ( llt + st ) ){
510 x3 = ( lt + 1 - x1 - x2 ) / std::s
511 }
512 else{
513 x3 = 1.0;
514 }
515 }
516 else x3 = 1.0;
517 reduced_mat *= x3;
518
519 return( reduced_mat );
520 }
521 /*
522 ==============================================
523 */
524 static int parity( int x ) {
525
526 return( ( ( x / 2 ) % 2 == 0 ) ? 1 : -1 );
527 }
528 /*
529 ==============================================
530 */
531 static int max3( int a, int b, int c ) {
532
533 if( a < b ) a = b;
534 if( a < c ) a = c;
535 return( a );
536 }
537 /*
538 ==============================================
539 */
540 static int max4( int a, int b, int c, int d )
541
542 if( a < b ) a = b;
543 if( a < c ) a = c;
544 if( a < d ) a = d;
545 return( a );
546 }
547 /*
548 ==============================================
549 */
550 static int min3( int a, int b, int c ) {
551
552 if( a > b ) a = b;
553 if( a > c ) a = c;
554 return( a );
555 }
556
557 #if defined __cplusplus
558 }
559 #endif
560