wignerSymbols-fortran.f

C******************************************************-/
C This source code is subject to the terms of the GNU  -/
C Lesser Public License. If a copy of the LGPL was not -/
C distributed with this file, you can obtain one at    -/
C https://www.gnu.org/licenses/lgpl.html.              -/
C ******************************************************/ 


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