src/math/Wigner3jm.F90

subroutine WIGNER3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, &
     IER)
 ! use status
!
!! DRC3JM evaluates the 3j symbol g(M2) for all allowed values of M2.
!
!***PURPOSE  Evaluate the 3j symbol g(M2) = (L1 L2   L3  )
!                                           (M1 M2 -M1-M2)
!            for all allowed values of M2, the other parameters
!            being held fixed.
!
!***LIBRARY   SLATEC
!***CATEGORY  C19
!***TYPE      DOUBLE PRECISION (RC3JM-S, DRC3JM-D)
!***KEYWORDS  3J COEFFICIENTS, 3J SYMBOLS, CLEBSCH-GORDAN COEFFICIENTS,
!             RACAH COEFFICIENTS, VECTOR ADDITION COEFFICIENTS,
!             WIGNER COEFFICIENTS
!***AUTHOR  Gordon, R. G., Harvard University
!           Schulten, K., Max Planck Institute
!***DESCRIPTION
!
! *Usage:
!
!        DOUBLE PRECISION L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM)
!        INTEGER NDIM, IER
!
!        call DRC3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, IER)
!
! *Arguments:
!
!     L1 :IN      Parameter in 3j symbol.
!
!     L2 :IN      Parameter in 3j symbol.
!
!     L3 :IN      Parameter in 3j symbol.
!
!     M1 :IN      Parameter in 3j symbol.
!
!     M2MIN :OUT  Smallest allowable M2 in 3j symbol.
!
!     M2MAX :OUT  Largest allowable M2 in 3j symbol.
!
!     THRCOF :OUT Set of 3j coefficients generated by evaluating the
!                 3j symbol for all allowed values of M2.  THRCOF(I)
!                 will contain g(M2MIN+I-1), I=1,2,...,M2MAX-M2MIN+1.
!
!     NDIM :IN    Declared length of THRCOF in calling program.
!
!     IER :OUT    Error flag.
!                 IER=0 No errors.
!                 IER=1 Either L1 < ABS(M1) or L1+ABS(M1) non-integer.
!                 IER=2 ABS(L1-L2) <= L3 <= L1+L2 not satisfied.
!                 IER=3 L1+L2+L3 not an integer.
!                 IER=4 M2MAX-M2MIN not an integer.
!                 IER=5 M2MAX less than M2MIN.
!                 IER=6 NDIM less than M2MAX-M2MIN+1.
!
! *Description:
!
!     Although conventionally the parameters of the vector addition
!  coefficients satisfy certain restrictions, such as being integers
!  or integers plus 1/2, the restrictions imposed on input to this
!  subroutine are somewhat weaker. See, for example, Section 27.9 of
!  Abramowitz and Stegun or Appendix C of Volume II of A. Messiah.
!  The restrictions imposed by this subroutine are
!       1. L1 >= ABS(M1) and L1+ABS(M1) must be an integer;
!       2. ABS(L1-L2) <= L3 <= L1+L2;
!       3. L1+L2+L3 must be an integer;
!       4. M2MAX-M2MIN must be an integer, where
!          M2MAX=MIN(L2,L3-M1) and M2MIN=MAX(-L2,-L3-M1).
!  If the conventional restrictions are satisfied, then these
!  restrictions are met.
!
!     The user should be cautious in using input parameters that do
!  not satisfy the conventional restrictions. For example, the
!  the subroutine produces values of
!       g(M2) = (0.751.50   1.75  )
!               (0.25  M2  -0.25-M2)
!  for M2=-1.5,-0.5,0.5,1.5 but none of the symmetry properties of the
!  3j symbol, set forth on page 1056 of Messiah, is satisfied.
!
!     The subroutine generates g(M2MIN), g(M2MIN+1), ..., g(M2MAX)
!  where M2MIN and M2MAX are defined above. The sequence g(M2) is
!  generated by a three-term recurrence algorithm with scaling to
!  control overflow. Both backward and forward recurrence are used to
!  maintain numerical stability. The two recurrence sequences are
!  matched at an interior point and are normalized from the unitary
!  property of 3j coefficients and Wigner's phase convention.
!
!    The algorithm is suited to applications in which large quantum
!  numbers arise, such as in molecular dynamics.
!
!***REFERENCES  1. Abramowitz, M., and Stegun, I. A., Eds., Handbook
!                  of Mathematical Functions with Formulas, Graphs
!                  and Mathematical Tables, NBS Applied Mathematics
!                  Series 55, June 1964 and subsequent printings.
!               2. Messiah, Albert., Quantum Mechanics, Volume II,
!                  North-Holland Publishing Company, 1963.
!               3. Schulten, Klaus and Gordon, Roy G., Exact recursive
!                  evaluation of 3j and 6j coefficients for quantum-
!                  mechanical coupling of angular momenta, J Math
!                  Phys, v 16, no. 10, October 1975, pp. 1961-1970.
!               4. Schulten, Klaus and Gordon, Roy G., Semiclassical
!                  approximations to 3j and 6j coefficients for
!                  quantum-mechanical coupling of angular momenta,
!                  J Math Phys, v 16, no. 10, October 1975,
!                  pp. 1971-1988.
!               5. Schulten, Klaus and Gordon, Roy G., Recursive
!                  evaluation of 3j and 6j coefficients, Computer
!                  Phys Comm, v 11, 1976, pp. 269-278.
!***ROUTINES CALLED  D1MACH, XERMSG
!***REVISION HISTORY  (YYMMDD)
!   750101  DATE WRITTEN
!   880515  SLATEC prologue added by G. C. Nielson, NBS; parameters
!           HUGE and TINY revised to depend on D1MACH.
!   891229  Prologue description rewritten; other prologue sections
!           revised; MMATCH (location of match point for recurrences)
!           removed from argument list; argument IER changed to serve
!           only as an error flag (previously, in cases without error,
!           it returned the number of scalings); number of error codes
!           increased to provide more precise error information;
!           program comments revised; SLATEC error handler calls
!           introduced to enable printing of error messages to meet
!           SLATEC standards. These changes were done by D. W. Lozier,
!           M. A. McClain and J. M. Smith of the National Institute
!           of Standards and Technology, formerly NBS.
!   910415  Mixed type expressions eliminated; variable C1 initialized;
!           description of THRCOF expanded. These changes were done by
!           D. W. Lozier.
!***END PROLOGUE  DRC3JM
!
  INTEGER NDIM, IER
  DOUBLE PRECISION L1, L2, L3, M1, M2MIN, M2MAX, THRCOF(NDIM)
