| 1 |
chuckv |
2869 |
subroutine WIGNER3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, & |
| 2 |
|
|
IER) |
| 3 |
|
|
use status |
| 4 |
|
|
! |
| 5 |
|
|
!! DRC3JM evaluates the 3j symbol g(M2) for all allowed values of M2. |
| 6 |
|
|
! |
| 7 |
|
|
!***PURPOSE Evaluate the 3j symbol g(M2) = (L1 L2 L3 ) |
| 8 |
|
|
! (M1 M2 -M1-M2) |
| 9 |
|
|
! for all allowed values of M2, the other parameters |
| 10 |
|
|
! being held fixed. |
| 11 |
|
|
! |
| 12 |
|
|
!***LIBRARY SLATEC |
| 13 |
|
|
!***CATEGORY C19 |
| 14 |
|
|
!***TYPE DOUBLE PRECISION (RC3JM-S, DRC3JM-D) |
| 15 |
|
|
!***KEYWORDS 3J COEFFICIENTS, 3J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS, |
| 16 |
|
|
! RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS, |
| 17 |
|
|
! WIGNER COEFFICIENTS |
| 18 |
|
|
!***AUTHOR Gordon, R. G., Harvard University |
| 19 |
|
|
! Schulten, K., Max Planck Institute |
| 20 |
|
|
!***DESCRIPTION |
| 21 |
|
|
! |
| 22 |
|
|
! *Usage: |
| 23 |
|
|
! |
| 24 |
|
|
! DOUBLE PRECISION L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM) |
| 25 |
|
|
! INTEGER NDIM, IER |
| 26 |
|
|
! |
| 27 |
|
|
! call DRC3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, IER) |
| 28 |
|
|
! |
| 29 |
|
|
! *Arguments: |
| 30 |
|
|
! |
| 31 |
|
|
! L1 :IN Parameter in 3j symbol. |
| 32 |
|
|
! |
| 33 |
|
|
! L2 :IN Parameter in 3j symbol. |
| 34 |
|
|
! |
| 35 |
|
|
! L3 :IN Parameter in 3j symbol. |
| 36 |
|
|
! |
| 37 |
|
|
! M1 :IN Parameter in 3j symbol. |
| 38 |
|
|
! |
| 39 |
|
|
! M2MIN :OUT Smallest allowable M2 in 3j symbol. |
| 40 |
|
|
! |
| 41 |
|
|
! M2MAX :OUT Largest allowable M2 in 3j symbol. |
| 42 |
|
|
! |
| 43 |
|
|
! THRCOF :OUT Set of 3j coefficients generated by evaluating the |
| 44 |
|
|
! 3j symbol for all allowed values of M2. THRCOF(I) |
| 45 |
|
|
! will contain g(M2MIN+I-1), I=1,2,...,M2MAX-M2MIN+1. |
| 46 |
|
|
! |
| 47 |
|
|
! NDIM :IN Declared length of THRCOF in calling program. |
| 48 |
|
|
! |
| 49 |
|
|
! IER :OUT Error flag. |
| 50 |
|
|
! IER=0 No errors. |
| 51 |
|
|
! IER=1 Either L1 < ABS(M1) or L1+ABS(M1) non-integer. |
| 52 |
|
|
! IER=2 ABS(L1-L2) <= L3 <= L1+L2 not satisfied. |
| 53 |
|
|
! IER=3 L1+L2+L3 not an integer. |
| 54 |
|
|
! IER=4 M2MAX-M2MIN not an integer. |
| 55 |
|
|
! IER=5 M2MAX less than M2MIN. |
| 56 |
|
|
! IER=6 NDIM less than M2MAX-M2MIN+1. |
| 57 |
|
|
! |
| 58 |
|
|
! *Description: |
| 59 |
|
|
! |
| 60 |
|
|
! Although conventionally the parameters of the vector addition |
| 61 |
|
|
! coefficients satisfy certain restrictions, such as being integers |
| 62 |
|
|
! or integers plus 1/2, the restrictions imposed on input to this |
| 63 |
|
|
! subroutine are somewhat weaker. See, for example, Section 27.9 of |
| 64 |
|
|
! Abramowitz and Stegun or Appendix C of Volume II of A. Messiah. |
| 65 |
|
|
! The restrictions imposed by this subroutine are |
| 66 |
|
|
! 1. L1 >= ABS(M1) and L1+ABS(M1) must be an integer; |
| 67 |
|
|
! 2. ABS(L1-L2) <= L3 <= L1+L2; |
| 68 |
|
|
! 3. L1+L2+L3 must be an integer; |
| 69 |
|
|
! 4. M2MAX-M2MIN must be an integer, where |
| 70 |
|
|
! M2MAX=MIN(L2,L3-M1) and M2MIN=MAX(-L2,-L3-M1). |
| 71 |
|
|
! If the conventional restrictions are satisfied, then these |
| 72 |
|
|
! restrictions are met. |
| 73 |
|
|
! |
| 74 |
|
|
! The user should be cautious in using input parameters that do |
| 75 |
|
|
! not satisfy the conventional restrictions. For example, the |
| 76 |
|
|
! the subroutine produces values of |
| 77 |
|
|
! g(M2) = (0.751.50 1.75 ) |
| 78 |
|
|
! (0.25 M2 -0.25-M2) |
| 79 |
|
|
! for M2=-1.5,-0.5,0.5,1.5 but none of the symmetry properties of the |
| 80 |
|
|
! 3j symbol, set forth on page 1056 of Messiah, is satisfied. |
| 81 |
|
|
! |
| 82 |
|
|
! The subroutine generates g(M2MIN), g(M2MIN+1), ..., g(M2MAX) |
| 83 |
|
|
! where M2MIN and M2MAX are defined above. The sequence g(M2) is |
| 84 |
|
|
! generated by a three-term recurrence algorithm with scaling to |
| 85 |
|
|
! control overflow. Both backward and forward recurrence are used to |
| 86 |
|
|
! maintain numerical stability. The two recurrence sequences are |
| 87 |
|
|
! matched at an interior point and are normalized from the unitary |
| 88 |
|
|
! property of 3j coefficients and Wigner's phase convention. |
| 89 |
|
|
! |
| 90 |
|
|
! The algorithm is suited to applications in which large quantum |
| 91 |
|
|
! numbers arise, such as in molecular dynamics. |
| 92 |
|
|
! |
| 93 |
|
|
!***REFERENCES 1. Abramowitz, M., and Stegun, I. A., Eds., Handbook |
| 94 |
|
|
! of Mathematical Functions with Formulas, Graphs |
| 95 |
|
|
! and Mathematical Tables, NBS Applied Mathematics |
| 96 |
|
|
! Series 55, June 1964 and subsequent printings. |
| 97 |
|
|
! 2. Messiah, Albert., Quantum Mechanics, Volume II, |
| 98 |
|
|
! North-Holland Publishing Company, 1963. |
| 99 |
|
|
! 3. Schulten, Klaus and Gordon, Roy G., Exact recursive |
| 100 |
|
|
! evaluation of 3j and 6j coefficients for quantum- |
| 101 |
|
|
! mechanical coupling of angular momenta, J Math |
| 102 |
|
|
! Phys, v 16, no. 10, October 1975, pp. 1961-1970. |
| 103 |
|
|
! 4. Schulten, Klaus and Gordon, Roy G., Semiclassical |
| 104 |
|
|
! approximations to 3j and 6j coefficients for |
| 105 |
|
|
! quantum-mechanical coupling of angular momenta, |
| 106 |
|
|
! J Math Phys, v 16, no. 10, October 1975, |
| 107 |
|
|
! pp. 1971-1988. |
| 108 |
|
|
! 5. Schulten, Klaus and Gordon, Roy G., Recursive |
| 109 |
|
|
! evaluation of 3j and 6j coefficients, Computer |
| 110 |
|
|
! Phys Comm, v 11, 1976, pp. 269-278. |
| 111 |
|
|
!***ROUTINES CALLED D1MACH, XERMSG |
| 112 |
|
|
!***REVISION HISTORY (YYMMDD) |
| 113 |
|
|
! 750101 DATE WRITTEN |
| 114 |
|
|
! 880515 SLATEC prologue added by G. C. Nielson, NBS; parameters |
| 115 |
|
|
! HUGE and TINY revised to depend on D1MACH. |
| 116 |
|
|
! 891229 Prologue description rewritten; other prologue sections |
| 117 |
|
|
! revised; MMATCH (location of match point for recurrences) |
| 118 |
|
|
! removed from argument list; argument IER changed to serve |
| 119 |
|
|
! only as an error flag (previously, in cases without error, |
| 120 |
|
|
! it returned the number of scalings); number of error codes |
| 121 |
|
|
! increased to provide more precise error information; |
| 122 |
|
|
! program comments revised; SLATEC error handler calls |
| 123 |
|
|
! introduced to enable printing of error messages to meet |
| 124 |
|
|
! SLATEC standards. These changes were done by D. W. Lozier, |
| 125 |
|
|
! M. A. McClain and J. M. Smith of the National Institute |
| 126 |
|
|
! of Standards and Technology, formerly NBS. |
| 127 |
|
|
! 910415 Mixed type expressions eliminated; variable C1 initialized; |
| 128 |
|
|
! description of THRCOF expanded. These changes were done by |
| 129 |
|
|
! D. W. Lozier. |
| 130 |
|
|
!***END PROLOGUE DRC3JM |
| 131 |
|
|
! |
| 132 |
|
|
INTEGER NDIM, IER |
| 133 |
|
|
DOUBLE PRECISION L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM) |
| 134 |
|
|
! |
| 135 |
|
|
INTEGER I, INDEX, LSTEP, N, NFIN, NFINP1, NFINP2, NFINP3, NLIM, & |
| 136 |
|
|
NSTEP2 |
| 137 |
|
|
DOUBLE PRECISION A1, A1S, C1, C1OLD, C2, CNORM, D1MACH, DV, EPS, & |
| 138 |
|
|
HUGE, M2, M3, NEWFAC, OLDFAC, ONE, RATIO, SIGN1, & |
| 139 |
|
|
SIGN2, SRHUGE, SRTINY, SUM1, SUM2, SUMBAC, & |
| 140 |
|
|
SUMFOR, SUMUNI, THRESH, TINY, TWO, X, X1, X2, X3, & |
| 141 |
|
|
Y, Y1, Y2, Y3, ZERO |
| 142 |
|
|
! |
| 143 |
|
|
DATA ZERO,EPS,ONE,TWO /0.0D0,0.01D0,1.0D0,2.0D0/ |
| 144 |
|
|
! |
| 145 |
|
|
!***FIRST EXECUTABLE STATEMENT DRC3JM |
| 146 |
|
|
IER=0 |
| 147 |
|
|
! HUGE is the square root of one twentieth of the largest floating |
| 148 |
|
|
! point number, approximately. |
| 149 |
|
|
HUGE = SQRT(D1MACH(2)/20.0D0) |
| 150 |
|
|
SRHUGE = SQRT(HUGE) |
| 151 |
|
|
TINY = 1.0D0/HUGE |
| 152 |
|
|
SRTINY = 1.0D0/SRHUGE |
| 153 |
|
|
! |
| 154 |
|
|
! MMATCH = ZERO |
| 155 |
|
|
! |
| 156 |
|
|
! |
| 157 |
|
|
! Check error conditions 1, 2, and 3. |
| 158 |
|
|
if ( (L1-ABS(M1)+EPS < ZERO).OR. & |
| 159 |
|
|
(MOD(L1+ABS(M1)+EPS,ONE) >= EPS+EPS))THEN |
| 160 |
|
|
IER=1 |
| 161 |
|
|
call handleError('Wigner3jm','L1-ABS(M1) less than zero or '// & |
| 162 |
|
|
'L1+ABS(M1) not integer.') |
| 163 |
|
|
return |
| 164 |
|
|
ELSEIF((L1+L2-L3 < -EPS).OR.(L1-L2+L3 < -EPS).OR. & |
| 165 |
|
|
(-L1+L2+L3 < -EPS))THEN |
| 166 |
|
|
IER=2 |
| 167 |
|
|
call handleError('Wigner3jm','L1, L2, L3 do not satisfy '// & |
| 168 |
|
|
'triangular condition.') |
| 169 |
|
|
return |
| 170 |
|
|
ELSEIF(MOD(L1+L2+L3+EPS,ONE) >= EPS+EPS)THEN |
| 171 |
|
|
IER=3 |
| 172 |
|
|
call handleError('Wigner3jm','L1+L2+L3 not integer.') |
| 173 |
|
|
return |
| 174 |
|
|
end if |
| 175 |
|
|
! |
| 176 |
|
|
! |
| 177 |
|
|
! Limits for M2 |
| 178 |
|
|
M2MIN = MAX(-L2,-L3-M1) |
| 179 |
|
|
M2MAX = MIN(L2,L3-M1) |
| 180 |
|
|
! |
| 181 |
|
|
! Check error condition 4. |
| 182 |
|
|
if ( MOD(M2MAX-M2MIN+EPS,ONE) >= EPS+EPS)THEN |
| 183 |
|
|
IER=4 |
| 184 |
|
|
call handleError('Wigner3jm','M2MAX-M2MIN not integer.') |
| 185 |
|
|
return |
| 186 |
|
|
end if |
| 187 |
|
|
if ( M2MIN < M2MAX-EPS) go to 20 |
| 188 |
|
|
if ( M2MIN < M2MAX+EPS) go to 10 |
| 189 |
|
|
! |
| 190 |
|
|
! Check error condition 5. |
| 191 |
|
|
IER=5 |
| 192 |
|
|
call handleError('Wigner3jm','M2MIN greater than M2MAX.') |
| 193 |
|
|
return |
| 194 |
|
|
! |
| 195 |
|
|
! |
| 196 |
|
|
! This is reached in case that M2 and M3 can take only one value. |
| 197 |
|
|
10 CONTINUE |
| 198 |
|
|
! MSCALE = 0 |
| 199 |
|
|
THRCOF(1) = (-ONE) ** INT(ABS(L2-L3-M1)+EPS) / & |
| 200 |
|
|
SQRT(L1+L2+L3+ONE) |
| 201 |
|
|
return |
| 202 |
|
|
! |
| 203 |
|
|
! This is reached in case that M1 and M2 take more than one value. |
| 204 |
|
|
20 CONTINUE |
| 205 |
|
|
! MSCALE = 0 |
| 206 |
|
|
NFIN = INT(M2MAX-M2MIN+ONE+EPS) |
| 207 |
|
|
if ( NDIM-NFIN) 21, 23, 23 |
| 208 |
|
|
! |
| 209 |
|
|
! Check error condition 6. |
| 210 |
|
|
21 IER = 6 |
| 211 |
|
|
call handleError('Wigner3jm','Dimension of result array for '// & |
| 212 |
|
|
'3j coefficients too small.') |
| 213 |
|
|
return |
| 214 |
|
|
! |
| 215 |
|
|
! |
| 216 |
|
|
! |
| 217 |
|
|
! Start of forward recursion from M2 = M2MIN |
| 218 |
|
|
! |
| 219 |
|
|
23 M2 = M2MIN |
| 220 |
|
|
THRCOF(1) = SRTINY |
| 221 |
|
|
NEWFAC = 0.0D0 |
| 222 |
|
|
C1 = 0.0D0 |
| 223 |
|
|
SUM1 = TINY |
| 224 |
|
|
! |
| 225 |
|
|
! |
| 226 |
|
|
LSTEP = 1 |
| 227 |
|
|
30 LSTEP = LSTEP + 1 |
| 228 |
|
|
M2 = M2 + ONE |
| 229 |
|
|
M3 = - M1 - M2 |
| 230 |
|
|
! |
| 231 |
|
|
! |
| 232 |
|
|
OLDFAC = NEWFAC |
| 233 |
|
|
A1 = (L2-M2+ONE) * (L2+M2) * (L3+M3+ONE) * (L3-M3) |
| 234 |
|
|
NEWFAC = SQRT(A1) |
| 235 |
|
|
! |
| 236 |
|
|
! |
| 237 |
|
|
DV = (L1+L2+L3+ONE)*(L2+L3-L1) - (L2-M2+ONE)*(L3+M3+ONE) & |
| 238 |
|
|
- (L2+M2-ONE)*(L3-M3-ONE) |
| 239 |
|
|
! |
| 240 |
|
|
if ( LSTEP-2) 32, 32, 31 |
| 241 |
|
|
! |
| 242 |
|
|
31 C1OLD = ABS(C1) |
| 243 |
|
|
32 C1 = - DV / NEWFAC |
| 244 |
|
|
! |
| 245 |
|
|
if ( LSTEP > 2) go to 60 |
| 246 |
|
|
! |
| 247 |
|
|
! |
| 248 |
|
|
! If M2 = M2MIN + 1, the third term in the recursion equation vanishes, |
| 249 |
|
|
! hence |
| 250 |
|
|
! |
| 251 |
|
|
X = SRTINY * C1 |
| 252 |
|
|
THRCOF(2) = X |
| 253 |
|
|
SUM1 = SUM1 + TINY * C1*C1 |
| 254 |
|
|
if ( LSTEP == NFIN) go to 220 |
| 255 |
|
|
go to 30 |
| 256 |
|
|
! |
| 257 |
|
|
! |
| 258 |
|
|
60 C2 = - OLDFAC / NEWFAC |
| 259 |
|
|
! |
| 260 |
|
|
! Recursion to the next 3j coefficient |
| 261 |
|
|
X = C1 * THRCOF(LSTEP-1) + C2 * THRCOF(LSTEP-2) |
| 262 |
|
|
THRCOF(LSTEP) = X |
| 263 |
|
|
SUMFOR = SUM1 |
| 264 |
|
|
SUM1 = SUM1 + X*X |
| 265 |
|
|
if ( LSTEP == NFIN) go to 100 |
| 266 |
|
|
! |
| 267 |
|
|
! See if last unnormalized 3j coefficient exceeds SRHUGE |
| 268 |
|
|
! |
| 269 |
|
|
if ( ABS(X) < SRHUGE) go to 80 |
| 270 |
|
|
! |
| 271 |
|
|
! This is reached if last 3j coefficient larger than SRHUGE, |
| 272 |
|
|
! so that the recursion series THRCOF(1), ... , THRCOF(LSTEP) |
| 273 |
|
|
! has to be rescaled to prevent overflow |
| 274 |
|
|
! |
| 275 |
|
|
! MSCALE = MSCALE + 1 |
| 276 |
|
|
DO 70 I=1,LSTEP |
| 277 |
|
|
if ( ABS(THRCOF(I)) < SRTINY) THRCOF(I) = ZERO |
| 278 |
|
|
70 THRCOF(I) = THRCOF(I) / SRHUGE |
| 279 |
|
|
SUM1 = SUM1 / HUGE |
| 280 |
|
|
SUMFOR = SUMFOR / HUGE |
| 281 |
|
|
X = X / SRHUGE |
| 282 |
|
|
! |
| 283 |
|
|
! |
| 284 |
|
|
! As long as ABS(C1) is decreasing, the recursion proceeds towards |
| 285 |
|
|
! increasing 3j values and, hence, is numerically stable. Once |
| 286 |
|
|
! an increase of ABS(C1) is detected, the recursion direction is |
| 287 |
|
|
! reversed. |
| 288 |
|
|
! |
| 289 |
|
|
80 if ( C1OLD-ABS(C1)) 100, 100, 30 |
| 290 |
|
|
! |
| 291 |
|
|
! |
| 292 |
|
|
! Keep three 3j coefficients around MMATCH for comparison later |
| 293 |
|
|
! with backward recursion values. |
| 294 |
|
|
! |
| 295 |
|
|
100 CONTINUE |
| 296 |
|
|
! MMATCH = M2 - 1 |
| 297 |
|
|
NSTEP2 = NFIN - LSTEP + 3 |
| 298 |
|
|
X1 = X |
| 299 |
|
|
X2 = THRCOF(LSTEP-1) |
| 300 |
|
|
X3 = THRCOF(LSTEP-2) |
| 301 |
|
|
! |
| 302 |
|
|
! Starting backward recursion from M2MAX taking NSTEP2 steps, so |
| 303 |
|
|
! that forwards and backwards recursion overlap at the three points |
| 304 |
|
|
! M2 = MMATCH+1, MMATCH, MMATCH-1. |
| 305 |
|
|
! |
| 306 |
|
|
NFINP1 = NFIN + 1 |
| 307 |
|
|
NFINP2 = NFIN + 2 |
| 308 |
|
|
NFINP3 = NFIN + 3 |
| 309 |
|
|
THRCOF(NFIN) = SRTINY |
| 310 |
|
|
SUM2 = TINY |
| 311 |
|
|
! |
| 312 |
|
|
! |
| 313 |
|
|
! |
| 314 |
|
|
M2 = M2MAX + TWO |
| 315 |
|
|
LSTEP = 1 |
| 316 |
|
|
110 LSTEP = LSTEP + 1 |
| 317 |
|
|
M2 = M2 - ONE |
| 318 |
|
|
M3 = - M1 - M2 |
| 319 |
|
|
OLDFAC = NEWFAC |
| 320 |
|
|
A1S = (L2-M2+TWO) * (L2+M2-ONE) * (L3+M3+TWO) * (L3-M3-ONE) |
| 321 |
|
|
NEWFAC = SQRT(A1S) |
| 322 |
|
|
DV = (L1+L2+L3+ONE)*(L2+L3-L1) - (L2-M2+ONE)*(L3+M3+ONE) & |
| 323 |
|
|
- (L2+M2-ONE)*(L3-M3-ONE) |
| 324 |
|
|
C1 = - DV / NEWFAC |
| 325 |
|
|
if ( LSTEP > 2) go to 120 |
| 326 |
|
|
! |
| 327 |
|
|
! If M2 = M2MAX + 1 the third term in the recursion equation vanishes |
| 328 |
|
|
! |
| 329 |
|
|
Y = SRTINY * C1 |
| 330 |
|
|
THRCOF(NFIN-1) = Y |
| 331 |
|
|
if ( LSTEP == NSTEP2) go to 200 |
| 332 |
|
|
SUMBAC = SUM2 |
| 333 |
|
|
SUM2 = SUM2 + Y*Y |
| 334 |
|
|
go to 110 |
| 335 |
|
|
! |
| 336 |
|
|
120 C2 = - OLDFAC / NEWFAC |
| 337 |
|
|
! |
| 338 |
|
|
! Recursion to the next 3j coefficient |
| 339 |
|
|
! |
| 340 |
|
|
Y = C1 * THRCOF(NFINP2-LSTEP) + C2 * THRCOF(NFINP3-LSTEP) |
| 341 |
|
|
! |
| 342 |
|
|
if ( LSTEP == NSTEP2) go to 200 |
| 343 |
|
|
! |
| 344 |
|
|
THRCOF(NFINP1-LSTEP) = Y |
| 345 |
|
|
SUMBAC = SUM2 |
| 346 |
|
|
SUM2 = SUM2 + Y*Y |
| 347 |
|
|
! |
| 348 |
|
|
! |
| 349 |
|
|
! See if last 3j coefficient exceeds SRHUGE |
| 350 |
|
|
! |
| 351 |
|
|
if ( ABS(Y) < SRHUGE) go to 110 |
| 352 |
|
|
! |
| 353 |
|
|
! This is reached if last 3j coefficient larger than SRHUGE, |
| 354 |
|
|
! so that the recursion series THRCOF(NFIN), ... , THRCOF(NFIN-LSTEP+1) |
| 355 |
|
|
! has to be rescaled to prevent overflow. |
| 356 |
|
|
! |
| 357 |
|
|
! MSCALE = MSCALE + 1 |
| 358 |
|
|
DO 111 I=1,LSTEP |
| 359 |
|
|
INDEX = NFIN - I + 1 |
| 360 |
|
|
if ( ABS(THRCOF(INDEX)) < SRTINY) & |
| 361 |
|
|
THRCOF(INDEX) = ZERO |
| 362 |
|
|
111 THRCOF(INDEX) = THRCOF(INDEX) / SRHUGE |
| 363 |
|
|
SUM2 = SUM2 / HUGE |
| 364 |
|
|
SUMBAC = SUMBAC / HUGE |
| 365 |
|
|
! |
| 366 |
|
|
go to 110 |
| 367 |
|
|
! |
| 368 |
|
|
! |
| 369 |
|
|
! |
| 370 |
|
|
! The forward recursion 3j coefficients X1, X2, X3 are to be matched |
| 371 |
|
|
! with the corresponding backward recursion values Y1, Y2, Y3. |
| 372 |
|
|
! |
| 373 |
|
|
200 Y3 = Y |
| 374 |
|
|
Y2 = THRCOF(NFINP2-LSTEP) |
| 375 |
|
|
Y1 = THRCOF(NFINP3-LSTEP) |
| 376 |
|
|
! |
| 377 |
|
|
! |
| 378 |
|
|
! Determine now RATIO such that YI = RATIO * XI (I=1,2,3) holds |
| 379 |
|
|
! with minimal error. |
| 380 |
|
|
! |
| 381 |
|
|
RATIO = ( X1*Y1 + X2*Y2 + X3*Y3 ) / ( X1*X1 + X2*X2 + X3*X3 ) |
| 382 |
|
|
NLIM = NFIN - NSTEP2 + 1 |
| 383 |
|
|
! |
| 384 |
|
|
if ( ABS(RATIO) < ONE) go to 211 |
| 385 |
|
|
! |
| 386 |
|
|
DO 210 N=1,NLIM |
| 387 |
|
|
210 THRCOF(N) = RATIO * THRCOF(N) |
| 388 |
|
|
SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC |
| 389 |
|
|
go to 230 |
| 390 |
|
|
! |
| 391 |
|
|
211 NLIM = NLIM + 1 |
| 392 |
|
|
RATIO = ONE / RATIO |
| 393 |
|
|
DO 212 N=NLIM,NFIN |
| 394 |
|
|
212 THRCOF(N) = RATIO * THRCOF(N) |
| 395 |
|
|
SUMUNI = SUMFOR + RATIO*RATIO*SUMBAC |
| 396 |
|
|
go to 230 |
| 397 |
|
|
! |
| 398 |
|
|
220 SUMUNI = SUM1 |
| 399 |
|
|
! |
| 400 |
|
|
! |
| 401 |
|
|
! Normalize 3j coefficients |
| 402 |
|
|
! |
| 403 |
|
|
230 CNORM = ONE / SQRT((L1+L1+ONE) * SUMUNI) |
| 404 |
|
|
! |
| 405 |
|
|
! Sign convention for last 3j coefficient determines overall phase |
| 406 |
|
|
! |
| 407 |
|
|
SIGN1 = SIGN(ONE,THRCOF(NFIN)) |
| 408 |
|
|
SIGN2 = (-ONE) ** INT(ABS(L2-L3-M1)+EPS) |
| 409 |
|
|
if ( SIGN1*SIGN2) 235,235,236 |
| 410 |
|
|
235 CNORM = - CNORM |
| 411 |
|
|
! |
| 412 |
|
|
236 if ( ABS(CNORM) < ONE) go to 250 |
| 413 |
|
|
! |
| 414 |
|
|
DO 240 N=1,NFIN |
| 415 |
|
|
240 THRCOF(N) = CNORM * THRCOF(N) |
| 416 |
|
|
return |
| 417 |
|
|
! |
| 418 |
|
|
250 THRESH = TINY / ABS(CNORM) |
| 419 |
|
|
DO 251 N=1,NFIN |
| 420 |
|
|
if ( ABS(THRCOF(N)) < THRESH) THRCOF(N) = ZERO |
| 421 |
|
|
251 THRCOF(N) = CNORM * THRCOF(N) |
| 422 |
|
|
! |
| 423 |
|
|
! |
| 424 |
|
|
! |
| 425 |
|
|
return |
| 426 |
|
|
end |
