Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/lend/src/nf_angularMomentumCoupling.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 /*
  2 *   calculate coupling coefficients of angular momenta
  3 *
  4 *   Author:
  5 *                 Kawano, T <kawano@mailaps.org>
  6 *
  7 *   Modified by David Brown <dbrown@bnl.gov>
  8 *       No longer must precompute the logarithm of the factorials.  
  9 *       Also renamed things to make more Python friendly.
 10 *       Finally, fixed a bunch of bugs & confusing conventions
 11 *
 12 *   Functions:
 13 *
 14 *   Note that arguments of those functions must be doubled, namely 1/2 is 1, etc.
 15 *
 16 *   wigner_3j(j1,j2,j3,j4,j5,j6)
 17 *       Wigner's 3J symbol (similar to Clebsh-Gordan)
 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 j3 }
 35 *                                      { l1 l2 l3 }
 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::sqrt(2*j3+1) * / j1 j2   j3   \
 41 *                                                    \ m1 m2 -m1-m2 /
 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 Operator
 49 *               = < l1j1 || T(YL,sigma_S)J || l0j0 >
 50 *
 51 * References:
 52 *   A. R. Edmonds, Angular Momentum in Quantum Mechanics, Princeton University Press 1974.
 53 *   E. Condon, and G. Shortley, The Theory of Atomic Spectra, Cambridge, 1935.
 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;   // maximal factorial n!  (2 x Lmax)
 70 /*static const double ARRAY_OVER  = 1.0e+300;   // force overflow */
 71 static const double nf_amc_log_fact[] = {0.0, 0.0, 0.69314718056, 1.79175946923, 3.17805383035, 4.78749174278, 6.57925121201, 8.52516136107, 10.6046029027, 12.8018274801, 15.1044125731, 17.5023078459, 19.9872144957, 22.5521638531, 25.1912211827, 27.8992713838, 30.6718601061, 33.5050734501, 36.395445208, 39.3398841872, 42.3356164608, 45.3801388985, 48.4711813518, 51.6066755678, 54.7847293981, 58.003605223, 61.261701761, 64.557538627, 67.8897431372, 71.2570389672, 74.6582363488, 78.0922235533, 81.5579594561, 85.0544670176, 88.5808275422, 92.1361756037, 95.7196945421, 99.3306124548, 102.968198615, 106.631760261, 110.320639715, 114.034211781, 117.7718814, 121.533081515, 125.317271149, 129.123933639, 132.952575036, 136.802722637, 140.673923648, 144.565743946, 148.477766952, 152.409592584, 156.360836303, 160.331128217, 164.320112263, 168.327445448, 172.352797139, 176.395848407, 180.456291418, 184.533828861, 188.628173424, 192.739047288, 196.866181673, 201.009316399, 205.168199483, 209.342586753, 213.532241495, 217.736934114, 221.956441819, 226.190548324, 230.439043566, 234.701723443, 238.978389562, 243.268849003, 247.572914096, 251.89040221, 256.22113555, 260.564940972, 264.921649799, 269.291097651, 273.673124286, 278.06757344, 282.474292688, 286.893133295, 291.323950094, 295.766601351, 300.220948647, 304.686856766, 309.16419358, 313.65282995, 318.15263962, 322.663499127, 327.185287704, 331.717887197, 336.261181979, 340.815058871, 345.379407062, 349.954118041, 354.539085519, 359.13420537, 363.739375556, 368.354496072, 372.979468886, 377.614197874, 382.258588773, 386.912549123, 391.575988217, 396.248817052, 400.930948279, 405.622296161, 410.322776527, 415.032306728, 419.7508056, 424.478193418, 429.214391867, 433.959323995, 438.712914186, 443.475088121, 448.245772745, 453.024896238, 457.812387981, 462.608178527, 467.412199572, 472.224383927, 477.044665493, 481.87297923, 486.709261137, 491.553448223, 496.405478487, 501.265290892, 506.132825342, 511.008022665, 515.890824588, 520.781173716, 525.679013516, 530.584288294, 535.49694318, 540.416924106, 545.344177791, 550.278651724, 555.220294147, 560.169054037, 565.124881095, 570.087725725, 575.057539025, 580.034272767, 585.017879389, 590.008311976, 595.005524249, 600.009470555, 605.020105849, 610.037385686, 615.061266207, 620.091704128, 625.128656731, 630.172081848, 635.221937855, 640.27818366, 645.340778693, 650.409682896, 655.484856711, 660.566261076, 665.653857411, 670.747607612, 675.84747404, 680.953419514, 686.065407302, 691.183401114, 696.307365094, 701.437263809, 706.573062246, 711.714725802, 716.862220279, 722.015511874, 727.174567173, 732.339353147, 737.509837142, 742.685986874, 747.867770425, 753.05515623, 758.248113081, 763.446610113, 768.6506168, 773.860102953, 779.07503871, 784.295394535, 789.521141209, 794.752249826, 799.988691789, 805.230438804, 810.477462876, 815.729736304, 820.987231676, 826.249921865, 831.517780024, 836.790779582, 842.068894242, 847.35209797, 852.640365001, 857.933669826, 863.231987192};
 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, int, int, int );
 81 static double cg3( int, int, int, int, int, int );
 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. INFINITY is return if n is negative or too large.