!
  INTEGER I, INDEX, LSTEP, N, NFIN, NFINP1, NFINP2, NFINP3, NLIM, &
          NSTEP2
  DOUBLE PRECISION A1, A1S, C1, C1OLD, C2, CNORM, D1MACH, DV, EPS, &
                   HUGE, M2, M3, NEWFAC, OLDFAC, ONE, RATIO, SIGN1, &
                   SIGN2, SRHUGE, SRTINY, SUM1, SUM2, SUMBAC, &
                   SUMFOR, SUMUNI, THRESH, TINY, TWO, X, X1, X2, X3, &
                   Y, Y1, Y2, Y3, ZERO
!
  DATA  ZERO,EPS,ONE,TWO /0.0D0,0.01D0,1.0D0,2.0D0/
!
!***FIRST EXECUTABLE STATEMENT  DRC3JM
  IER=0
!  HUGE is the square root of one twentieth of the largest floating
!  point number, approximately.
  HUGE = SQRT(D1MACH(2)/20.0D0)
  SRHUGE = SQRT(HUGE)
  TINY = 1.0D0/HUGE
  SRTINY = 1.0D0/SRHUGE
!
!     MMATCH = ZERO
!
!
!  Check error conditions 1, 2, and 3.
  if ( (L1-ABS(M1)+EPS < ZERO).OR. &
     (MOD(L1+ABS(M1)+EPS,ONE) >= EPS+EPS))THEN
     IER=1
!     call handleError('Wigner3jm','L1-ABS(M1) less than zero or '// &
!        'L1+ABS(M1) not integer.')
     return
  ELSEIF((L1+L2-L3 < -EPS).OR.(L1-L2+L3 < -EPS).OR. &
     (-L1+L2+L3 < -EPS))THEN
     IER=2
     write(*,*) eps
     write(*,*) L1,L2,L3
     write(*,*) "L1,L2,L3 do not satisfy triangular condition"
!     call handleError('Wigner3jm','L1, L2, L3 do not satisfy '// &
!        'triangular condition.')
     return
  ELSEIF(MOD(L1+L2+L3+EPS,ONE) >= EPS+EPS)THEN
     IER=3
!     call handleError('Wigner3jm','L1+L2+L3 not integer.')
     return
  end if
!
!
!  Limits for M2
  M2MIN = MAX(-L2,-L3-M1)
  M2MAX = MIN(L2,L3-M1)
!
!  Check error condition 4.
  if ( MOD(M2MAX-M2MIN+EPS,ONE) >= EPS+EPS)THEN
     IER=4
