61 #if defined __cplusplus
70 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};
73 static int max3(
int a,
int b,
int c );
74 static int max4(
int a,
int b,
int c,
int d );
75 static int min3(
int a,
int b,
int c );
76 static double w6j0(
int,
int * );
77 static double w6j1(
int * );
78 static double cg1(
int,
int,
int );
79 static double cg2(
int,
int,
int,
int,
int,
int,
int,
int );
80 static double cg3(
int,
int,
int,
int,
int,
int );
90 if( n < 0 )
return( INFINITY );
113 if( ( j4 + j5 + j6 ) != 0 )
return( 0.0 );
115 if( cg == INFINITY )
return( cg );
116 return( ( ( ( j1 - j2 - j6 ) % 4 == 0 ) ? 1.0 : -1.0 ) * cg / std::sqrt( j3 + 1.0 ) );
129 x[0] = j1; x[1] = j2; x[2] = j3; x[3] = j4; x[4] = j5; x[5] = j6;
130 for( i = 0; i < 6; i++ )
if ( x[i] == 0 )
return(
w6j0( i, x ) );
137 static double w6j0(
int i,
int *
x ) {
140 case 0:
if ( ( x[1] != x[2] ) || ( x[4] != x[5] ) )
return( 0.0 );
141 x[5] = x[3]; x[0] = x[1]; x[3] = x[4];
break;
142 case 1:
if ( ( x[0] != x[2] ) || ( x[3] != x[5] ) )
return( 0.0 );
144 case 2:
if ( ( x[0] != x[1] ) || ( x[3] != x[4] ) )
return( 0.0 );
148 case 3:
if ( ( x[1] != x[5] ) || ( x[2] != x[4] ) )
return( 0.0 );
149 x[5] = x[0]; x[0] = x[4]; x[3] = x[1];
break;
150 case 4:
if ( ( x[0] != x[5] ) || ( x[2] != x[3] ) )
return( 0.0 );
152 case 5:
if ( ( x[0] != x[4] ) || ( x[1] != x[3] ) )
return( 0.0 );
156 if( ( x[5] > ( x[0] + x[3] ) ) || ( x[5] <
std::abs( x[0] - x[3] ) ) )
return( 0.0 );
161 return( 1.0 / std::sqrt( (
double) ( ( x[0] + 1 ) * ( x[3] + 1 ) ) ) * ( ( ( x[0] + x[3] + x[5] ) / 2 ) % 2 != 0 ? -1 : 1 ) );
169 int i,
k,
k1,
k2,
n, l1, l2, l3, l4,
n1, n2, n3, m1,
m2,
m3,
x1,
x2,
x3,
y[4];
170 static int a[3][4] = { { 0, 0, 3, 3},
176 for ( k = 0; k < 4; k++ ){
177 x1 = x[ ( a[0][
k] ) ];
178 x2 = x[ ( a[1][
k] ) ];
179 x3 = x[ ( a[2][
k] ) ];
181 n = ( x1 + x2 +
x3 ) / 2;
183 return( INFINITY ); }
188 if ( ( n1 = n - x3 ) < 0 )
return( 0.0 );
189 if ( ( n2 = n - x2 ) < 0 )
return( 0.0 );
190 if ( ( n3 = n - x1 ) < 0 )
return( 0.0 );
196 n1 = ( x[0] + x[1] + x[3] + x[4] ) / 2;
197 n2 = ( x[0] + x[2] + x[3] + x[5] ) / 2;
198 n3 = ( x[1] + x[2] + x[4] + x[5] ) / 2;
200 k1 =
max4( y[0], y[1], y[2], y[3] ) - 1;
201 k2 =
min3( n1, n2, n3 ) + 1;
203 l1 = k1 - y[0] + 1; m1 = n1 - k1 + 1;
204 l2 = k1 - y[1] + 1; m2 = n2 - k1 + 1;
205 l3 = k1 - y[2] + 1; m3 = n3 - k1 + 1;
210 if( w6j == INFINITY )
return( INFINITY );
214 m1 -= k-1; m2 -= k-1; m3 -= k-1;
215 l1 +=
k ; l2 +=
k ; l3 +=
k ; l4 +=
k;
217 for ( i = 0; i <
k; i++ )
218 w6j = w - w6j * ( ( k2 - i ) * ( m1 + i ) * ( m2 + i ) * ( m3 + i ) )
219 / ( ( l1 - i ) * ( l2 - i ) * ( l3 - i ) * ( l4 - i ) );
226 double nf_amc_wigner_9j(
int j1,
int j2,
int j3,
int j4,
int j5,
int j6,
int j7,
int j8,
int j9 ) {
238 i1 =
min3( ( j1 + j9 ), ( j2 + j6 ), ( j4 + j8 ) );
241 for ( i = i0; i <= i1; i += 2 ){
245 if( rac == INFINITY )
return( INFINITY );
248 return( ( ( (
int)( ( j1 + j3 + j5 + j8 ) / 2 + j2 + j4 + j9 ) % 4 == 0 ) ? 1.0 : -1.0 ) * rac );
253 double nf_amc_racah(
int j1,
int j2,
int l2,
int l1,
int j3,
int l3 ) {
265 sig = ( ( ( j1 + j2 + l1 + l2 ) % 4 == 0 ) ? 1.0 : -1.0 );
301 if ( j1 < 0 || j2 < 0 || j3 < 0)
return( 0.0 );
306 if ( ( x1 = ( j1 + m1 ) / 2 + 1 ) <= 0 )
return( 0.0 );
307 if ( ( x2 = ( j2 + m2 ) / 2 + 1 ) <= 0 )
return( 0.0 );
308 if ( ( x3 = ( j3 - m3 ) / 2 + 1 ) <= 0 )
return( 0.0 );
310 if ( ( y1 = x1 - m1 ) <= 0 )
return( 0.0 );
311 if ( ( y2 = x2 - m2 ) <= 0 )
return( 0.0 );
312 if ( ( y3 = x3 + m3 ) <= 0 )
return( 0.0 );
315 if ( j1 == j2 ) cg = ( 1.0 / std::sqrt( (
double)j1 + 1.0 ) * ( ( y1 % 2 == 0 ) ? -1:1 ) );
317 else if ( (j1 == 0 || j2 == 0 ) ){
318 if ( ( j1 + j2 ) == j3 ) cg = 1.0;
321 if( m3 == 0 &&
std::abs( m1 ) <= 1 ){
322 if( m1 == 0 ) cg =
cg1( x1, x2, x3 );
323 else cg =
cg2( x1 + y1 - y2, x3 - 1, x1 + x2 - 2, x1 - y2, j1, j2, j3, m2 );
325 else if ( m2 == 0 &&
std::abs( m1 ) <=1 ){
326 cg =
cg2( x1 - y2 + y3, x2 - 1, x1 + x3 - 2, x3 - y1, j1, j3, j3, m1 );
328 else if ( m1 == 0 &&
std::abs( m3 ) <= 1 ){
329 cg =
cg2( x1, x1 - 1, x2 + x3 - 2, x2 - y3, j2, j3, j3, -m3 );
331 else cg =
cg3( x1, x2, x3, y1, y2, y3 );
341 int p1, p2, p3, p4, q1, q2, q3, q4;
344 p1 = x1 + x2 + x3 - 1;
if ( ( p1 % 2 ) != 0 )
return( 0.0 );
348 if ( p2 <= 0 || p3 <= 0 || p4 <= 0 )
return( 0.0 );
351 q1 = ( p1 + 1 ) / 2 - 1; p1--;
352 q2 = ( p2 + 1 ) / 2 - 1; p2--;
353 q3 = ( p3 + 1 ) / 2 - 1; p3--;
354 q4 = ( p4 + 1 ) / 2 - 1; p4--;
360 return( ( ( ( q1 + x1 - x2 ) % 2 == 0 ) ? 1.0 : -1.0 ) *
G4Exp( a ) );
365 static double cg2(
int k,
int q0,
int z1,
int z2,
int w1,
int w2,
int w3,
int mm ) {
367 int q1, q2, q3, q4, p1, p2, p3, p4;
374 if ( p2 <= 0 || p3 <= 0 || p4 <= 0)
return( 0.0 );
377 q1 = ( p1 + 1 ) / 2 - 1; p1--;
378 q2 = ( p2 + 1 ) / 2 - 1; p2--;
379 q3 = ( p3 + 1 ) / 2 - 1; p3--;
380 q4 = ( p4 + 1 ) / 2 - 1; p4--;
388 return( ( ( ( q4 + k + ( mm > 0 ) * ( p1 + 2 ) ) % 2 == 0 ) ? -1.0 : 1.0 ) * 2.0 *
G4Exp( a ) );
395 int nx, i,
k1,
k2, q1, q2, q3, q4, p1, p2, p3,
z1,
z2, z3;
398 nx = x1 + x2 + x3 - 1;
399 if ( ( z1 = nx - x1 - y1 ) < 0 )
return( 0.0 );
400 if ( ( z2 = nx - x2 - y2 ) < 0 )
return( 0.0 );
401 if ( ( z3 = nx - x3 - y3 ) < 0 )
return( 0.0 );
406 q1 =
max3( k1, k2, 0 );
407 q2 =
min3( y1, x2, z3 + 1 ) - 1;
421 if( cg == INFINITY )
return( INFINITY );
429 for( i = 0; i < ( q2 - q1 ); i++ )
430 cg = a - cg * ( ( p1 + i ) * ( p2 + i ) * ( p3 + i ) ) / ( ( q2 - i ) * ( q3 - i ) * ( q4 - i ) );
442 double z, clebsh_gordan =
nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah =
nf_amc_racah( l1, j1, l2, j2, s, ll );
444 if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) )
return( INFINITY );
445 z = ( ( ( -l1 + l2 + ll ) % 8 == 0 ) ? 1.0 : -1.0 )
446 * std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
463 double zbar, clebsh_gordan =
nf_amc_clebsh_gordan( l1, l2, 0, 0, ll ), racah =
nf_amc_racah( l1, j1, l2, j2, s, ll );
465 if( ( clebsh_gordan == INFINITY ) || ( racah == INFINITY ) )
return( INFINITY );
466 zbar = std::sqrt( l1 + 1.0 ) * std::sqrt( l2 + 1.0 ) * std::sqrt( j1 + 1.0 ) * std::sqrt( j2 + 1.0 ) * clebsh_gordan * racah;
483 double x1,
x2,
x3, reduced_mat, clebsh_gordan;
486 if (
std::abs( l0 - l1 ) > lt || ( l0 + l1 ) < lt )
return( 0.0 );
487 if (
std::abs( ( j0 - j1 ) / 2 ) > jt || ( ( j0 + j1 ) / 2 ) < jt )
return( 0.0 );
493 if( ( clebsh_gordan =
nf_amc_clebsh_gordan( j1, j0, 1, -1, jt ) ) == INFINITY )
return( INFINITY );
495 reduced_mat = 1.0 / std::sqrt( 4 *
M_PI ) * clebsh_gordan / std::sqrt( jt + 1.0 )
496 * std::sqrt( ( j0 + 1.0 ) * ( j1 + 1.0 ) * ( llt + 1.0 ) )
497 *
parity( ( j1 - j0 ) / 2 ) *
parity( ( -l0 + l1 + lt ) / 2 ) *
parity( ( j0 - 1 ) / 2 );
500 x1 = ( l0 - j0 / 2.0 ) * ( j0 + 1.0 );
501 x2 = ( l1 - j1 / 2.0 ) * ( j1 + 1.0 );
503 x3 = ( lt == 0 ) ? 0 : ( x1 - x2 ) / std::sqrt( lt * ( lt + 1.0 ) );
505 else if ( jt == ( llt - st ) ){
506 x3 = ( lt == 0 ) ? 0 : -( lt + x1 + x2 ) / std::sqrt( lt * ( 2.0 * lt + 1.0 ) );
508 else if ( jt == ( llt + st ) ){
509 x3 = ( lt + 1 - x1 -
x2 ) / std::sqrt( ( 2.0 * lt + 1.0 ) * ( lt + 1.0 ) );
518 return( reduced_mat );
525 return( ( ( x / 2 ) % 2 == 0 ) ? 1 : -1 );
556 #if defined __cplusplus