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