!     call handleError('Wigner3jm','M2MAX-M2MIN not integer.')
     return
  end if
  if ( M2MIN < M2MAX-EPS)   go to 20
  if ( M2MIN < M2MAX+EPS)   go to 10
!
!  Check error condition 5.
  IER=5
!  call handleError('Wigner3jm','M2MIN greater than M2MAX.')
  return
!
!
!  This is reached in case that M2 and M3 can take only one value.
   10 CONTINUE
!     MSCALE = 0
  THRCOF(1) = (-ONE) ** INT(ABS(L2-L3-M1)+EPS) / &
   SQRT(L1+L2+L3+ONE)
  return
!
!  This is reached in case that M1 and M2 take more than one value.
   20 CONTINUE
!     MSCALE = 0
  NFIN = INT(M2MAX-M2MIN+ONE+EPS)
  if ( NDIM-NFIN)   21, 23, 23
!
!  Check error condition 6.
   21 IER = 6
!  call handleError('Wigner3jm','Dimension of result array for '// &
!              '3j coefficients too small.')
  return
!
!
!
!  Start of forward recursion from M2 = M2MIN
!
   23 M2 = M2MIN
  THRCOF(1) = SRTINY
  NEWFAC = 0.0D0
  C1 = 0.0D0
  SUM1 = TINY
!
!
  LSTEP = 1
   30 LSTEP = LSTEP + 1
  M2 = M2 + ONE
  M3 = - M1 - M2
!
!
  OLDFAC = NEWFAC
  A1 = (L2-M2+ONE) * (L2+M2) * (L3+M3+ONE) * (L3-M3)
  NEWFAC = SQRT(A1)
!
!
  DV = (L1+L2+L3+ONE)*(L2+L3-L1) - (L2-M2+ONE)*(L3+M3+ONE) &
                                 - (L2+M2-ONE)*(L3-M3-ONE)
!
  if ( LSTEP-2)  32, 32, 31
!
   31 C1OLD = ABS(C1)
   32 C1 = - DV / NEWFAC
!
  if ( LSTEP > 2)   go to 60
!
!
!  If M2 = M2MIN + 1, the third term in the recursion equation vanishes,
!  hence
!
  X = SRTINY * C1
  THRCOF(2) = X
  SUM1 = SUM1 + TINY * C1*C1
  if ( LSTEP == NFIN)   go to 220
  go to 30
!
!
   60 C2 = - OLDFAC / NEWFAC
!
!  Recursion to the next 3j coefficient
  X = C1 * THRCOF(LSTEP-1) + C2 * THRCOF(LSTEP-2)
  THRCOF(LSTEP) = X
  SUMFOR = SUM1
  SUM1 = SUM1 + X*X
  if ( LSTEP == NFIN)   go to 100
!
!  See if last unnormalized 3j coefficient exceeds SRHUGE
!
  if ( ABS(X) < SRHUGE)   go to 80
!
!  This is reached if last 3j coefficient larger than SRHUGE,
!  so that the recursion series THRCOF(1), ... , THRCOF(LSTEP)
!  has to be rescaled to prevent overflow
!
!     MSCALE = MSCALE + 1
  DO 70 I=1,LSTEP
  if ( ABS(THRCOF(I)) < SRTINY)   THRCOF(I) = ZERO
   70 THRCOF(I) = THRCOF(I) / SRHUGE
  SUM1 = SUM1 / HUGE
  SUMFOR = SUMFOR / HUGE
  X = X / SRHUGE
!
!
!  As long as ABS(C1) is decreasing, the recursion proceeds towards
!  increasing 3j values and, hence, is numerically stable.  Once
!  an increase of ABS(C1) is detected, the recursion direction is
!  reversed.
!
   80 if ( C1OLD-ABS(C1))   100, 100, 30
!
!
!  Keep three 3j coefficients around MMATCH for comparison later
!  with backward recursion values.
!
  100 CONTINUE
!     MMATCH = M2 - 1
  NSTEP2 = NFIN - LSTEP + 3
  X1 = X
  X2 = THRCOF(LSTEP-1)
  X3 = THRCOF(LSTEP-2)
!
!  Starting backward recursion from M2MAX taking NSTEP2 steps, so
!  that forwards and backwards recursion overlap at the three points
!  M2 = MMATCH+1, MMATCH, MMATCH-1.
!
  NFINP1 = NFIN + 1
  NFINP2 = NFIN + 2
  NFINP3 = NFIN + 3
  THRCOF(NFIN) = SRTINY
  SUM2 = TINY
