| 1 |
gezelter |
1600 |
/* |
| 2 |
|
|
* Matpack Wigner3jm special function imported and modified for use in |
| 3 |
|
|
* OpenMD |
| 4 |
|
|
* |
| 5 |
|
|
* Matpack Library Release 1.9.0 |
| 6 |
|
|
* Copyright (C) 1991-2003 by Berndt M. Gammel. All rights reserved. |
| 7 |
|
|
* |
| 8 |
|
|
* Permission to use, copy, and distribute Matpack in its entirety |
| 9 |
|
|
* and its documentation for non-commercial purpose and without fee |
| 10 |
|
|
* is hereby granted, provided that this license information and |
| 11 |
|
|
* copyright notice appear unmodified in all copies. This software |
| 12 |
|
|
* is provided 'as is' without express or implied warranty. In no |
| 13 |
|
|
* event will the author be held liable for any damages arising from |
| 14 |
|
|
* the use of this software. |
| 15 |
|
|
* |
| 16 |
|
|
* Note that distributing Matpack 'bundled' in with any product is |
| 17 |
|
|
* considered to be a 'commercial purpose'. |
| 18 |
|
|
* |
| 19 |
|
|
* The software may be modified for your own purposes, but modified |
| 20 |
|
|
* versions may not be distributed without prior consent of the |
| 21 |
|
|
* author. |
| 22 |
|
|
* |
| 23 |
|
|
* Read the COPYRIGHT and README files in this distribution about |
| 24 |
|
|
* registration and installation of Matpack. |
| 25 |
|
|
*/ |
| 26 |
|
|
|
| 27 |
|
|
#include "Wigner3jm.hpp" |
| 28 |
|
|
#include <cmath> |
| 29 |
|
|
#include <cfloat> |
| 30 |
|
|
#include <cstdio> |
| 31 |
|
|
#include "utils/simError.h" |
| 32 |
|
|
|
| 33 |
|
|
namespace MATPACK { |
| 34 |
|
|
|
| 35 |
|
|
//-----------------------------------------------------------------------------// |
| 36 |
|
|
// |
| 37 |
|
|
// void ThreeJSymbolM (RealType l1, RealType l2, RealType l3, RealType m1, |
| 38 |
|
|
// RealType &m2min, RealType &m2max, RealType *thrcof, int ndim, |
| 39 |
|
|
// int &errflag) |
| 40 |
|
|
// |
| 41 |
|
|
// Evaluate the Wigner 3j symbol |
| 42 |
|
|
// |
| 43 |
|
|
// g(m2) = ( l1 l2 l3 ) |
| 44 |
|
|
// ( m1 m2 -m1-m2 ) |
| 45 |
|
|
// |
| 46 |
|
|
// for all allowed values of m2, the other parameters being held fixed. |
| 47 |
|
|
// |
| 48 |
|
|
// Input Arguments: |
| 49 |
|
|
// ---------------- |
| 50 |
|
|
// |
| 51 |
|
|
// RealType l1 |
| 52 |
|
|
// RealType l2 |
| 53 |
|
|
// RealType l3 |
| 54 |
|
|
// RealType m1 Parameters in 3j symbol. |
| 55 |
|
|
// |
| 56 |
|
|
// int ndim Declared length of thrcof in calling program. |
| 57 |
|
|
// |
| 58 |
|
|
// Output Arguments: |
| 59 |
|
|
// ----------------- |
| 60 |
|
|
// |
| 61 |
|
|
// RealType &m2min Smallest allowable m2 in 3j symbol. |
| 62 |
|
|
// RealType &m2max Largest allowable m2 in 3j symbol. |
| 63 |
|
|
// RealType *thrcof Set of 3j coefficients generated by evaluating the |
| 64 |
|
|
// 3j symbol for all allowed values of m2. thrcof(i) |
| 65 |
|
|
// will contain g(m2min+i), i=0,2,...,m2max-m2min. |
| 66 |
|
|
// |
| 67 |
|
|
// int &errflag Error flag. |
| 68 |
|
|
// errflag=0 No errors. |
| 69 |
|
|
// errflag=1 Either l1 < abs(m1) or l1+abs(m1) non-integer. |
| 70 |
|
|
// errflag=2 abs(l1-l2)<= l3 <= l1+l2 not satisfied. |
| 71 |
|
|
// errflag=3 l1+l2+l3 not an integer. |
| 72 |
|
|
// errflag=4 m2max-m2min not an integer. |
| 73 |
|
|
// errflag=5 m2max less than m2min. |
| 74 |
|
|
// errflag=6 ndim less than m2max-m2min+1. |
| 75 |
|
|
// Description: |
| 76 |
|
|
// ------------ |
| 77 |
|
|
// |
| 78 |
|
|
// Although conventionally the parameters of the vector addition |
| 79 |
|
|
// coefficients satisfy certain restrictions, such as being integers |
| 80 |
|
|
// or integers plus 1/2, the restrictions imposed on input to this |
| 81 |
|
|
// subroutine are somewhat weaker. See, for example, Section 27.9 of |
| 82 |
|
|
// Abramowitz and Stegun or Appendix C of Volume II of A. Messiah. |
| 83 |
|
|
// |
| 84 |
|
|
// The restrictions imposed by this subroutine are |
| 85 |
|
|
// |
| 86 |
|
|
// 1. l1 >= abs(m1) and l1+abs(m1) must be an integer |
| 87 |
|
|
// 2. abs(l1-l2) <= l3 <= l1+l2 |
| 88 |
|
|
// 3. l1+l2+l3 must be an integer |
| 89 |
|
|
// 4. m2max-m2min must be an integer, where |
| 90 |
|
|
// m2max=min(l2,l3-m1) and m2min=max(-l2,-l3-m1) |
| 91 |
|
|
// |
| 92 |
|
|
// If the conventional restrictions are satisfied, then these |
| 93 |
|
|
// restrictions are also met. |
| 94 |
|
|
// |
| 95 |
|
|
// The user should be cautious in using input parameters that do |
| 96 |
|
|
// not satisfy the conventional restrictions. For example, the |
| 97 |
|
|
// the subroutine produces values of |
| 98 |
|
|
// g(m2) = (0.75 1.50 1.75 ) |
| 99 |
|
|
// (0.25 m2 -0.25-m2) |
| 100 |
|
|
// for m2=-1.5,-0.5,0.5,1.5 but none of the symmetry properties of the |
| 101 |
|
|
// 3j symbol, set forth on page 1056 of Messiah, is satisfied. |
| 102 |
|
|
// |
| 103 |
|
|
// The subroutine generates g(m2min), g(m2min+1), ..., g(m2max) |
| 104 |
|
|
// where m2min and m2max are defined above. The sequence g(m2) is |
| 105 |
|
|
// generated by a three-term recurrence algorithm with scaling to |
| 106 |
|
|
// control overflow. Both backward and forward recurrence are used to |
| 107 |
|
|
// maintain numerical stability. The two recurrence sequences are |
| 108 |
|
|
// matched at an interior point and are normalized from the unitary |
| 109 |
|
|
// property of 3j coefficients and Wigner's phase convention. |
| 110 |
|
|
// |
| 111 |
|
|
// The algorithm is suited to applications in which large quantum |
| 112 |
|
|
// numbers arise, such as in molecular dynamics. |
| 113 |
|
|
// |
| 114 |
|
|
// References: |
| 115 |
|
|
// ----------- |
| 116 |
|
|
// 1. Abramowitz, M., and Stegun, I. A., Eds., Handbook |
| 117 |
|
|
// of Mathematical Functions with Formulas, Graphs |
| 118 |
|
|
// and Mathematical Tables, NBS Applied Mathematics |
| 119 |
|
|
// Series 55, June 1964 and subsequent printings. |
| 120 |
|
|
// 2. Messiah, Albert., Quantum Mechanics, Volume II, |
| 121 |
|
|
// North-Holland Publishing Company, 1963. |
| 122 |
|
|
// 3. Schulten, Klaus and Gordon, Roy G., Exact recursive |
| 123 |
|
|
// evaluation of 3j and 6j coefficients for quantum- |
| 124 |
|
|
// mechanical coupling of angular momenta, J Math |
| 125 |
|
|
// Phys, v 16, no. 10, October 1975, pp. 1961-1970. |
| 126 |
|
|
// 4. Schulten, Klaus and Gordon, Roy G., Semiclassical |
| 127 |
|
|
// approximations to 3j and 6j coefficients for |
| 128 |
|
|
// quantum-mechanical coupling of angular momenta, |
| 129 |
|
|
// J Math Phys, v 16, no. 10, October 1975, pp. 1971-1988. |
| 130 |
|
|
// 5. Schulten, Klaus and Gordon, Roy G., Recursive |
| 131 |
|
|
// evaluation of 3j and 6j coefficients, Computer |
| 132 |
|
|
// Phys Comm, v 11, 1976, pp. 269-278. |
| 133 |
|
|
// 6. SLATEC library, category C19, |
| 134 |
|
|
// double precision algorithm DRC3JM.F |
| 135 |
|
|
// Keywords: 3j coefficients, 3j symbols, Clebsch-Gordan coefficients, |
| 136 |
|
|
// Racah coefficients, vector addition coefficients, |
| 137 |
|
|
// Wigner coefficients |
| 138 |
|
|
// Author: Gordon, R. G., Harvard University |
| 139 |
|
|
// Schulten, K., Max Planck Institute |
| 140 |
|
|
// Revision history (YYMMDD) |
| 141 |
|
|
// 750101 DATE WRITTEN |
| 142 |
|
|
// 880515 SLATEC prologue added by G. C. Nielson, NBS; parameters |
| 143 |
|
|
// HUGE and TINY revised to depend on D1MACH. |
| 144 |
|
|
// 891229 Prologue description rewritten; other prologue sections |
| 145 |
|
|
// revised; MMATCH (location of match point for recurrences) |
| 146 |
|
|
// removed from argument list; argument IER changed to serve |
| 147 |
|
|
// only as an error flag (previously, in cases without error, |
| 148 |
|
|
// it returned the number of scalings); number of error codes |
| 149 |
|
|
// increased to provide more precise error information; |
| 150 |
|
|
// program comments revised; SLATEC error handler calls |
| 151 |
|
|
// introduced to enable printing of error messages to meet |
| 152 |
|
|
// SLATEC standards. These changes were done by D. W. Lozier, |
| 153 |
|
|
// M. A. McClain and J. M. Smith of the National Institute |
| 154 |
|
|
// of Standards and Technology, formerly NBS. |
| 155 |
|
|
// 910415 Mixed type expressions eliminated; variable C1 initialized; |
| 156 |
|
|
// description of THRCOF expanded. These changes were done by |
| 157 |
|
|
// D. W. Lozier. |
| 158 |
|
|
// 7. Rewritting of the SLATEX algorithm in C++ and adaption to the |
| 159 |
|
|
// Matpack C++ Numerics and Graphics Library by Berndt M. Gammel |
| 160 |
|
|
// in June 1997. |
| 161 |
|
|
// |
| 162 |
|
|
//-----------------------------------------------------------------------------// |
| 163 |
|
|
|
| 164 |
|
|
|
| 165 |
|
|
void Wigner3jm(RealType l1, RealType l2, RealType l3, RealType m1, |
| 166 |
|
|
RealType &m2min, RealType &m2max, RealType *thrcof, int ndim, |
| 167 |
|
|
int &errflag) { |
| 168 |
|
|
|
| 169 |
|
|
// In single precision, the largest floating point number is not |
| 170 |
|
|
// the same as in double precision: |
| 171 |
|
|
#ifdef SINGLE_PRECISION |
| 172 |
|
|
RealType MaxFloat = FLT_MAX; |
| 173 |
|
|
#else |
| 174 |
|
|
RealType MaxFloat = DBL_MAX; |
| 175 |
|
|
#endif |
| 176 |
|
|
|
| 177 |
|
|
const RealType zero = 0.0, eps = 0.01, one = 1.0, two = 2.0; |
| 178 |
|
|
|
| 179 |
|
|
int nfin, nlim, i, n, index, lstep, nfinp1, nfinp2, nfinp3, nstep2; |
| 180 |
|
|
RealType oldfac, dv, newfac, sumbac = 0.0, thresh, a1s, sumfor, sumuni, |
| 181 |
|
|
sum1, sum2, x, y, m2, m3, x1, x2, x3, y1, y2, y3, cnorm, |
| 182 |
|
|
ratio, a1, c1, c2, c1old = 0.0, sign1, sign2; |
| 183 |
|
|
|
| 184 |
|
|
// Parameter adjustments |
| 185 |
|
|
--thrcof; |
| 186 |
|
|
|
| 187 |
|
|
errflag = 0; |
| 188 |
|
|
|
| 189 |
|
|
// "hugeRealType" is the square root of one twentieth of the |
| 190 |
|
|
// largest floating point number, approximately. |
| 191 |
|
|
|
| 192 |
|
|
RealType hugeRealType = sqrt(MaxFloat / 20.0), |
| 193 |
|
|
srhuge = sqrt(hugeRealType), |
| 194 |
|
|
tiny = one / hugeRealType, |
| 195 |
|
|
srtiny = one / srhuge; |
| 196 |
|
|
|
| 197 |
|
|
// lmatch = zero |
| 198 |
|
|
|
| 199 |
|
|
// Check error conditions 1, 2, and 3. |
| 200 |
|
|
if (l1 - fabs(m1) + eps < zero |
| 201 |
|
|
|| fmod(l1 + fabs(m1) + eps, one) >= eps + eps) { |
| 202 |
|
|
errflag = 1; |
| 203 |
|
|
|
| 204 |
|
|
sprintf( painCave.errMsg, "%s: %s", "ThreeJSymbolM", |
| 205 |
|
|
"l1-abs(m1) less than zero or l1+abs(m1) not integer."); |
| 206 |
|
|
painCave.isFatal = 1; |
| 207 |
|
|
simError(); |
| 208 |
|
|
|
| 209 |
|
|
return; |
| 210 |
|
|
} else if (l1+l2-l3 < -eps || l1-l2+l3 < -eps || -(l1) + l2+l3 < -eps) { |
| 211 |
|
|
errflag = 2; |
| 212 |
|
|
|
| 213 |
|
|
sprintf( painCave.errMsg, "%s: %s", "ThreeJSymbolM", |
| 214 |
|
|
"l1, l2, l3 do not satisfy triangular condition."); |
| 215 |
|
|
painCave.isFatal = 1; |
| 216 |
|
|
simError(); |
| 217 |
|
|
|
| 218 |
|
|
return; |
| 219 |
|
|
} else if (fmod(l1 + l2 + l3 + eps, one) >= eps + eps) { |
| 220 |
|
|
errflag = 3; |
| 221 |
|
|
|
| 222 |
|
|
sprintf( painCave.errMsg, "%s: %s", "ThreeJSymbolM", |
| 223 |
|
|
"l1+l2+l3 not integer."); |
| 224 |
|
|
painCave.isFatal = 1; |
| 225 |
|
|
simError(); |
| 226 |
|
|
|
| 227 |
|
|
return; |
| 228 |
|
|
} |
| 229 |
|
|
|
| 230 |
|
|
// limits for m2 |
| 231 |
|
|
m2min = MpMax(-l2,-l3-m1); |
| 232 |
|
|
m2max = MpMin(l2,l3-m1); |
| 233 |
|
|
|
| 234 |
|
|
// Check error condition 4. |
| 235 |
|
|
if (fmod(m2max - m2min + eps, one) >= eps + eps) { |
| 236 |
|
|
errflag = 4; |
| 237 |
|
|
|
| 238 |
|
|
sprintf( painCave.errMsg, "%s: %s", "ThreeJSymbolM", |
| 239 |
|
|
"m2max-m2min not integer."); |
| 240 |
|
|
painCave.isFatal = 1; |
| 241 |
|
|
simError(); |
| 242 |
|
|
|
| 243 |
|
|
return; |
| 244 |
|
|
} |
| 245 |
|
|
if (m2min < m2max - eps) goto L20; |
| 246 |
|
|
if (m2min < m2max + eps) goto L10; |
| 247 |
|
|
|
| 248 |
|
|
// Check error condition 5. |
| 249 |
|
|
errflag = 5; |
| 250 |
|
|
|
| 251 |
|
|
sprintf( painCave.errMsg, "%s: %s", "ThreeJSymbolM", |
| 252 |
|
|
"m2min greater than m2max."); |
| 253 |
|
|
painCave.isFatal = 1; |
| 254 |
|
|
simError(); |
| 255 |
|
|
|
| 256 |
|
|
return; |
| 257 |
|
|
|
| 258 |
|
|
// This is reached in case that m2 and m3 can take only one value. |
| 259 |
|
|
L10: |
| 260 |
|
|
// mscale = 0 |
| 261 |
|
|
thrcof[1] = (odd(int(fabs(l2-l3-m1)+eps)) ? -one : one) / sqrt(l1+l2+l3+one); |
| 262 |
|
|
return; |
| 263 |
|
|
|
| 264 |
|
|
// This is reached in case that M1 and M2 take more than one |
| 265 |
|
|
// value. |
| 266 |
|
|
L20: |
| 267 |
|
|
// mscale = 0 |
| 268 |
|
|
nfin = int(m2max - m2min + one + eps); |
| 269 |
|
|
if (ndim - nfin >= 0) goto L23; |
| 270 |
|
|
|
| 271 |
|
|
// Check error condition 6. |
| 272 |
|
|
|
| 273 |
|
|
errflag = 6; |
| 274 |
|
|
sprintf( painCave.errMsg, "%s: %s", "ThreeJSymbolM", |
| 275 |
|
|
"Dimension of result array for 3j coefficients too small."); |
| 276 |
|
|
painCave.isFatal = 1; |
| 277 |
|
|
simError(); |
| 278 |
|
|
|
| 279 |
|
|
return; |
| 280 |
|
|
|
| 281 |
|
|
// Start of forward recursion from m2 = m2min |
| 282 |
|
|
|
| 283 |
|
|
L23: |
| 284 |
|
|
m2 = m2min; |
| 285 |
|
|
thrcof[1] = srtiny; |
| 286 |
|
|
newfac = 0.0; |
| 287 |
|
|
c1 = 0.0; |
| 288 |
|
|
sum1 = tiny; |
| 289 |
|
|
|
| 290 |
|
|
lstep = 1; |
| 291 |
|
|
L30: |
| 292 |
|
|
++lstep; |
| 293 |
|
|
m2 += one; |
| 294 |
|
|
m3 = -m1 - m2; |
| 295 |
|
|
|
| 296 |
|
|
oldfac = newfac; |
| 297 |
|
|
a1 = (l2 - m2 + one) * (l2 + m2) * (l3 + m3 + one) * (l3 - m3); |
| 298 |
|
|
newfac = sqrt(a1); |
| 299 |
|
|
|
| 300 |
|
|
dv = (l1+l2+l3+one) * (l2+l3-l1) - (l2-m2+one) * (l3+m3+one) |
| 301 |
|
|
- (l2+m2-one) * (l3-m3-one); |
| 302 |
|
|
|
| 303 |
|
|
if (lstep - 2 > 0) c1old = fabs(c1); |
| 304 |
|
|
|
| 305 |
|
|
// L32: |
| 306 |
|
|
c1 = -dv / newfac; |
| 307 |
|
|
|
| 308 |
|
|
if (lstep > 2) goto L60; |
| 309 |
|
|
|
| 310 |
|
|
// If m2 = m2min + 1, the third term in the recursion equation |
| 311 |
|
|
// vanishes, hence |
| 312 |
|
|
|
| 313 |
|
|
x = srtiny * c1; |
| 314 |
|
|
thrcof[2] = x; |
| 315 |
|
|
sum1 += tiny * c1 * c1; |
| 316 |
|
|
if (lstep == nfin) goto L220; |
| 317 |
|
|
goto L30; |
| 318 |
|
|
|
| 319 |
|
|
L60: |
| 320 |
|
|
c2 = -oldfac / newfac; |
| 321 |
|
|
|
| 322 |
|
|
// Recursion to the next 3j coefficient |
| 323 |
|
|
x = c1 * thrcof[lstep-1] + c2 * thrcof[lstep-2]; |
| 324 |
|
|
thrcof[lstep] = x; |
| 325 |
|
|
sumfor = sum1; |
| 326 |
|
|
sum1 += x * x; |
| 327 |
|
|
if (lstep == nfin) goto L100; |
| 328 |
|
|
|
| 329 |
|
|
// See if last unnormalized 3j coefficient exceeds srhuge |
| 330 |
|
|
|
| 331 |
|
|
if (fabs(x) < srhuge) goto L80; |
| 332 |
|
|
|
| 333 |
|
|
// This is reached if last 3j coefficient larger than srhuge, so |
| 334 |
|
|
// that the recursion series thrcof(1), ... , thrcof(lstep) has to |
| 335 |
|
|
// be rescaled to prevent overflow |
| 336 |
|
|
|
| 337 |
|
|
// mscale = mscale + 1 |
| 338 |
|
|
for (i = 1; i <= lstep; ++i) { |
| 339 |
|
|
if (fabs(thrcof[i]) < srtiny) thrcof[i] = zero; |
| 340 |
|
|
thrcof[i] /= srhuge; |
| 341 |
|
|
} |
| 342 |
|
|
sum1 /= hugeRealType; |
| 343 |
|
|
sumfor /= hugeRealType; |
| 344 |
|
|
x /= srhuge; |
| 345 |
|
|
|
| 346 |
|
|
// As long as abs(c1) is decreasing, the recursion proceeds |
| 347 |
|
|
// towards increasing 3j values and, hence, is numerically stable. |
| 348 |
|
|
// Once an increase of abs(c1) is detected, the recursion |
| 349 |
|
|
// direction is reversed. |
| 350 |
|
|
|
| 351 |
|
|
L80: |
| 352 |
|
|
if (c1old - fabs(c1) > 0.0) goto L30; |
| 353 |
|
|
|
| 354 |
|
|
// Keep three 3j coefficients around mmatch for comparison later |
| 355 |
|
|
// with backward recursion values. |
| 356 |
|
|
|
| 357 |
|
|
L100: |
| 358 |
|
|
// mmatch = m2 - 1 |
| 359 |
|
|
nstep2 = nfin - lstep + 3; |
| 360 |
|
|
x1 = x; |
| 361 |
|
|
x2 = thrcof[lstep-1]; |
| 362 |
|
|
x3 = thrcof[lstep-2]; |
| 363 |
|
|
|
| 364 |
|
|
// Starting backward recursion from m2max taking nstep2 steps, so |
| 365 |
|
|
// that forwards and backwards recursion overlap at the three |
| 366 |
|
|
// points m2 = mmatch+1, mmatch, mmatch-1. |
| 367 |
|
|
|
| 368 |
|
|
nfinp1 = nfin + 1; |
| 369 |
|
|
nfinp2 = nfin + 2; |
| 370 |
|
|
nfinp3 = nfin + 3; |
| 371 |
|
|
thrcof[nfin] = srtiny; |
| 372 |
|
|
sum2 = tiny; |
| 373 |
|
|
|
| 374 |
|
|
m2 = m2max + two; |
| 375 |
|
|
lstep = 1; |
| 376 |
|
|
L110: |
| 377 |
|
|
++lstep; |
| 378 |
|
|
m2 -= one; |
| 379 |
|
|
m3 = -m1 - m2; |
| 380 |
|
|
oldfac = newfac; |
| 381 |
|
|
a1s = (l2-m2+two) * (l2+m2-one) * (l3+m3+two) * (l3-m3-one); |
| 382 |
|
|
newfac = sqrt(a1s); |
| 383 |
|
|
dv = (l1+l2+l3+one) * (l2+l3-l1) - (l2-m2+one) * (l3+m3+one) |
| 384 |
|
|
- (l2+m2-one) * (l3-m3-one); |
| 385 |
|
|
c1 = -dv / newfac; |
| 386 |
|
|
if (lstep > 2) goto L120; |
| 387 |
|
|
|
| 388 |
|
|
// if m2 = m2max + 1 the third term in the recursion equation |
| 389 |
|
|
// vanishes |
| 390 |
|
|
|
| 391 |
|
|
y = srtiny * c1; |
| 392 |
|
|
thrcof[nfin - 1] = y; |
| 393 |
|
|
if (lstep == nstep2) goto L200; |
| 394 |
|
|
sumbac = sum2; |
| 395 |
|
|
sum2 += y * y; |
| 396 |
|
|
goto L110; |
| 397 |
|
|
|
| 398 |
|
|
L120: |
| 399 |
|
|
c2 = -oldfac / newfac; |
| 400 |
|
|
|
| 401 |
|
|
// Recursion to the next 3j coefficient |
| 402 |
|
|
|
| 403 |
|
|
y = c1 * thrcof[nfinp2 - lstep] + c2 * thrcof[nfinp3 - lstep]; |
| 404 |
|
|
|
| 405 |
|
|
if (lstep == nstep2) goto L200; |
| 406 |
|
|
|
| 407 |
|
|
thrcof[nfinp1 - lstep] = y; |
| 408 |
|
|
sumbac = sum2; |
| 409 |
|
|
sum2 += y * y; |
| 410 |
|
|
|
| 411 |
|
|
// See if last 3j coefficient exceeds SRHUGE |
| 412 |
|
|
|
| 413 |
|
|
if (fabs(y) < srhuge) goto L110; |
| 414 |
|
|
|
| 415 |
|
|
// This is reached if last 3j coefficient larger than srhuge, so |
| 416 |
|
|
// that the recursion series |
| 417 |
|
|
// thrcof(nfin), ... , thrcof(nfin-lstep+1) |
| 418 |
|
|
// has to be rescaled to prevent overflow. |
| 419 |
|
|
|
| 420 |
|
|
// mscale = mscale + 1 |
| 421 |
|
|
for (i = 1; i <= lstep; ++i) { |
| 422 |
|
|
index = nfin - i + 1; |
| 423 |
|
|
if (fabs(thrcof[index]) < srtiny) thrcof[index] = zero; |
| 424 |
|
|
thrcof[index] /= srhuge; |
| 425 |
|
|
} |
| 426 |
|
|
sum2 /= hugeRealType; |
| 427 |
|
|
sumbac /= hugeRealType; |
| 428 |
|
|
|
| 429 |
|
|
goto L110; |
| 430 |
|
|
|
| 431 |
|
|
// The forward recursion 3j coefficients x1, x2, x3 are to be |
| 432 |
|
|
// matched with the corresponding backward recursion values y1, |
| 433 |
|
|
// y2, y3. |
| 434 |
|
|
|
| 435 |
|
|
L200: |
| 436 |
|
|
y3 = y; |
| 437 |
|
|
y2 = thrcof[nfinp2-lstep]; |
| 438 |
|
|
y1 = thrcof[nfinp3-lstep]; |
| 439 |
|
|
|
| 440 |
|
|
// Determine now ratio such that yi = ratio * xi (i=1,2,3) holds |
| 441 |
|
|
// with minimal error. |
| 442 |
|
|
|
| 443 |
|
|
ratio = (x1*y1 + x2*y2 + x3*y3) / (x1*x1 + x2*x2 + x3*x3); |
| 444 |
|
|
nlim = nfin - nstep2 + 1; |
| 445 |
|
|
|
| 446 |
|
|
if (fabs(ratio) < one) goto L211; |
| 447 |
|
|
for (n = 1; n <= nlim; ++n) |
| 448 |
|
|
thrcof[n] = ratio * thrcof[n]; |
| 449 |
|
|
sumuni = ratio * ratio * sumfor + sumbac; |
| 450 |
|
|
goto L230; |
| 451 |
|
|
|
| 452 |
|
|
L211: |
| 453 |
|
|
++nlim; |
| 454 |
|
|
ratio = one / ratio; |
| 455 |
|
|
for (n = nlim; n <= nfin; ++n) |
| 456 |
|
|
thrcof[n] = ratio * thrcof[n]; |
| 457 |
|
|
sumuni = sumfor + ratio * ratio * sumbac; |
| 458 |
|
|
goto L230; |
| 459 |
|
|
|
| 460 |
|
|
L220: |
| 461 |
|
|
sumuni = sum1; |
| 462 |
|
|
|
| 463 |
|
|
// Normalize 3j coefficients |
| 464 |
|
|
|
| 465 |
|
|
L230: |
| 466 |
|
|
cnorm = one / sqrt((l1+l1+one) * sumuni); |
| 467 |
|
|
|
| 468 |
|
|
// Sign convention for last 3j coefficient determines overall |
| 469 |
|
|
// phase |
| 470 |
|
|
|
| 471 |
|
|
sign1 = sign(thrcof[nfin]); |
| 472 |
|
|
sign2 = odd(int(fabs(l2-l3-m1)+eps)) ? -one : one; |
| 473 |
|
|
if (sign1 * sign2 <= 0.0) goto L235; |
| 474 |
|
|
else goto L236; |
| 475 |
|
|
|
| 476 |
|
|
L235: |
| 477 |
|
|
cnorm = -cnorm; |
| 478 |
|
|
|
| 479 |
|
|
L236: |
| 480 |
|
|
if (fabs(cnorm) < one) goto L250; |
| 481 |
|
|
|
| 482 |
|
|
for (n = 1; n <= nfin; ++n) |
| 483 |
|
|
thrcof[n] = cnorm * thrcof[n]; |
| 484 |
|
|
return; |
| 485 |
|
|
|
| 486 |
|
|
L250: |
| 487 |
|
|
thresh = tiny / fabs(cnorm); |
| 488 |
|
|
for (n = 1; n <= nfin; ++n) { |
| 489 |
|
|
if (fabs(thrcof[n]) < thresh) thrcof[n] = zero; |
| 490 |
|
|
thrcof[n] = cnorm * thrcof[n]; |
| 491 |
|
|
} |
| 492 |
|
|
} |
| 493 |
|
|
} |