100 */
101     return G4Exp( nf_amc_log_factorial( n ) );
102 }
103 /*
104 ============================================================
105 */
106 double nf_amc_wigner_3j( int j1, int j2, int j3, int j4, int j5, int j6 ) {
107 /*
108 *       Wigner's 3J symbol (similar to Clebsh-Gordan)
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, j4, j5, j3 ) ) == 0.0 ) return ( 0.0 );
116     if( cg == INFINITY ) return( cg );
117     return( ( ( ( j1 - j2 - j6 ) % 4 == 0 ) ?  1.0 : -1.0 ) * cg / std::sqrt( j3 + 1.0 ) );          /* BRB j3 + 1 <= 0? */
118 }
119 /*
120 ============================================================
121 */
122 double nf_amc_wigner_6j( int j1, int j2, int j3, int j4, int j5, int j6 ) {
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; x[4] = j5; x[5] = j6;
131     for( i = 0; i < 6; i++ ) if ( x[i] == 0 ) return( w6j0( i, x ) );
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] != x[5] ) ) return( 0.0 );
142                x[5] = x[3]; x[0] = x[1]; x[3] = x[4];      break;
143        case 1: if ( ( x[0] != x[2] ) || ( x[3] != x[5] ) ) return( 0.0 );
144                x[5] = x[4];                                break;
145        case 2: if ( ( x[0] != x[1] ) || ( x[3] != x[4] ) ) return( 0.0 );
146                break;
147          //TK fix bug and add comment on 17-05-23
148                //This is the case of 6.3.2 of A. R. Edmonds, Angular Momentum in Quantum Mechanics, Princeton University Press 1974.
149        case 3: if ( ( x[1] != x[5] ) || ( x[2] != x[4] ) ) return( 0.0 );
150                x[5] = x[0]; x[0] = x[4]; x[3] = x[1];      break;
151        case 4: if ( ( x[0] != x[5] ) || ( x[2] != x[3] ) ) return( 0.0 );
152                x[5] = x[1];                                break;
153        case 5: if ( ( x[0] != x[4] ) || ( x[1] != x[3] ) ) return( 0.0 );
154                x[5] = x[2];                                break;
155     }
156 
157     if( ( x[5] > ( x[0] + x[3] ) ) || ( x[5] < std::abs( x[0] - x[3] ) ) ) return( 0.0 );
158     if( x[0] > MAX_FACTORIAL || x[3] > MAX_FACTORIAL ) {        /* BRB Why this test? Why not x[5]? */
159         return( INFINITY );
160     }
161 
162     return( 1.0 / std::sqrt( (double) ( ( x[0] + 1 ) * ( x[3] + 1 ) ) ) * ( ( ( x[0] + x[3] + x[5] ) / 2 ) % 2 != 0 ? -1 : 1 ) );
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, n2, n3, m1, m2, m3, x1, x2, x3, y[4];
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_log_fact[n2] + nf_amc_log_fact[n3] - nf_amc_log_fact[n+1];
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_fact[k1] - nf_amc_log_fact[l1] - nf_amc_log_fact[l2] - nf_amc_log_fact[l3] - nf_amc_log_fact[l4]
210                              - nf_amc_log_fact[m1] - nf_amc_log_fact[m2] - nf_amc_log_fact[m3] ) * ( ( k1 % 2 ) == 0 ? -1: 1 );
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 += k;
217 
218         for ( i = 0; i < k; i++ )
219             w6j = w - w6j * ( ( k2 - i ) * ( m1 + i ) * ( m2 + i ) * ( m3 + i ) )
220                     / ( ( l1 - i ) * ( l2 - i ) * ( l3 - i ) * ( l4 - i ) );
221     }
222     return( w6j );
223 }
224 /*
225 ============================================================
226 */
227 double nf_amc_wigner_9j( int j1, int j2, int j3, int j4, int j5, int j6, int j7, int j8, int j9 ) {
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( j2 - j6 ), std::abs( j4 - j8 ) );
239     i1 = min3(     ( j1 + j9 ),     ( j2 + j6 ),     ( j4 + j8 ) );
240 
241     rac = 0.0;
242     for ( i = i0; i <= i1; i += 2 ){ 
243         rac += nf_amc_racah( j1, j4, j9, j8, j7,  i )
244             *  nf_amc_racah( j2, j5,  i, j4, j8, j6 )
245             *  nf_amc_racah( j9,  i, j3, j2, j1, j6 ) * ( i + 1 );
246         if( rac == INFINITY ) return( INFINITY );
247     }
248 
249     return( ( ( (int)( ( j1 + j3 + j5 + j8 ) / 2 + j2 + j4 + j9 ) % 4 == 0 ) ?  1.0 : -1.0 ) * rac );
250 }
251 /*
252 ============================================================
253 */
254 double nf_amc_racah( int j1, int j2, int l2, int l1, int j3, int l3 ) {
255 /*
256 *       Racah coefficient definition in Edmonds (AR Edmonds, "Angular Momentum in Quantum Mechanics", Princeton (1980) is
257 *       W(j1, j2, l2, l1 ; j3, l3) = (-1)^(j1+j2+l1+l2) * { j1 j2 j3 }
258 *                                                         { l1 l2 l3 }
259 *       The call signature of W(...) appears jumbled, but hey, that's the convention.
260 *
261 *       This convention is exactly that used by Blatt-Biedenharn (Rev. Mod. Phys. 24, 258 (1952)) too
262 */
263 
264     double sig;
265 
266     sig = ( ( ( j1 + j2 + l1 + l2 ) % 4 == 0 ) ?  1.0 : -1.0 );
267     return sig * nf_amc_wigner_6j( j1, j2, j3, l1, l2, l3 );
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 ) return( 0.0 );
279    if ( ( j2 = (  a - b + c ) / 2 ) < 0 ) return( 0.0 );
280    if ( ( j3 = ( -a + b + c ) / 2 ) < 0 ) return( 0.0 );
281    j4 = ( a + b + c ) / 2 + 1;
282 
283    return( std::exp( 0.5 * ( nf_amc_log_fact[j1] + nf_amc_log_fact[j2] + nf_amc_log_fact[j3] - nf_amc_log_fact[j4] ) ) );
284 }
285 */
286 /*
287 ============================================================
288 */
289 double nf_amc_clebsh_gordan( int j1, int j2, int m1, int m2, int j3 ) {
290 /*
291 *       Clebsh-Gordan coefficient 
292 *           = <j1,j2,m1,m2|j3,m1+m2>
293 *           = (-)^(j1-j2+m1+m2) * std::sqrt(2*j3+1) * / j1 j2   j3   \
294 *                                                \ m1 m2 -m1-m2 /
295 *
296 *   Note: Last value m3 is preset to m1+m2.  Any other value will evaluate to 0.0.
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.0 );
303     if ( j1 + j2 + j3 > 2 * MAX_FACTORIAL ) return( INFINITY );
304 
305     m3 = m1 + m2;
306 
307     if ( ( x1 = ( j1 + m1 ) / 2 + 1 ) <= 0 ) return( 0.0 );
308     if ( ( x2 = ( j2 + m2 ) / 2 + 1 ) <= 0 ) return( 0.0 );
309     if ( ( x3 = ( j3 - m3 ) / 2 + 1 ) <= 0 ) return( 0.0 );
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( (double)j1 + 1.0 ) * ( ( y1 % 2 == 0 ) ? -1:1 ) );
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 - y2, x3 - 1, x1 + x2 - 2, x1 - y2, j1, j2, j3, m2 );
325         }
326         else if ( m2 == 0 && std::abs( m1 ) <=1 ){
327                           cg = cg2( x1 - y2 + y3, x2 - 1, x1 + x3 - 2, x3 - y1, j1, j3, j3, m1 );
328         }
329         else if ( m1 == 0 && std::abs( m3 ) <= 1 ){
330                           cg = cg2( x1, x1 - 1, x2 + x3 - 2, x2 - y3, j2, j3, j3, -m3 );
331         }
332         else              cg = cg3( x1, x2, x3, y1, y2, y3 );
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 ) != 0 ) return( 0.0 );
346     p2 = x1 + x2 - x3;
347     p3 =-x1 + x2 + x3;
348     p4 = x1 - x2 + x3;
349     if ( p2 <= 0 || p3 <= 0 || p4 <= 0 ) return( 0.0 );
350     if ( p1 >= MAX_FACTORIAL ) return( INFINITY );
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[q2] + nf_amc_log_fact[q3] + nf_amc_log_fact[q4] )
358         + 0.5 * ( nf_amc_log_fact[ 2 * x3 - 1 ] - nf_amc_log_fact[ 2 * x3 - 2 ]
359         + nf_amc_log_fact[p2] + nf_amc_log_fact[p3] + nf_amc_log_fact[p4] - nf_amc_log_fact[p1] );
360 
361     return( ( ( ( q1 + x1 - x2 ) % 2 == 0 ) ? 1.0 : -1.0 ) * G4Exp( a ) );
362 }
363 /*
364 ============================================================
365 */
366 static double cg2( int k, int q0, int z1, int z2, int w1, int w2, int w3, int mm ) {
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( 0.0 );
376     if ( p1 >= MAX_FACTORIAL ) return( INFINITY );
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_fact[ q2 ] + nf_amc_log_fact[ q3 ] + nf_amc_log_fact[ q4 ] )
384         + 0.5 * ( nf_amc_log_fact[ w3 + 1 ] - nf_amc_log_fact[ w3  ]
385                 + nf_amc_log_fact[ w1  ]    - nf_amc_log_fact[ w1 + 1 ]
386                 + nf_amc_log_fact[ w2  ]    - nf_amc_log_fact[ w2 + 1 ]
387                 + nf_amc_log_fact[ p2 ]     + nf_amc_log_fact[ p3 ] + nf_amc_log_fact[ p4 ] - nf_amc_log_fact[ p1 ] );
388 
389     return( ( ( ( q4 + k + ( mm > 0 ) * ( p1 + 2 ) ) % 2 == 0 ) ? -1.0 : 1.0 ) * 2.0 * G4Exp( a ) );
390 }
391 /*
392 ============================================================
393 */
394 static double cg3( int x1, int x2, int x3, int y1, int y2, int y3 ) {
395 
396     int nx, i, k1, k2, q1, q2, q3, q4, p1, p2, p3, z1, z2, z3;
397     double a, cg;
398 
399     nx = x1 + x2 + x3 - 1;
400     if ( ( z1 = nx - x1 - y1 ) < 0 ) return( 0.0 );
401     if ( ( z2 = nx - x2 - y2 ) < 0 ) return( 0.0 );
402     if ( ( z3 = nx - x3 - y3 ) < 0 ) return( 0.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[ x3 + y3 - 1 ] - nf_amc_log_fact[ x3 + y3 - 2 ] - nf_amc_log_fact[ nx - 1 ]
417                  + nf_amc_log_fact[ z1  ]    + nf_amc_log_fact[ z2  ]    + nf_amc_log_fact[ z3  ]
418                  + nf_amc_log_fact[ x1 - 1 ] + nf_amc_log_fact[ x2 - 1 ] + nf_amc_log_fact[ x3 - 1 ]
419                  + nf_amc_log_fact[ y1 - 1 ] + nf_amc_log_fact[ y2 - 1 ] + nf_amc_log_fact[ y3 - 1 ] )
420                  - nf_amc_log_fact[ p1  ]    - nf_amc_log_fact[ p2  ]    - nf_amc_log_fact[ p3  ]
421                  - nf_amc_log_fact[ q1  ]    - nf_amc_log_fact[ q3  ]    - nf_amc_log_fact[ q4  ] ) * ( ( ( q1 % 2 ) == 0 ) ? 1 : -1 );
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 + i ) * ( p3 + i ) ) / ( ( q2 - i ) * ( q3 - i ) * ( q4 - i ) );
432     }
433     return( cg );
434 }
435 /*
436 ============================================================
437 */
438 double nf_amc_z_coefficient( int l1, int j1, int l2, int j2, int s, int ll ) {
439 /*
440 *       Biedenharn's Z-coefficient coefficient
441 *           =  Z(l1  j1  l2  j2 | S L )
442 */
443     double z, clebsh_gordan = nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah = nf_amc_racah( l1, j1, l2, j2, s, ll );
444 
445     if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) ) return( INFINITY );
446     z = ( ( ( -l1 + l2 + ll ) % 8 == 0 ) ? 1.0 : -1.0 )
447         * std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
448 
449    return( z );
450 }
451 /*
452 ============================================================
453 */
454 double nf_amc_zbar_coefficient( int l1, int j1, int l2, int j2, int s, int ll ) {
455 /*
456 *       Lane & Thomas's Zbar-coefficient coefficient
457 *           = Zbar(l1  j1  l2  j2 | S L )
458 *           = (-i)^( -l1 + l2 + ll ) * Z(l1  j1  l2  j2 | S L )
459 *
460 *       Lane & Thomas Rev. Mod. Phys. 30, 257-353 (1958).
461 *       Note, Lane & Thomas define this because they did not like the different phase convention in Blatt & Biedenharn's Z coefficient.  They changed it to get better time-reversal behavior.
462 *       Froehner uses Lane & Thomas convention as does T. Kawano.
463 */
464     double zbar, clebsh_gordan = nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah = nf_amc_racah( l1, j1, l2, j2, s, ll );
465 
466     if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) ) return( INFINITY );
467     zbar = std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
468 
469    return( zbar );
470 }
471 /*
472 ============================================================
473 */
474 double nf_amc_reduced_matrix_element( int lt, int st, int jt, int l0, int j0, int l1, int j1 ) {
475 /*
476 *       Reduced Matrix Element for Tensor Operator
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 factor sqrt(2j1+1)
482 */
483     int    llt;
484     double x1, x2, x3, reduced_mat, clebsh_gordan;
485 
486     if ( parity( lt ) != parity( l0 ) * parity( l1 ) ) return( 0.0 );
487     if ( std::abs( l0 - l1 ) > lt || ( l0 + l1 ) < lt ) return( 0.0 );
488     if ( std::abs( ( j0 - j1 ) / 2 ) > jt || ( ( j0 + j1 ) / 2 ) < jt ) return( 0.0 );
489 
490     llt = 2 * lt;
491     jt *= 2;
492     st *= 2;
493 
494     if( ( clebsh_gordan = nf_amc_clebsh_gordan( j1, j0, 1, -1, jt ) ) == INFINITY ) return( INFINITY );
495 
496     reduced_mat = 1.0 / std::sqrt( 4 * M_PI ) * clebsh_gordan / std::sqrt( jt + 1.0 )         /* BRB jt + 1 <= 0? */
497                 * std::sqrt( ( j0 + 1.0 ) * ( j1 + 1.0 ) * ( llt + 1.0 ) )
498                 * parity( ( j1 - j0 ) / 2 ) * parity( ( -l0 + l1 + lt ) / 2 ) * parity( ( j0 - 1 ) / 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 ) / std::sqrt( lt * ( lt + 1.0 ) );
505         }
506         else if ( jt == ( llt - st ) ){
507             x3 = ( lt == 0 ) ? 0 : -( lt + x1 + x2 ) / std::sqrt( lt * ( 2.0 * lt + 1.0 ) );
508         }
509         else if ( jt == ( llt + st ) ){
510             x3 = ( lt + 1 - x1 - x2 ) / std::sqrt( ( 2.0 * lt + 1.0 ) * ( lt + 1.0 ) );
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