!
!
!
  M2 = M2MAX + TWO
  LSTEP = 1
  110 LSTEP = LSTEP + 1
  M2 = M2 - ONE
  M3 = - M1 - M2
  OLDFAC = NEWFAC
  A1S = (L2-M2+TWO) * (L2+M2-ONE) * (L3+M3+TWO) * (L3-M3-ONE)
  NEWFAC = SQRT(A1S)
  DV = (L1+L2+L3+ONE)*(L2+L3-L1) - (L2-M2+ONE)*(L3+M3+ONE) &
                                 - (L2+M2-ONE)*(L3-M3-ONE)
  C1 = - DV / NEWFAC
  if ( LSTEP > 2)   go to 120
!
!  If M2 = M2MAX + 1 the third term in the recursion equation vanishes
!
  Y = SRTINY * C1
  THRCOF(NFIN-1) = Y
  if ( LSTEP == NSTEP2)   go to 200
  SUMBAC = SUM2
  SUM2 = SUM2 + Y*Y
  go to 110
!
  120 C2 = - OLDFAC / NEWFAC
!
!  Recursion to the next 3j coefficient
!
  Y = C1 * THRCOF(NFINP2-LSTEP) + C2 * THRCOF(NFINP3-LSTEP)
!
  if ( LSTEP == NSTEP2)   go to 200
!
  THRCOF(NFINP1-LSTEP) = Y
  SUMBAC = SUM2
  SUM2 = SUM2 + Y*Y
!
!
!  See if last 3j coefficient exceeds SRHUGE
!
  if ( ABS(Y) < SRHUGE)   go to 110
!
!  This is reached if last 3j coefficient larger than SRHUGE,
!  so that the recursion series THRCOF(NFIN), ... , THRCOF(NFIN-LSTEP+1)
!  has to be rescaled to prevent overflow.
!
!     MSCALE = MSCALE + 1
  DO 111 I=1,LSTEP
  INDEX = NFIN - I + 1
  if ( ABS(THRCOF(INDEX)) < SRTINY) &
    THRCOF(INDEX) = ZERO
  111 THRCOF(INDEX) = THRCOF(INDEX) / SRHUGE
  SUM2 = SUM2 / HUGE
  SUMBAC = SUMBAC / HUGE
!
  go to 110
!
!
!
!  The forward recursion 3j coefficients X1, X2, X3 are to be matched
!  with the corresponding backward recursion values Y1, Y2, Y3.
!
  200 Y3 = Y
  Y2 = THRCOF(NFINP2-LSTEP)
  Y1 = THRCOF(NFINP3-LSTEP)
!
!
!  Determine now RATIO such that YI = RATIO * XI  (I=1,2,3) holds
!  with minimal error.
!
  RATIO = ( X1*Y1 + X2*Y2 + X3*Y3 ) / ( X1*X1 + X2*X2 + X3*X3 )
  NLIM = NFIN - NSTEP2 + 1
!
  if ( ABS(RATIO) < ONE)   go to 211
!
  DO 210 N=1,NLIM
  210 THRCOF(N) = RATIO * THRCOF(N)
  SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC
  go to 230
!
  211 NLIM = NLIM + 1
  RATIO = ONE / RATIO
  DO 212 N=NLIM,NFIN
  212 THRCOF(N) = RATIO * THRCOF(N)
  SUMUNI = SUMFOR + RATIO*RATIO*SUMBAC
  go to 230
!
  220 SUMUNI = SUM1
!
!
!  Normalize 3j coefficients
!
  230 CNORM = ONE / SQRT((L1+L1+ONE) * SUMUNI)
!
!  Sign convention for last 3j coefficient determines overall phase
!
  SIGN1 = SIGN(ONE,THRCOF(NFIN))
  SIGN2 = (-ONE) ** INT(ABS(L2-L3-M1)+EPS)
  if ( SIGN1*SIGN2)  235,235,236
  235 CNORM = - CNORM
!
  236 if ( ABS(CNORM) < ONE)   go to 250
!
  DO 240 N=1,NFIN
  240 THRCOF(N) = CNORM * THRCOF(N)
  return
!
  250 THRESH = TINY / ABS(CNORM)
  DO 251 N=1,NFIN
  if ( ABS(THRCOF(N)) < THRESH)   THRCOF(N) = ZERO
  251 THRCOF(N) = CNORM * THRCOF(N)
!
!
!
  return
end


!DECK D1MACH
      DOUBLE PRECISION FUNCTION D1MACH (I)
      IMPLICIT NONE
      INTEGER :: I
      DOUBLE PRECISION :: B, X
!***BEGIN PROLOGUE  D1MACH
!***PURPOSE  Return floating point machine dependent constants.
!***LIBRARY   SLATEC
!***CATEGORY  R1
!***TYPE      SINGLE PRECISION (D1MACH-S, D1MACH-D)
!***KEYWORDS  MACHINE CONSTANTS
!***AUTHOR  Fox, P. A., (Bell Labs)
!           Hall, A. D., (Bell Labs)
!           Schryer, N. L., (Bell Labs)
!***DESCRIPTION
!
!   D1MACH can be used to obtain machine-dependent parameters for the
!   local machine environment.  It is a function subprogram with one
!   (input) argument, and can be referenced as follows:
!
!        A = D1MACH(I)
!
!   where I=1,...,5.  The (output) value of A above is determined by
!   the (input) value of I.  The results for various values of I are
!   discussed below.
!
!   D1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
!   D1MACH(2) = B**EMAX*(1 - B**(-T)), the largest magnitude.
!   D1MACH(3) = B**(-T), the smallest relative spacing.
!   D1MACH(4) = B**(1-T), the largest relative spacing.
!   D1MACH(5) = LOG10(B)
!
!   Assume single precision numbers are represented in the T-digit,
!   base-B form
!
!              sign (B**E)*( (X(1)/B) + ... + (X(T)/B**T) )
!
!   where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
!   EMIN .LE. E .LE. EMAX.
!
!   The values of B, T, EMIN and EMAX are provided in I1MACH as
!   follows:
!   I1MACH(10) = B, the base.
!   I1MACH(11) = T, the number of base-B digits.
!   I1MACH(12) = EMIN, the smallest exponent E.
!   I1MACH(13) = EMAX, the largest exponent E.
!
!
!***REFERENCES  P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
!                 a portable library, ACM Transactions on Mathematical
!                 Software 4, 2 (June 1978), pp. 177-188.
!***ROUTINES CALLED  XERMSG
!***REVISION HISTORY  (YYMMDD)
!   790101  DATE WRITTEN
!   960329  Modified for Fortran 90 (BE after suggestions by EHG)      
!***END PROLOGUE  D1MACH
!      
      X = 1.0D0
      B = RADIX(X)
      SELECT CASE (I)
        CASE (1)
          D1MACH = B**(MINEXPONENT(X)-1) ! the smallest positive magnitude.
        CASE (2)
          D1MACH = HUGE(X)               ! the largest magnitude.
        CASE (3)
          D1MACH = B**(-DIGITS(X))       ! the smallest relative spacing.
        CASE (4)
          D1MACH = B**(1-DIGITS(X))      ! the largest relative spacing.
        CASE (5)
          D1MACH = LOG10(B)
        CASE DEFAULT
          WRITE (*, FMT = 9000)
 9000     FORMAT ('1ERROR    1 IN D1MACH - I OUT OF BOUNDS')
          STOP
      END SELECT
      RETURN
      END
Revision:	1000
Committed:	Mon Jul 3 19:40:52 2006 UTC (18 years, 10 months ago) by chuckv
File size:	15770 byte(s)
Log Message:	Added utility function to winger3jm from slatac.
#	User	Rev	Content
1	chuckv	991	subroutine WIGNER3JM (L1, L2, L3, M1, M2MIN, M2MAX, THRCOF, NDIM, &
2			IER)
3	chuckv	1000	! use status
4	chuckv	991	!
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	chuckv	1000	! call handleError('Wigner3jm','L1-ABS(M1) less than zero or '// &
162			! 'L1+ABS(M1) not integer.')
163	chuckv	991	return
164			ELSEIF((L1+L2-L3 < -EPS).OR.(L1-L2+L3 < -EPS).OR. &
165			(-L1+L2+L3 < -EPS))THEN
166			IER=2
167	chuckv	1000	write(,) eps
168			write(,) L1,L2,L3
169			write(,) "L1,L2,L3 do not satisfy triangular condition"
170			! call handleError('Wigner3jm','L1, L2, L3 do not satisfy '// &
171			! 'triangular condition.')
172	chuckv	991	return
173			ELSEIF(MOD(L1+L2+L3+EPS,ONE) >= EPS+EPS)THEN
174			IER=3
175	chuckv	1000	! call handleError('Wigner3jm','L1+L2+L3 not integer.')
176	chuckv	991	return
177			end if
178			!
179			!
180			! Limits for M2
181			M2MIN = MAX(-L2,-L3-M1)
182			M2MAX = MIN(L2,L3-M1)
183			!
184			! Check error condition 4.
185			if ( MOD(M2MAX-M2MIN+EPS,ONE) >= EPS+EPS)THEN
186			IER=4
187	chuckv	1000	! call handleError('Wigner3jm','M2MAX-M2MIN not integer.')
188	chuckv	991	return
189			end if
190			if ( M2MIN < M2MAX-EPS) go to 20
191			if ( M2MIN < M2MAX+EPS) go to 10
192			!
193			! Check error condition 5.
194			IER=5
195	chuckv	1000	! call handleError('Wigner3jm','M2MIN greater than M2MAX.')
196	chuckv	991	return
197			!
198			!
199			! This is reached in case that M2 and M3 can take only one value.
200			10 CONTINUE
201			! MSCALE = 0
202			THRCOF(1) = (-ONE) ** INT(ABS(L2-L3-M1)+EPS) / &
203			SQRT(L1+L2+L3+ONE)
204			return
205			!
206			! This is reached in case that M1 and M2 take more than one value.
207			20 CONTINUE
208			! MSCALE = 0
209			NFIN = INT(M2MAX-M2MIN+ONE+EPS)
210			if ( NDIM-NFIN) 21, 23, 23
211			!
212			! Check error condition 6.
213			21 IER = 6
214	chuckv	1000	! call handleError('Wigner3jm','Dimension of result array for '// &
215			! '3j coefficients too small.')
216	chuckv	991	return
217			!
218			!
219			!
220			! Start of forward recursion from M2 = M2MIN
221			!
222			23 M2 = M2MIN
223			THRCOF(1) = SRTINY
224			NEWFAC = 0.0D0
225			C1 = 0.0D0
226			SUM1 = TINY
227			!
228			!
229			LSTEP = 1
230			30 LSTEP = LSTEP + 1
231			M2 = M2 + ONE
232			M3 = - M1 - M2
233			!
234			!
235			OLDFAC = NEWFAC
236			A1 = (L2-M2+ONE) * (L2+M2) * (L3+M3+ONE) * (L3-M3)
237			NEWFAC = SQRT(A1)
238			!
239			!
240			DV = (L1+L2+L3+ONE)(L2+L3-L1) - (L2-M2+ONE)(L3+M3+ONE) &
241			- (L2+M2-ONE)*(L3-M3-ONE)
242			!
243			if ( LSTEP-2) 32, 32, 31
244			!
245			31 C1OLD = ABS(C1)
246			32 C1 = - DV / NEWFAC
247			!
248			if ( LSTEP > 2) go to 60
249			!
250			!
251			! If M2 = M2MIN + 1, the third term in the recursion equation vanishes,
252			! hence
253			!
254			X = SRTINY * C1
255			THRCOF(2) = X
256			SUM1 = SUM1 + TINY * C1*C1
257			if ( LSTEP == NFIN) go to 220
258			go to 30
259			!
260			!
261			60 C2 = - OLDFAC / NEWFAC
262			!
263			! Recursion to the next 3j coefficient
264			X = C1 * THRCOF(LSTEP-1) + C2 * THRCOF(LSTEP-2)
265			THRCOF(LSTEP) = X
266			SUMFOR = SUM1
267			SUM1 = SUM1 + X*X
268			if ( LSTEP == NFIN) go to 100
269			!
270			! See if last unnormalized 3j coefficient exceeds SRHUGE
271			!
272			if ( ABS(X) < SRHUGE) go to 80
273			!
274			! This is reached if last 3j coefficient larger than SRHUGE,
275			! so that the recursion series THRCOF(1), ... , THRCOF(LSTEP)
276			! has to be rescaled to prevent overflow
277			!
278			! MSCALE = MSCALE + 1
279			DO 70 I=1,LSTEP
280			if ( ABS(THRCOF(I)) < SRTINY) THRCOF(I) = ZERO
281			70 THRCOF(I) = THRCOF(I) / SRHUGE
282			SUM1 = SUM1 / HUGE
283			SUMFOR = SUMFOR / HUGE
284			X = X / SRHUGE
285			!
286			!
287			! As long as ABS(C1) is decreasing, the recursion proceeds towards
288			! increasing 3j values and, hence, is numerically stable. Once
289			! an increase of ABS(C1) is detected, the recursion direction is
290			! reversed.
291			!
292			80 if ( C1OLD-ABS(C1)) 100, 100, 30
293			!
294			!
295			! Keep three 3j coefficients around MMATCH for comparison later
296			! with backward recursion values.
297			!
298			100 CONTINUE
299			! MMATCH = M2 - 1
300			NSTEP2 = NFIN - LSTEP + 3
301			X1 = X
302			X2 = THRCOF(LSTEP-1)
303			X3 = THRCOF(LSTEP-2)
304			!
305			! Starting backward recursion from M2MAX taking NSTEP2 steps, so
306			! that forwards and backwards recursion overlap at the three points
307			! M2 = MMATCH+1, MMATCH, MMATCH-1.
308			!
309			NFINP1 = NFIN + 1
310			NFINP2 = NFIN + 2
311			NFINP3 = NFIN + 3
312			THRCOF(NFIN) = SRTINY
313			SUM2 = TINY
314			!
315			!
316			!
317			M2 = M2MAX + TWO
318			LSTEP = 1
319			110 LSTEP = LSTEP + 1
320			M2 = M2 - ONE
321			M3 = - M1 - M2
322			OLDFAC = NEWFAC
323			A1S = (L2-M2+TWO) * (L2+M2-ONE) * (L3+M3+TWO) * (L3-M3-ONE)
324			NEWFAC = SQRT(A1S)
325			DV = (L1+L2+L3+ONE)(L2+L3-L1) - (L2-M2+ONE)(L3+M3+ONE) &
326			- (L2+M2-ONE)*(L3-M3-ONE)
327			C1 = - DV / NEWFAC
328			if ( LSTEP > 2) go to 120
329			!
330			! If M2 = M2MAX + 1 the third term in the recursion equation vanishes
331			!
332			Y = SRTINY * C1
333			THRCOF(NFIN-1) = Y
334			if ( LSTEP == NSTEP2) go to 200
335			SUMBAC = SUM2
336			SUM2 = SUM2 + Y*Y
337			go to 110
338			!
339			120 C2 = - OLDFAC / NEWFAC
340			!
341			! Recursion to the next 3j coefficient
342			!
343			Y = C1 * THRCOF(NFINP2-LSTEP) + C2 * THRCOF(NFINP3-LSTEP)
344			!
345			if ( LSTEP == NSTEP2) go to 200
346			!
347			THRCOF(NFINP1-LSTEP) = Y
348			SUMBAC = SUM2
349			SUM2 = SUM2 + Y*Y
350			!
351			!
352			! See if last 3j coefficient exceeds SRHUGE
353			!
354			if ( ABS(Y) < SRHUGE) go to 110
355			!
356			! This is reached if last 3j coefficient larger than SRHUGE,
357			! so that the recursion series THRCOF(NFIN), ... , THRCOF(NFIN-LSTEP+1)
358			! has to be rescaled to prevent overflow.
359			!
360			! MSCALE = MSCALE + 1
361			DO 111 I=1,LSTEP
362			INDEX = NFIN - I + 1
363			if ( ABS(THRCOF(INDEX)) < SRTINY) &
364			THRCOF(INDEX) = ZERO
365			111 THRCOF(INDEX) = THRCOF(INDEX) / SRHUGE
366			SUM2 = SUM2 / HUGE
367			SUMBAC = SUMBAC / HUGE
368			!
369			go to 110
370			!
371			!
372			!
373			! The forward recursion 3j coefficients X1, X2, X3 are to be matched
374			! with the corresponding backward recursion values Y1, Y2, Y3.
375			!
376			200 Y3 = Y
377			Y2 = THRCOF(NFINP2-LSTEP)
378			Y1 = THRCOF(NFINP3-LSTEP)
379			!
380			!
381			! Determine now RATIO such that YI = RATIO * XI (I=1,2,3) holds
382			! with minimal error.
383			!
384			RATIO = ( X1Y1 + X2Y2 + X3Y3 ) / ( X1X1 + X2X2 + X3X3 )
385			NLIM = NFIN - NSTEP2 + 1
386			!
387			if ( ABS(RATIO) < ONE) go to 211
388			!
389			DO 210 N=1,NLIM
390			210 THRCOF(N) = RATIO * THRCOF(N)
391			SUMUNI = RATIO * RATIO * SUMFOR + SUMBAC
392			go to 230
393			!
394			211 NLIM = NLIM + 1
395			RATIO = ONE / RATIO
396			DO 212 N=NLIM,NFIN
397			212 THRCOF(N) = RATIO * THRCOF(N)
398			SUMUNI = SUMFOR + RATIORATIOSUMBAC
399			go to 230
400			!
401			220 SUMUNI = SUM1
402			!
403			!
404			! Normalize 3j coefficients
405			!
406			230 CNORM = ONE / SQRT((L1+L1+ONE) * SUMUNI)
407			!
408			! Sign convention for last 3j coefficient determines overall phase
409			!
410			SIGN1 = SIGN(ONE,THRCOF(NFIN))
411			SIGN2 = (-ONE) ** INT(ABS(L2-L3-M1)+EPS)
412			if ( SIGN1*SIGN2) 235,235,236
413			235 CNORM = - CNORM
414			!
415			236 if ( ABS(CNORM) < ONE) go to 250
416			!
417			DO 240 N=1,NFIN
418			240 THRCOF(N) = CNORM * THRCOF(N)
419			return
420			!
421			250 THRESH = TINY / ABS(CNORM)
422			DO 251 N=1,NFIN
423			if ( ABS(THRCOF(N)) < THRESH) THRCOF(N) = ZERO
424			251 THRCOF(N) = CNORM * THRCOF(N)
425			!
426			!
427			!
428			return
429			end
430	chuckv	1000
431
432			!DECK D1MACH
433			DOUBLE PRECISION FUNCTION D1MACH (I)
434			IMPLICIT NONE
435			INTEGER :: I
436			DOUBLE PRECISION :: B, X
437			!***BEGIN PROLOGUE D1MACH
438			!***PURPOSE Return floating point machine dependent constants.
439			!***LIBRARY SLATEC
440			!***CATEGORY R1
441			!***TYPE SINGLE PRECISION (D1MACH-S, D1MACH-D)
442			!***KEYWORDS MACHINE CONSTANTS
443			!***AUTHOR Fox, P. A., (Bell Labs)
444			! Hall, A. D., (Bell Labs)
445			! Schryer, N. L., (Bell Labs)
446			!***DESCRIPTION
447			!
448			! D1MACH can be used to obtain machine-dependent parameters for the
449			! local machine environment. It is a function subprogram with one
450			! (input) argument, and can be referenced as follows:
451			!
452			! A = D1MACH(I)
453			!
454			! where I=1,...,5. The (output) value of A above is determined by
455			! the (input) value of I. The results for various values of I are
456			! discussed below.
457			!
458			! D1MACH(1) = B**(EMIN-1), the smallest positive magnitude.
459			! D1MACH(2) = B*EMAX(1 - B**(-T)), the largest magnitude.
460			! D1MACH(3) = B**(-T), the smallest relative spacing.
461			! D1MACH(4) = B**(1-T), the largest relative spacing.
462			! D1MACH(5) = LOG10(B)
463			!
464			! Assume single precision numbers are represented in the T-digit,
465			! base-B form
466			!
467			! sign (B*E)( (X(1)/B) + ... + (X(T)/B**T) )
468			!
469			! where 0 .LE. X(I) .LT. B for I=1,...,T, 0 .LT. X(1), and
470			! EMIN .LE. E .LE. EMAX.
471			!
472			! The values of B, T, EMIN and EMAX are provided in I1MACH as
473			! follows:
474			! I1MACH(10) = B, the base.
475			! I1MACH(11) = T, the number of base-B digits.
476			! I1MACH(12) = EMIN, the smallest exponent E.
477			! I1MACH(13) = EMAX, the largest exponent E.
478			!
479			!
480			!***REFERENCES P. A. Fox, A. D. Hall and N. L. Schryer, Framework for
481			! a portable library, ACM Transactions on Mathematical
482			! Software 4, 2 (June 1978), pp. 177-188.
483			!***ROUTINES CALLED XERMSG
484			!***REVISION HISTORY (YYMMDD)
485			! 790101 DATE WRITTEN
486			! 960329 Modified for Fortran 90 (BE after suggestions by EHG)
487			!***END PROLOGUE D1MACH
488			!
489			X = 1.0D0
490			B = RADIX(X)
491			SELECT CASE (I)
492			CASE (1)
493			D1MACH = B**(MINEXPONENT(X)-1) ! the smallest positive magnitude.
494			CASE (2)
495			D1MACH = HUGE(X) ! the largest magnitude.
496			CASE (3)
497			D1MACH = B**(-DIGITS(X)) ! the smallest relative spacing.
498			CASE (4)
499			D1MACH = B**(1-DIGITS(X)) ! the largest relative spacing.
500			CASE (5)
501			D1MACH = LOG10(B)
502			CASE DEFAULT
503			WRITE (*, FMT = 9000)
504			9000 FORMAT ('1ERROR 1 IN D1MACH - I OUT OF BOUNDS')
505			STOP
506			END SELECT
507			RETURN
508			END