C N. R. BADNELL                                                 22/02/21
C
      PROGRAM TWIG3NJ
C
C-----------------------------------------------------------------------
C TEST
C EVALUATION OF WIGNER 3N-J SYMBOL BY SUMMATION OF 6-J N-PRODUCT -WIG6JR
C SEE E.G. A. R. EDMONDS "ANGULAR MOMENTUM IN QUANTUM MECHANICS" (1957)
C FOR N=2,3 AND ORD-SMITH PHYS.REV.94 1227 (1954) FOR N>3 (& BY ANALOGY)
C
C THESE FORMS OF 3N-J COEFFICIENTS ARE ALSO EQUIVALENT TO THE FIRST KIND
C OF 3N-J COEFFICIENTS DEFINED  BY YUTSIS, LEVINSON & VANAGAS (YLV)
C EQS (17.1) AND (17.2) "MATHEMATICAL APPARATUS OF THE THEORY OF ANGULAR
C MOMENTUM" ISRAEL PROGRAM FOR SCIENTIFIC TRANSLATIONS: JERUSALEM (1962)
C THEY JUST HAVE A DIFFERENT ARRANGEMENT OF ARGUMENTS (EXCEPT FOR N=2
C WHERE YLV HAS THE PHASE OF THE RACAH W-COEFFICIENT NOT WIGNER 6-J.)
C
C
C EXAMPLES:
C
C N=2::                      TWIG6JR EG 3
C  20  32  28
C  26  30  30
C
C N=3:                       TWIG9J  EG 2
C 200 160 100
C 100 200 140
C 120 100 200
C
C N=4:                       WEI J.PHYS.A40 229501 (2013) WIGNER ORDER
C  36  18  30  40
C  20  30  24  18
C  16  20  28  30
C
C N=5                        YU PH.D. THESIS BERKELEY (2010)
C 197 187 148 150 148
C 173 205 176 192 194
C 190 180   2   2 160
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (IWRITE=0)              !=0 TO SCREEN, >0 TO WIG3NJ.out
C
      ALLOCATABLE :: A(:),B(:),C(:)
C
C
      IF(IWRITE.GT.0)OPEN(IWRITE,FILE='C3NJ.out')      !STATUS='UNKNOWN'
C
C ASK USER TO ENTER N AND ALL 3*N PREDEFINED ARGUMENTS (ROWWISE)
C
      WRITE(*,*)"DETERMINES WIGNER 3N-J SYMBOL FOR GIVEN "
     X         ,"3*N ARGUMENTS (rowwise)"
C
   1  WRITE(IWRITE,*)
      WRITE(*,*)"ENTER N (.eq.0 exits)"
C
      READ(*,*)N0
C
      N=ABS(N0)
      IF(N.EQ.0.OR.N0.EQ.-1)THEN
        WRITE(*,*)'*** EXITING AS N.eq.0,-1'
        IF(IWRITE.GT.0)WRITE(*,*)'*** 3N-J OUTPUT IN WIG3NJ.out'
        STOP
      ELSE
        N3=MAX(N,3)
        ALLOCATE (A(N3),B(N3),C(N3))
      ENDIF
C
   2  WRITE(IWRITE,*)
      IF(N.EQ.1)THEN
        WRITE(*,*)"ENTER THE 6 INTEGERS (ROWWISE) AT *TWICE* "
     X           ,"ACTUAL VALUE (any.lt.0 exits):"
      ELSE
        WRITE(*,*)"ENTER THE 3*N INTEGERS (ROWWISE) AT *TWICE* "
     X           ,"ACTUAL VALUE (any.lt.0 exits):"
      ENDIF
      IF(N0.GT.0)THEN
C        WRITE(*,*)"(DO *NOT* USE YLV ORDERING!)"
      ELSE
        WRITE(*,*)"(N.B. YOU HAVE FLAGGED INPUT AS YLV ORDERED)"
      ENDIF
C
      IF(N.LT.3)THEN
        READ(*,*)(A(I),I=1,3),(B(I),I=1,3)
        DO I=1,3
          C(I)=0
        ENDDO
      ELSE
        READ(*,*)(A(I),I=1,N),(B(I),I=1,N),(C(I),I=1,N)
      ENDIF
C
      DO I=1,N3
        IF(A(I).LT.0)THEN
          WRITE(*,*)'*** EXITING AS ARGUMENTS ARE NOT ALL .GE. ZERO'
          IF(IWRITE.GT.0)WRITE(*,*)'*** 3N-J OUTPUT IN WIG3NJ.out'
          DEALLOCATE (A,B,C)
          GO TO 1
        ENDIF
        IF(N.GT.1)THEN
          IF(B(I).LT.0.OR.C(I).LT.0)THEN
            WRITE(*,*)'*** EXITING AS ARGUMENTS ARE NOT ALL .GE. ZERO'
            IF(IWRITE.GT.0)WRITE(*,*)'*** 3N-J OUTPUT IN WIG3NJ.out'
            DEALLOCATE (A,B,C)
            GO TO 1
          ENDIF
        ENDIF
        A(I)=A(I)/2
        B(I)=B(I)/2
        C(I)=C(I)/2
      ENDDO
C
C CONVERT FROM YLV62 ORDERED ARGUMENTS TO WIGNER (FOR N.GT.1)
C NOTE PHASE FOR CASE N=2, WIGNER 6-J IS OUTPUT.
C
      IF(N0.LT.0)THEN
        IF(N.EQ.2)THEN                    !ALIGN FOR YLV2WIG TO RE-ORDER
          C(2)=B(3)
          C(1)=B(2)
          B(2)=B(1)
          B(1)=A(3)
        ENDIF
        CALL YLV2WIG(N,A,B,C)
      ENDIF
c
ct       call cpu_time(timei)
C
C
ct      do i=1,10000
      CALL WIG3NJ(W3NJ,N,A,B,C,KMIN,KMAX,IH)
ct      enddo
c
ct       call cpu_time(timef)
ct       times=timef-timei
ct       write(*,"(' recursion time=',1pe9.1,' msecs')")times*1000
C
      IF(KMIN.LT.0)THEN
        IF(KMIN.EQ.-3)THEN                            !SHOULD NOT HAPPEN
          WRITE(*,*)'*** DIMENSION FAILURE, INCREASE NDIMW TO ',KMAX
        ELSE
          WRITE(*,*)'*** SELECTION RULE FAILURE, ALL 3N-J EQUAL ZERO!'
          IF(KMIN.EQ.-1)WRITE(*,*)'*** non-integer triangle sum'
          IF(KMIN.EQ.-2)WRITE(*,*)'*** triangle rule failure'
        ENDIF
        GO TO 2
      ENDIF
C
      IF(IWRITE.GT.0)WRITE(*,"('W3NJ=',1P2D25.16)")W3NJ
C
      IF(IWRITE.GE.0)THEN
        WRITE(IWRITE,"(/4X,'W3NJ(A,B,C)',3X,'N=',I3,', 2A 2B 2C =')")N
        IF(N0.LT.0)
     X WRITE(IWRITE,"(25X,'OUTPUT ARGUMENTS ARE STILL WIGNER-ORDERED')")
        WRITE(IWRITE,"(32X,10I7)")(NINT(2*A(I)),I=1,N3)
        WRITE(IWRITE,"(32X,10I7)")(NINT(2*B(I)),I=1,N3)
        IF(N.GT.2)WRITE(IWRITE,"(32X,10I7)")(NINT(2*C(I)),I=1,N3)
        WRITE(IWRITE,"(10X,1P2D25.16)")W3NJ
        WRITE(IWRITE,*)
      ENDIF
c
      if(iwrite.gt.0)call flush(IWRITE)
C
      GO TO 2
C
      END PROGRAM TWIG3NJ
C
C***********************************************************************
C
      SUBROUTINE WIG3NJ(W3NJ,NMAX,A,B,C,KMIN,KMAX,IH)
C
C-----------------------------------------------------------------------
C
C N.R.BADNELL
C
C CALCULATES THE WIGNER 3N-J SYMBOL ( A(N), N=1,NMAX )
C                                   ( B(N), N=1,NMAX )
C                                   ( C(N), N=1,NMAX )
C
C CORRESPONDING TO THE COUPLING OF N+1 ANGULAR MOMENTA:
C
C (  J1  J2 J12 J125 J1256 ... J125...NMAX+1 )
C (  J3  J4 J34 J135 J1356 ... J135...NMAX+1 )
C ( J13 J24  J5   J6    J7 ...       JNMAX+2 )
C
C FOR N=NMAX.GE.3
C
C IN TERMS OF A SUM OF AN N-PRODUCT OF 6-J SYMBOLS USING SR.WIG6JR.
C SEE E.G. A. R. EDMONDS "ANGULAR MOMENTUM IN QUANTUM MECHANICS" (1957)
C FOR N=2,3 AND ORD-SMITH PHYS.REV.94 1227 (1954) FOR N>3 (& BY ANALOGY)
C SEE ALSO YUTSIS, LEVINSON & VANAGAS (1962) EQS (17.1) AND (17.2)
C USES WIGNER 6-J RECURRENCE OF BADNELL AT AL MNRAS (2021).
C
C FOR NMAX=1 OR 2, JUST CALLS THE APPROPRIATE WIG3JR/6JR.
C
C INPUT:
C  NMAX, THE N-VALUE FOR WIGNER 3N-J COEFFICIENTS TO BE EVALUATED.
C  A(N),B(N),C(N), N=1,NMAX: ALL 3*N ARGUMENTS OF THE 3N-J SYMBOL, AS
C  DEFINED ABOVE, USING THEIR ACTUAL VALUE, I.E. *NOT* TWICE  (REAL*8).
C  (THESE MUST BE IN WIGNER ORDERING, NOT YUTSIS, LEVINSON & VANAGAS.)
C  (MAKE A CALL TO SR.YLVWIG BEFOREHAND, IF NECESSARY, TO CONVERT.)
C
C OUTPUT:
C  W3NJ: THE REQUIRED WIGNER 3N-J SYMBOL (REAL*8)
C  KMIN,KMAX,IH: THE RANGE (KMIN+IH/2,KMAX+IH/2)OF THE K-SUM OVER THE
C       6-J SYMBOL N-PRODUCT (INTEGER)
C
C  KMIN.LT.0 FLAGS AN ERROR:
C         -1 NON-INTEGER TRIANGLE SUM.
C         -2 TRIANGLE RULE FAILURE.
C         -3 DIMENSION ERROR, KMAX THEN CONTAINS THE REQUIRED NDIMW.
C
C UPDATE LOG:
C 23/12/15 - PHASE WAS NOT DETERMINED IF W6J(JMAX)=0
C 11/12/15 - REMOVED TEMP ARRAY.
C 20/11/15 - INITIAL RELEASE.
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      logical debug1,debug2
C
      PARAMETER (DZERO=0.0D0)
      PARAMETER (DHALF=0.5D0)
      PARAMETER (DONE=1.0D0)
      PARAMETER (DEPS1=1.D-30)                  !SMALLEST NON-ZERO COEFF
      PARAMETER (DEPS2=1.D-5)
C
      DIMENSION A(*),B(*),C(*)               !* NOT NMAX, FOR CASE N=1,2
C
      ALLOCATABLE :: W6J(:),W1(:)
C
      data debug1/.true./,debug2/.true./
      data iwritd/0/     !=0 to screen; set < 0 for non-interactive use
C
      if(iwritd.lt.0)debug1=.false.
C
C INITIALIZE (LEAVE 3N-J SELECTION RULES TO WIG6JR)
C
      W3NJ=DZERO
      IH=0
C
C CATCH SIMPLE CASES (N=1 & 2):
C
      IF(NMAX.EQ.1)THEN
C
        IF(NINT(-B(2)-B(3)-B(2)-B(3)).NE.NINT(B(1)+B(1)))THEN
          if(debug2)write(*,*)'*** sr.wig3njr: triangle rule failure...'
          KMIN=-2
          RETURN
        ENDIF
C
        NDIMW=NINT(A(2)+A(3)-DEPS2)
C
        ALLOCATE (W6J(-1:NDIMW+1))
C
        CALL WIG3JRJ(W6J,A(2),A(3),B(2),B(3),JMIN,JMAX,IH,NDIMW)
C
        IF(JMIN.LT.0)GO TO 500
C
        KMIN=JMIN
        KMAX=JMAX
C
        IA1=NINT(A(1)-DEPS2)
        W3NJ=W6J(IA1)
C
        GO TO 100
C
      ENDIF
C
      IF(NMAX.EQ.2)THEN
C
        T=MIN(A(2)+A(3),B(2)+B(3))
        NDIMW=NINT(T-DEPS2)
C
        ALLOCATE (W6J(-1:NDIMW+1))
C
        CALL WIG6JR(W6J,A(2),A(3),B(1),B(2),B(3),JMIN,JMAX,IH,NDIMW)
C
        IF(JMIN.LT.0)GO TO 500
C
        KMIN=JMIN
        KMAX=JMAX
C
        IA1=NINT(A(1)-DEPS2)
        W3NJ=W6J(IA1)
C
        GO TO 100
C
      ENDIF
C
C GENERAL CASE N.GE.3
C
C DETERMINE 6J SUMMATION RANGE K=KMIN,KMAX (QUICK EXIT IF ZERO SUM)
C
      T=MAX(ABS(A(3)-C(1)),ABS(A(2)-B(1)),ABS(C(2)-B(3)))
      DO N=4,NMAX
        T=MAX(T,ABS(A(N)-B(N)))
      ENDDO
      KMIN=NINT(T-DEPS2)
C
      T=MIN(A(3)+C(1),A(2)+B(1),C(2)+B(3))
      DO N=4,NMAX
        T=MIN(T,A(N)+B(N))
      ENDDO
      KMAX=NINT(T-DEPS2)
C
      IF(KMIN.GT.KMAX)THEN
        KMIN=-2
        RETURN
      ENDIF
C
      ALLOCATE (W1(KMIN:KMAX))
C
      IF(KMAX.NE.NINT(T+DEPS2))IH=1                        !HALF-INTEGER
C
      th=dhalf*ih
      if(debug1)
     x  write(iwritd,"(/' 6j sum range:   ',2f9.1)")kmin+th,kmax+th
C
C KMAX IS THE COMMON RANGE, WE NEED THE MAX RANGE OF THE N SEPARATE W6J.
C (NOT WORTH N SEPARATE ALLOCATES WITH BESPOKE NDIMW.)
C
      T=DZERO
      T0=MIN(A(3)+C(1),B(1)+A(2))
      T=MAX(T,T0)
      T0=MIN(A(2)+B(1),B(3)+C(2))
      T=MAX(T,T0)
      IF(NMAX.EQ.3)THEN
        T0=MIN(C(2)+B(3),A(3)+C(1))
        T=MAX(T,T0)
      ELSE
        T0=MIN(C(2)+B(3),A(NMAX)+B(NMAX))
        T=MAX(T,T0)
        T0=MIN(A(4)+B(4),C(1)+A(3))
        T=MAX(T,T0)
        DO N=5,NMAX
          T0=MIN(A(N)+B(N),A(N-1)+B(N-1))
          T=MAX(T,T0)
        ENDDO
      ENDIF
C
      NDIMW=NINT(T-DEPS2)
C
      ALLOCATE (W6J(-1:NDIMW+1))
C
C INITIALIZE
C
      W1=DONE
      JMIN1=-1
      JMAX1=999999999
C
C INITIAL PAIR
C
      CALL WIG6JR(W6J,B(1),A(2),B(2),C(2),B(3),JMIN,JMAX,IH1,NDIMW)
C
      IF(JMIN.LT.0)GO TO 500
c
c should be caught already
      if(ih.ne.ih1)then
        write(*,*)'6j mixes integer and half integer...?',ih,ih1
        stop
      endif
C
      JMIN1=MAX(JMIN1,JMIN)                               !FOR INFO ONLY
      JMAX1=MIN(JMAX1,JMAX)                               !FOR INFO ONLY
C
      DO K=KMIN,KMAX
        W1(K)=W1(K)*W6J(K)
c        write(*,*)k,w6j(k)
      ENDDO
C
      CALL WIG6JR(W6J,C(1),A(3),A(1),A(2),B(1),JMIN,JMAX,IH1,NDIMW)
C
      IF(JMIN.LT.0)GO TO 500
c
c should be caught already
      if(ih.ne.ih1)then
        write(*,*)'6j mixes integer and half integer...?',ih,ih1
        stop
      endif
C
      JMIN1=MAX(JMIN1,JMIN)                               !FOR INFO ONLY
      JMAX1=MIN(JMAX1,JMAX)                               !FOR INFO ONLY
C
      DO K=KMIN,KMAX
        W1(K)=W1(K)*W6J(K)
c        write(*,*)k,w6j(k)
      ENDDO
C
C TO ARBITRARY NMAX
C
      DO N=4,NMAX
C
        B3=B(N-1)
        IF(N.EQ.4)B3=C(1)
C
        CALL WIG6JR(W6J,B(N),A(N),C(N-1),A(N-1),B3,JMIN,JMAX,IH1,NDIMW)
C
        IF(JMIN.LT.0)GO TO 500
c
c should be caught already
        if(ih.ne.ih1)then
          write(*,*)'6j mixes integer and half integer...?',ih,ih1
          stop
        endif
C
        JMIN1=MAX(JMIN1,JMIN)                             !FOR INFO ONLY
        JMAX1=MIN(JMAX1,JMAX)                             !FOR INFO ONLY
C
        DO K=KMIN,KMAX
          W1(K)=W1(K)*W6J(K)
c          write(*,*)k,w6j(k)
        ENDDO
C
      ENDDO
C
C COMPLETE
C
      B3=B(NMAX)
      IF(NMAX.EQ.3)B3=C(1)
C
      CALL WIG6JR(W6J,B(3),C(2),C(NMAX),B3,A(NMAX),JMIN,JMAX,IH1,NDIMW)
C
      IF(JMIN.LT.0)GO TO 500
c
c should be caught already
      if(ih.ne.ih1)then
        write(*,*)'6j mixes integer and half integer...?',ih,ih1
        stop
      endif
C
      JMIN1=MAX(JMIN1,JMIN)                               !FOR INFO ONLY
      JMAX1=MIN(JMAX1,JMAX)                               !FOR INFO ONLY
C
      DO K=KMIN,KMAX
        W1(K)=W1(K)*W6J(K)
c        write(*,*)k,w6j(k)
      ENDDO
C
C DETERMINE PHASES (might be useful to do so earlier for quicker exit)
C
      SUMR=DZERO
      DO N=1,NMAX
        SUMR=SUMR+A(N)+B(N)+C(N)
      ENDDO
      IR=NINT(SUMR-DEPS2)
      IF(IR.NE.NINT(SUMR+DEPS2))THEN                       !HALF-INTEGER
        IHR=1
      ELSE                                                      !INTEGER
        IHR=0
      ENDIF
      IS1=IR+NMAX*KMIN
      IS1=IS1-KMIN
      IS1=MOD(IS1,2)
      IS1=1-2*IS1
      if(ih.eq.0)then
        if(ihr.ne.0)then
          write(*,*)'nj sum mixes integer and half integer...?',ih,ihr
          stop
        endif
      else                  !1st & 2nd can't both exist for a given nmax
        IS=IHR+NMAX*IH
        IS=IS-IH
        if(2*(is/2).ne.is)then
          write(*,*)'nj sum mixes integer and half integer...?',ih,ihr
          stop
        else
          IS1=IS1+IS/2
          IS1=MOD(IS1,2)
          IS1=1-2*IS1
        endif
      endif
C
C PERFORM SUMMATION
C
      ISWAP=-1
      IHP=IH+1
      C1=DZERO
      wmax=0
      IF(MOD(NMAX,2).EQ.0)THEN
        DO K=KMIN,KMAX
          ISWAP=-ISWAP
          T=(2*K+IHP)*W1(K)
          C1=C1+ISWAP*T
          wmax=max(wmax,abs(t))
c          if(debug1)write(*,*)k,t,'wmax=',wmax
        ENDDO
      ELSE
        DO K=KMIN,KMAX
          T=(2*K+IHP)*W1(K)
          C1=C1+T
          wmax=max(wmax,abs(t))
c          if(debug1)write(*,*)k,t,'wmax=',wmax
        ENDDO
      ENDIF
C
      C1=C1*IS1                                        !APPLY KMIN PHASE
C
c      write(*,*)'wmax=',wmax
      IF(ABS(C1).LT.-DEPS1*wmax)C1=DZERO                        !ZEROIZE
C
      W3NJ=C1
C
  100 DEALLOCATE (W6J)
      IF(NMAX.GT.2)DEALLOCATE (W1)
C
      RETURN
C
  500 KMIN=JMIN
      KMAX=JMAX
      GO TO 100
C
      END SUBROUTINE WIG3NJ
C
C***********************************************************************
C
      SUBROUTINE WIG6JR(W6J,A2,A3,B1,B2,B3,JMIN,JMAX,IH,NDIMW)
C
C-----------------------------------------------------------------------
C
C N.R.BADNELL
C
C CALCULATES WIGNER 6-J SYMBOL ( A1 A2 A3 )
C                              ( B1 B2 B3 )
C FOR ALL ALLOWED A1 GIVEN BY JMIN+IH/2,...,JMAX+IH/2, WHERE IH=0 OR 1,
C VIA FORWARD AND BACKWARD RECURSION. OUTPUT IN W6J(J): J=JMIN,...,JMAX.
C VALUES LESS THAN DEPS1 *WITHIN THE CLASSICAL REGION* ARE SET TO ZERO.
C SEE E.G. A. R. EDMONDS "ANGULAR MOMENTUM IN QUANTUM MECHANICS" (1957).
C
C COMBINES THE ALGORITHMS (BUT NOT CODING) OF
C SCHULTEN K. AND GORDON R. G., 1975, J.MATH.PHYS., 16, 1961 AND 1971
C   AND
C LUSCOMBE J. H. AND LUBAN M., 1998, PHYS.REV.E, 57, 7274
C - SEE 
C BADNELL N. R., GUZMAN F., BRODIE S., WILLIAMS R. J. R., VAN HOOF
C P. A. M., CHATZIKOS M., & FERLAND G. J., 2021 MNRAS, 507, 2922
C
C INPUT:
C  A2,A3,B1,B2,B3: THE FIVE ARGUMENTS OF THE 6-J SYMBOL AS DEFINED ABOVE
C                  USING THEIR ACTUAL VALUE, I.E. *NOT* TWICE  (REAL*8).
C  NDIMW: THE UPPER DIMENSION OF THE ARRAY W6J(-1:NDIMW+1) (INTEGER)
C         REQUIRED IS AT LEAST JMAX.
C
C OUTPUT:
C  JMIN,JMAX,IH: THE RANGE OF ARGUMENT OF THE 6-J SYMBOL (INTEGER)
C  W6J(JMIN+IH/2,JMAX+IH/2): 6-J SYMBOL FOR ALL POSSIBLE A1. (REAL*8)
C
C  JMIN.LT.0 FLAGS AN ERROR:
C         -1 NON-INTEGER TRIANGLE SUM.
C         -2 TRIANGLE RULE FAILURE.
C         -3 DIMENSION ERROR, JMAX THEN CONTAINS THE REQUIRED NDIMW.
C
C UPDATE LOG:
C 23/12/15 - PHASE WAS NOT DETERMINED IF W6J(JMAX)=0
C 11/12/15 - REMOVED TEMP ARRAY.
C 20/11/15 - INITIAL RELEASE.
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      logical debug0,debug1,debug2
C
      PARAMETER (DZERO=0.0D0)
      PARAMETER (DHALF=0.5D0)
      PARAMETER (DONE=1.0D0)
      PARAMETER (DEPS1=1.D-16)                              !REAL*8 ZERO
      PARAMETER (DEPS2=1.D-5)
C
      PARAMETER (NDIMM=5)                       !MAX NO. MATCHING POINTS
C
      DIMENSION W6J(-1:NDIMW+1),W1(NDIMM),W2(NDIMM)
C
      data debug0/.false./,debug1/.true./,debug2/.true./
      data iwritd/0/     !=0 to screen; set < 0 for non-interactive use
C
      E(A1)=SQRT(
     X           (A1*A1-(A2-A3)*(A2-A3))*((A2+A3+1)*(A2+A3+1)-A1*A1)*
     X           (A1*A1-(B2-B3)*(B2-B3))*((B2+B3+1)*(B2+B3+1)-A1*A1)
     X           )
c     X           (A1-A2+A3)*(A1+A2-A3)*(A2+A3+1-A1)*(A2+A3+1+A1)*
c     X           (A1-B2+B3)*(A1+B2-B3)*(B2+B3+1-A1)*(B2+B3+1+A1)
c     X           )
C
      F(A1)=(A1+A1+1)*(
     X           A1*(A1+1)*(-A1*(A1+1)+A2*(A2+1)+A3*(A3+1)-2*B1*(B1+1))+
     X           B2*(B2+1)*(A1*(A1+1)+A2*(A2+1)-A3*(A3+1))+B3*(B3+1)*
     X          (A1*(A1+1)-A2*(A2+1)+A3*(A3+1))
     X                 )
c     X           A1*(A1+1)*(-A1*(A1+1)+A2*(A2+1)+A3*(A3+1)+
c     X           B2*(B2+1)+B3*(B3+1)-2*B1*(B1+1))+
c     X          (A2*(A2+1)-A3*(A3+1))*(B2*(B2+1)-B3*(B3+1))
c     X                 )
C
      X(A1)=A1*E(A1+1)
      Z(A1)=(A1+1)*E(A1)
      Y(A1)=F(A1)
C
      if(iwritd.lt.0)debug1=.false.
C
C CHECK 6J-SELECTION RULES (JMIN FLAGS FAILURE TYPE)
C
      W6J=0                                               !INITIALZE ALL
      IH=0
      JMAX=-1
C
      IF(
     X   NINT(A2+B1+B3-DEPS2).NE.NINT(A2+B1+B3+DEPS2)
     X                       .OR.
     X   NINT(A3+B1+B2-DEPS2).NE.NINT(A3+B1+B2+DEPS2)
     X                                               )THEN
        if(debug2)write(*,*)'*** sr.wig6jr: non-integer triangle sum...'
        JMIN=-1
        RETURN
      ENDIF
      IF(
     X   B3.GT.A2+B1+DEPS2.OR.B3.LT.ABS(A2-B1)-DEPS2
     X                    .OR.
     X   B2.GT.A3+B1+DEPS2.OR.B2.LT.ABS(A3-B1)-DEPS2
     X                                              )THEN
        if(debug2)write(*,*)'*** sr.wig6jr: triangle rule failure...'
        JMIN=-2
        RETURN
      ENDIF
C
C QUANTUM LIMITS
C
      J23M=NINT(ABS(A2-A3)-DEPS2)
      L23M=NINT(ABS(B2-B3)-DEPS2)
      JMIN=MAX(J23M,L23M)
      J23P=NINT(A2+A3-DEPS2)
      L23P=NINT(B2+B3-DEPS2)
      JMAX=MIN(J23P,L23P)
      IF(J23P.NE.NINT(A2+A3+DEPS2))THEN
        TH=DHALF
      ELSE
        TH=0
      ENDIF
      IH=NINT(2*TH)
      if(debug1)
     x  write(iwritd,"(/' quantum range:  ',2f9.1)")jmin+th,jmax+th
C
C DIMENSION CHECK (IN CASE USER HAS CALLED WITH HARD-WIRED DIMENSION)
C
      IF(JMAX.GT.NDIMW)THEN
        JMIN=-3
        if(debug2)write(*,*)'*** SR.WIG6JR: INCREASE NDIMW TO ',jmax
        RETURN
      ENDIF
C
C QUICK RETURN
C
      IF(JMIN.EQ.JMAX)THEN                                !GET FROM NORM
        W6J(JMIN)=1
        JMID1=JMIN
        JMD1M=JMID1-1
        JMID2=JMID1
        JMD2P=JMID2+1
        GO TO 800
      ENDIF
C
C DETERMINE NUMBER OF NODES (WOULD COUNT TO DETERMINE APPROACH TO
C CLASSICALLY FORBIDDEN REGION IF WE DID NOT HAVE THIS LIMIT.)
C
C      T=MIN(A2+B3,A3+B2)
C      NODES=NINT(T-B1)
C
C CLASSICAL LIMITS (TETRAHEDRON V=0)
C
      C2=A2+DHALF
      C2=C2*C2
      C3=A3+DHALF
      C3=C3*C3
      D1=B1+DHALF
      D1=D1*D1
      D2=B2+DHALF
      D2=D2*D2
      D3=B3+DHALF
      D3=D3*D3
      A=-D1
      B=D1*(D2+C3-D1)+C2*(D1+D2-C3)+D3*(D1-D2+C3)
      C=D3*C2*(C3+D2-D1)+D2*C3*(C2+D3-D1)-
     X  D3*C3*(C3+D3-D1)-D2*C2*(C2+D2-D1)
      B=B/(A+A)
      C=C/A
      D=SQRT(B*B-C)
      C1MIN=MAX(DZERO,-B-D)
      C1MIN=SQRT(C1MIN)-DHALF
      C1MAX=MAX(DZERO,-B+D)
      C1MAX=SQRT(C1MAX)-DHALF
      if(debug1)
     x  write(iwritd,"(' classical range:',2f9.1)")c1min,c1max
C
C MATCHINGS
C
      JMID1=INT(C1MIN+DEPS2)        !INT fallback -> JMID1=0 when JMIN=0
      JMID1=MAX(JMID1,JMIN)
      JMID1=MIN(JMID1,JMAX-1)
      JMD1M=JMID1-1
      JMID2=NINT(C1MAX+DEPS2)
      JMID2=MIN(JMID2,JMAX)
      JMID2=MAX(JMID2,JMIN+1)
      JMD2P=JMID2+1
c
      if(debug2.and.jmid1.gt.jmid2)then                !shouldn't happen
        write(mw6,*)'*** sr.wig6jr: jmid1,2=',jmid1,jmid2,jmin,jmax
        stop '*** sr.wig6jr: jmid1 .gt. jmid2'
      endif
C
      JMID=(JMID1+JMID2)/2
ct      jmid=jmid2-1                  !approx schulten & gordon matching
      JMID=MAX(JMID,JMIN)
      JMID=MIN(JMID,JMAX)
      NMATCH=1
C                                    .GT.2 COVERED BY ABOVE
      IF(JMID.GT.JMIN.AND.JMID.LT.JMAX
     X               .AND.JMAX-JMIN+1.GT.2
     X                                )NMATCH=3
C                                     BUT NEED IF NMATCH.GT.3
      NMID=NMATCH/2-1
      JMID0=JMID-NMID-2
      IF(NMATCH.GT.NDIMM)THEN                  !ONLY IF USER INCREASES 3
        JMIN=-4
        JMAX=NMATCH
        if(debug2)write(*,*)'*** SR.WIG6JR: INCREASE NDIMM TO ',nmatch
        RETURN
      ENDIF
C
C BEGIN MAIN RECURSIONS
C
C FORWARD
      IF(JMIN+IH.EQ.0)THEN                             !CASE X(0)=0=Y(0)
        if(jmid1.ne.jmin)stop '*** sr.wig6jr: jmin=0, jmid1>0...'
        JMID1=1
        JMD1M=JMID1-1
C
        E1=4*SQRT(A2*(A2+1)*B2*(B2+1))                 !E1=X(J)/J AT J=0
        F00=2*(A2*(A2+1)+B2*(B2+1)-B1*(B1+1))         !F00=Y(J)/J AT J=0
C
        W6J(JMIN)=1
        W6J(JMIN+1)=-W6J(JMIN)*F00/E1                !3-TERM all the way
        go to 100
C
c2        W6J(JMIN)=-E1/F00                             !2-TERM one step
c
c or could go backward all the way to zero, since c1min=0 here
cb        write(0,*)'going backwards'
cb        jmid=jmin
cb        nmatch=1
cb        nmid=-1
cb        jmid0=jmin-1
cb        w1(1)=1
cb        go to 200
      ELSEIF(JMID1.GT.JMIN)THEN
        T=TH+JMIN
        W6J(JMIN)=-X(T)/Y(T)
      ENDIF
C
C 2-TERM
C      W6J(JMIN-1)=0                                 !NOT USED CURRENTLY
      DO J=JMIN+1,JMD1M
        T=TH+J
        W6J(J)=-X(T)/(Y(T)+Z(T)*W6J(J-1))
      ENDDO
C
      W6J(JMID1)=1
      DO J=JMD1M,JMIN,-1
        W6J(J)=W6J(J+1)*W6J(J)
      ENDDO
C
C 3-TERM
 100  continue
      DO J=JMID1,JMID+NMID
        T=TH+J
        W6J(J+1)=-(Y(T)*W6J(J)+Z(T)*W6J(J-1))/X(T)
      ENDDO
C
      J=JMID0
      DO N=1,NMATCH
         J=J+1
         W1(N)=W6J(J)
      ENDDO
C
C BACKWARD
C
C 2-TERM
cb 200  continue
      W6J(JMAX+1)=0
      DO J=JMAX,JMD2P,-1
        T=TH+J
        W6J(J)=-Z(T)/(Y(T)+X(T)*W6J(J+1))
      ENDDO
C
      jsign=1
      W6J(JMID2)=1
      DO J=JMD2P,JMAX
        jsign=jsign*nint(sign(done,w6j(j)))!case w6j(jmax)=0 (underflow)
        W6J(J)=W6J(J-1)*W6J(J)
      ENDDO
C
C 3-TERM
      DO J=JMID2,JMID-NMID,-1
        T=TH+J
        W6J(J-1)=-(Y(T)*W6J(J)+X(T)*W6J(J+1))/Z(T)
      ENDDO
C
      J=JMID0
      DO N=1,NMATCH
         J=J+1
         W2(N)=W6J(J)
      ENDDO
C
C RELATIVE NORM
C
      T12=0
      T11=0
      DO N=1,NMATCH
        T12=T12+W1(N)*W2(N)
        T11=T11+W1(N)*W1(N)
      ENDDO
      WMATCH=T12/T11
      if(debug0)
     x  write(iwritd,"(' jmatch=',i6,'  wmatch=',f7.1)")jmid,wmatch
      DO J=JMIN,JMID0
        W6J(J)=W6J(J)*WMATCH
      ENDDO
C
C PHASE
C
 800  ISIGN=NINT(A2+A3+B2+B3)
      ISIGN=MOD(ISIGN,2)
      ISIGN=-2*ISIGN+1
C
      T=ISIGN*W6J(JMAX)
      IF(T.GT.DZERO)THEN
        ISIGN=1
      ELSEIF(T.LT.DZERO)THEN
        ISIGN=-1
      ELSE
cp        if(debug2)write(*,*)'*** sr.wig6jr: unable to determine phase'
cp        jmin=-5
cp        return
        isign=isign*jsign
      ENDIF
C
C ABSOLUTE NORM
C
      IHP=IH+1
      SUM=0
      DO J=JMIN,JMAX
        SUM=SUM+(J+J+IHP)*W6J(J)*W6J(J)
      ENDDO
      SUM=SUM*(B1+B1+1)
C
      SUM=DONE/SQRT(SUM)
      SUM=SUM*ISIGN
      if(debug0)
     x  write(iwritd,"(' wnorm=',1pe10.2)")sum
c
c some test code
c      so=sign(done,w6j(jmin))
c      nod=0
c      wmax=0
c      do j=jmin,jmax
c        if(so*w6j(j).lt.dzero)then
c          so=-so
c          nod=nod+1
c        endif
c        wmax=max(wmax,abs(w6j(j)))
c      enddo
c      if(nod.ne.nodes)write(*,*)'nodes expected/found=',nodes,nod
c      write(*,*)'wmax=',wmax                               !order unity
C
      DO J=JMIN,JMD1M
        W6J(J)=W6J(J)*SUM
      ENDDO
      DO J=JMID1,JMID2
        T=W6J(J)*SUM
        IF(ABS(T).LT.DEPS1)T=DZERO                        !*wmax ZEROIZE
        W6J(J)=T
      ENDDO
      DO J=JMD2P,JMAX
        W6J(J)=W6J(J)*SUM
      ENDDO
C
      RETURN
C
      END SUBROUTINE WIG6JR
C
C***********************************************************************
C
      SUBROUTINE WIG3JRJ(W3J,A2,A3,B2,B3,JMIN,JMAX,IH,NDIMW)
C
C-----------------------------------------------------------------------
C
C N. R. BADNELL                                                 02/12/21
C
C CALCULATES WIGNER 3-J SYMBOL ( A1 A2 A3 )
C                              ( B1 B2 B3 )
C FOR ALL ALLOWED A1 GIVEN BY JMIN+IH/2,...,JMAX+IH/2, WHERE IH=0 OR 1,
C VIA FORWARD AND BACKWARD RECURSION. OUTPUT IN W3J(J): J=JMIN,...,JMAX.
C VALUES LESS THAN DEPS1 *WITHIN THE CLASSICAL REGION* ARE SET TO ZERO.
C SEE E.G. A. R. EDMONDS "ANGULAR MOMENTUM IN QUANTUM MECHANICS" (1957).
C
C COMBINES THE ALGORITHMS (BUT NOT CODING) OF
C SCHULTEN K. AND GORDON R. G., 1975, J.MATH.PHYS., 16, 1961 AND 1971
C   AND
C LUSCOMBE J. H. AND LUBAN M., 1998, PHYS.REV.E, 57, 7274
C - SEE 
C BADNELL N. R., GUZMAN F., BRODIE S., WILLIAMS R. J. R., VAN HOOF
C P. A. M., CHATZIKOS M., & FERLAND G. J., 2021 MNRAS, 507, 2922
C
C INPUT:
C  A2,A3,B2,B3: THE FOUR ARGUMENTS OF THE 3-J SYMBOL AS DEFINED ABOVE
C               USING THEIR ACTUAL VALUE, I.E. *NOT* TWICE  (REAL*8).
C  NDIMW: THE UPPER DIMENSION OF THE ARRAY W3J(-1:NDIMW+1) (INTEGER)
C         REQUIRED IS AT LEAST JMAX.
C
C OUTPUT:
C  JMIN,JMAX,IH: THE RANGE OF ARGUMENT OF THE 3-J SYMBOL (INTEGER)
C  W3J(JMIN+IH/2,JMAX+IH/2): 3-J SYMBOL FOR ALL POSSIBLE A1. (REAL*8).
C
C  JMIN.LT.0 FLAGS AN ERROR:
C         -1 NON-INTEGER TRIANGLE SUM.
C         -2 TRIANGLE RULE FAILURE.
C         -3 DIMENSION ERROR, JMAX THEN CONTAINS THE REQUIRED NDIMW.
C
C UPDATE LOG:
C 02/12/21 - OVERALL PHASE NOT DETERMINED WHEN A2-A3-B1(=B2+B3).LT.0
C 19/02/21 - INITIAL RELEASE
C
C BELOW IS THE WIG6JR LOG, ON WHICH WIG3JRJ WAS BASED
C 23/12/15 - PHASE WAS NOT DETERMINED IF W6J(JMAX)=0
C 11/12/15 - REMOVED TEMP ARRAY.
C 20/11/15 - INITIAL RELEASE.
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      logical debug0,debug1,debug2
C
      PARAMETER (DZERO=0.0D0)
      PARAMETER (DHALF=0.5D0)
      PARAMETER (DONE=1.0D0)
      PARAMETER (DTWO=2.0D0)
      PARAMETER (DEPS1=1.D-16)                              !REAL*8 ZERO
      PARAMETER (DEPS2=1.D-5)
C
      PARAMETER (NDIMM=5)                       !MAX NO. MATCHING POINTS
C
      DIMENSION W3J(-1:NDIMW+1),W1(NDIMM),W2(NDIMM)
C
      data debug0/.false./,debug1/.true./,debug2/.true./
      data iwritd/0/     !=0 to screen; set < 0 for non-interactive use
C
      A(A1)=SQRT(
     X           (A1*A1-(A2-A3)*(A2-A3))*((A2+A3+1)*(A2+A3+1)-A1*A1)*
     X           (A1*A1-(B2+B3)*(B2+B3))
     X           )
C
      B(A1)=(A1+A1+1)*(
     X           (B2+B3)*(A2*(A2+1)-A3*(A3+1))-
     X           (B2-B3)*A1*(A1+1)
     X                 )
C
      X(A1)=A1*A(A1+1)
      Z(A1)=(A1+1)*A(A1)
      Y(A1)=B(A1)
C
      if(iwritd.lt.0)debug1=.false.
C
C CHECK 3J-SELECTION RULES (JMIN FLAGS FAILURE TYPE)
C
      W3J=0                                               !INITIALZE ALL
      IH=0
      JMAX=-1
      B1=-B2-B3
C
      IF(
     X   NINT(A2+ABS(B2)-DEPS2).NE.NINT(A2+ABS(B2)+DEPS2)
     X                       .OR.
     X   NINT(A3+ABS(B3)-DEPS2).NE.NINT(A3+ABS(B3)+DEPS2)
     X                                               )THEN
       if(debug2)write(*,*)'*** sr.wig3jrj: non-integer triangle sum...'
        JMIN=-1
        RETURN
      ENDIF
      IF(
     X   ABS(B2).GT.A2+DEPS2.OR.ABS(B3).GT.A3+DEPS2
     X                                              )THEN
        if(debug2)write(*,*)'*** sr.wig3jrj: triangle rule failure...'
        JMIN=-2
        RETURN
      ENDIF
C
C QUANTUM LIMITS
C
      J23M=NINT(ABS(A2-A3)-DEPS2)
      L23M=NINT(ABS(B1)-DEPS2)
      JMIN=MAX(J23M,L23M)
      JMAX=NINT(A2+A3-DEPS2)
      IF(J23M.NE.NINT(ABS(A2-A3)+DEPS2))THEN
        TH=DHALF
      ELSE
        TH=0
      ENDIF
      IH=NINT(2*TH)
      if(debug1)
     x  write(iwritd,"(/' quantum range:  ',2f9.1)")jmin+th,jmax+th
C
C DIMENSION CHECK (IN CASE USER HAS CALLED WITH HARD-WIRED DIMENSION)
C
      IF(JMAX.GT.NDIMW)THEN
        JMIN=-3
        if(debug2)write(*,*)'*** SR.WIG3JRJ: INCREASE NDIMW TO ',jmax
        RETURN
      ENDIF
C
C QUICK RETURN
C
      IF(JMIN.EQ.JMAX)THEN                                !GET FROM NORM
        W3J(JMIN)=1
        JMID1=JMIN
        JMD1M=JMID1-1
        JMID2=JMID1
        JMD2P=JMID2+1
        GO TO 800
      ENDIF
C
C DETERMINE NUMBER OF NODES (WOULD COUNT TO DETERMINE APPROACH TO
C CLASSICALLY FORBIDDEN REGION IF WE DID NOT HAVE THIS LIMIT.)
C
C      T=MIN(A2-B2,A3+B3)
C      NODES=NINT(T)
C
C CLASSICAL LIMITS (TRIANGLE A=0)
C
      C2=A2+DHALF
      C2=C2*C2
      C3=A3+DHALF
      C3=C3*C3
      C23=C2+C3
      D2=B2*B2
      D3=B3*B3
      D23=B2*B3
      C=SQRT((C2-D2)*(C3-D3))
      C1MIN=C23-DTWO*(C-D23)
      C1MIN=SQRT(C1MIN)-DHALF
      C1MAX=C23+DTWO*(C+D23)
      C1MAX=SQRT(C1MAX)-DHALF
      if(debug1)
     x  write(iwritd,"(' classical range:',2f9.1)")c1min,c1max
C
C MATCHINGS
C
      JMID1=INT(C1MIN+DEPS2)        !INT fallback -> JMID1=0 when JMIN=0
      JMID1=MAX(JMID1,JMIN)
      JMID1=MIN(JMID1,JMAX-1)
      JMD1M=JMID1-1
      JMID2=NINT(C1MAX+DEPS2)
      JMID2=MIN(JMID2,JMAX)
      JMID2=MAX(JMID2,JMIN+1)
      JMD2P=JMID2+1
c
      if(debug2.and.jmid1.gt.jmid2)then                !shouldn't happen
        write(*,*)'*** sr.wig3jrj: jmid1,2=',jmid1,jmid2,jmin,jmax
        stop '*** sr.wig3jrj: jmid1 .gt. jmid2'
      endif
C
      JMID=(JMID1+JMID2)/2
ct      jmid=jmid2-1                  !approx schulten & gordon matching
      JMID=MAX(JMID,JMIN)
      JMID=MIN(JMID,JMAX)
      NMATCH=1
C                                    .GT.2 COVERED BY ABOVE
      IF(JMID.GT.JMIN.AND.JMID.LT.JMAX
     X               .AND.JMAX-JMIN+1.GT.2
     X                                )NMATCH=3
C                                     BUT NEED IF NMATCH.GT.3
      NMID=NMATCH/2-1
      JMID0=JMID-NMID-2
      IF(NMATCH.GT.NDIMM)THEN                  !ONLY IF USER INCREASES 3
        JMIN=-4
        JMAX=NMATCH
        if(debug2)write(*,*)'*** SR.WIG3JRJ: INCREASE NDIMM TO ',nmatch
        RETURN
      ENDIF
C
C BEGIN MAIN RECURSIONS
C
C FORWARD
C
      IF(JMIN+IH.EQ.0)THEN                             !CASE X(0)=0=Y(0)
        if(jmid1.ne.jmin)stop '*** sr.wig3jrj: jmin=0, jmid1>0...'
        JMID1=1
        JMD1M=JMID1-1
C
        A1=2*SQRT(A2*(A2+1))                           !A1=X(J)/J AT J=0
        B00=-2*B2                                     !B00=Y(J)/J AT J=0
C
        W3J(JMIN)=1
        W3J(JMIN+1)=-W3J(JMIN)*B00/A1                !3-TERM all the way
        go to 100
C
c2        W3J(JMIN)=-A1/B00                             !2-TERM one step
c
c or could go backward all the way to zero, since c1min=0 here
cb        write(0,*)'going backwards'
cb        jmid=jmin
cb        nmatch=1
cb        nmid=-1
cb        jmid0=jmin-1
cb        w1(1)=1
cb        go to 200
      ELSEIF(JMID1.GT.JMIN)THEN
        T=TH+JMIN
        W3J(JMIN)=-X(T)/Y(T)
      ENDIF
C
C 2-TERM
C      W3J(JMIN-1)=0                                 !NOT USED CURRENTLY
      DO J=JMIN+1,JMD1M
        T=TH+J
        W3J(J)=-X(T)/(Y(T)+Z(T)*W3J(J-1))
      ENDDO
C
      W3J(JMID1)=1
      DO J=JMD1M,JMIN,-1
        W3J(J)=W3J(J+1)*W3J(J)
      ENDDO
C
C 3-TERM
 100  continue
      DO J=JMID1,JMID+NMID
        T=TH+J
        W3J(J+1)=-(Y(T)*W3J(J)+Z(T)*W3J(J-1))/X(T)
      ENDDO
C
      J=JMID0
      DO N=1,NMATCH
         J=J+1
         W1(N)=W3J(J)
      ENDDO
C
C BACKWARD
C
C 2-TERM
cb 200  continue
      W3J(JMAX+1)=0
      DO J=JMAX,JMD2P,-1
        T=TH+J
        W3J(J)=-Z(T)/(Y(T)+X(T)*W3J(J+1))
      ENDDO
C
      jsign=1
      W3J(JMID2)=1
      DO J=JMD2P,JMAX
        jsign=jsign*nint(sign(done,w3j(j)))!case w3j(jmax)=0 (underflow)
        W3J(J)=W3J(J-1)*W3J(J)
      ENDDO
C
C 3-TERM
      DO J=JMID2,JMID-NMID,-1
        T=TH+J
        W3J(J-1)=-(Y(T)*W3J(J)+X(T)*W3J(J+1))/Z(T)
      ENDDO
C
      J=JMID0
      DO N=1,NMATCH
         J=J+1
         W2(N)=W3J(J)
      ENDDO
C
C RELATIVE NORM
C
      T12=0
      T11=0
      DO N=1,NMATCH
        T12=T12+W1(N)*W2(N)
        T11=T11+W1(N)*W1(N)
      ENDDO
      WMATCH=T12/T11
      if(debug0)
     x  write(iwritd,"(' jmatch=',i6,'  wmatch=',f7.1)")jmid,wmatch
      DO J=JMIN,JMID0
        W3J(J)=W3J(J)*WMATCH
      ENDDO
C
C PHASE
C
 800  ISIGN=NINT(A2-A3+B2+B3)
      ISIGN=MOD(ISIGN,2)
      ISIGN=-2*ABS(ISIGN)+1
C
      T=ISIGN*W3J(JMAX)
      IF(T.GT.DZERO)THEN
        ISIGN=1
      ELSEIF(T.LT.DZERO)THEN
        ISIGN=-1
      ELSE
        if(debug2)write(*,*)'*** sr.wig3jrj: unable to determine phase'
cp        jmin=-5
cp        return
        isign=isign*jsign
      ENDIF
C
C ABSOLUTE NORM
C
      IHP=IH+1
      SUM=0
      DO J=JMIN,JMAX
        SUM=SUM+(J+J+IHP)*W3J(J)*W3J(J)
      ENDDO
C
      SUM=DONE/SQRT(SUM)
      SUM=SUM*ISIGN
      if(debug0)
     x  write(iwritd,"(' wnorm=',1pe10.2)")sum
c
c some test code
c      so=sign(done,w3j(jmin))
c      nod=0
c      wmax=0
c      do j=jmin,jmax
c        if(so*w3j(j).lt.dzero)then
c          so=-so
c          nod=nod+1
c        endif
c        wmax=max(wmax,abs(w3j(j)))
c      enddo
c      if(nod.ne.nodes)write(*,*)'nodes expected/found=',nodes,nod
c      write(*,*)'wmax=',wmax                               !order unity
C
      DO J=JMIN,JMD1M
        W3J(J)=W3J(J)*SUM
      ENDDO
      DO J=JMID1,JMID2
        T=W3J(J)*SUM
        IF(ABS(T).LT.DEPS1)T=DZERO                        !*wmax ZEROIZE
        W3J(J)=T
      ENDDO
      DO J=JMD2P,JMAX
        W3J(J)=W3J(J)*SUM
      ENDDO
C
      RETURN
C
      END SUBROUTINE WIG3JRJ
C
C***********************************************************************
C
      SUBROUTINE YLV2WIG(NMAX,A,B,C)
C
C RE-MAPS YUTSIS ET AL ORDERED 3N-J ARGUMENTS TO WIGNER-ORDERED.
C SEE E.G. A. R. EDMONDS "ANGULAR MOMENTUM IN QUANTUM MECHANICS" (1957),
C ORD-SMITH PHYS.REV.94 1227 (1954) AND YUTSIS, LEVINSON & VANAGAS
C "MATHEMATICAL APPARATUS OF THE THEORY OF ANGULAR
C MOMENTUM" ISRAEL PROGRAM FOR SCIENTIFIC TRANSLATIONS: JERUSALEM (1962)
C
C INPUT A(N), B(N), C(N) THE YUTSIS ET AL 3N-J ARGUMENTS FOR N=1,...NMAX
C
C OUTPUT: A(N), B(N), C(N) ARE OVERWRITTEN BY THE WIGNER ORDERING,
C FOR NMAX.GT.2.
C
C FOR NMAX=2, A(3) & B(3) ARE UNDEFINED ON INPUT AND C(1) & C(2) ARE
C UNDEFINED ON OUTPUT. C(3) IS NEVER REFERENCED THEN.
C BEAR IN MIND WHEN USING THAT THE YLV62 PHASE IS THAT OF
C THE RACAH W-COEFFICIENT, NOT WIGNER 6-J.
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      DIMENSION A(*),B(*),C(*)                     !* FOR CASE OF NMAX=2
C
      ALLOCATABLE :: D(:),E(:),F(:)          !TEMP HOLD OF YLV ARGUMENTS
C
C QUICK RETURN
C
      IF(NMAX.LE.1)THEN
        WRITE(*,*)'*** NO REMAPPING OF 3N-J ARGUMENTS EXIST FOR N.LE.1'
        RETURN
      ENDIF
C
C HANDLE TRIVIAL CASE
C
      IF(NMAX.EQ.2)THEN
        A(3)=B(1)
        B(3)=B(2)
        B(1)=C(1)
        B(2)=C(2)
        RETURN
      ENDIF
C
      ALLOCATE (D(NMAX),E(NMAX),F(NMAX))
C
      D(NMAX)=A(1)
      E(NMAX)=C(1)
      F(NMAX)=B(NMAX)
C
      I=NMAX
      DO N=2,NMAX-3
        I=I-1
        D(I)=A(N)
        E(I)=C(N)
        F(I)=B(N-1)
      ENDDO
C
      IF(NMAX.GT.3)THEN
        I=I-1
        D(I)=A(NMAX-2)
        F(I)=B(NMAX-3)
      ENDIF
      E(I)=C(NMAX)
      I=I-1
      D(I)=A(NMAX-1)
      E(I)=B(NMAX-1)
      F(I)=A(NMAX)
      I=I-1
      D(I)=B(NMAX-2)
      E(I)=C(NMAX-1)
      F(I)=C(NMAX-2)
C
      DO N=1,NMAX
        A(N)=D(N)
        B(N)=E(N)
        C(N)=F(N)
      ENDDO
C
      DEALLOCATE (D,E,F)
C
      RETURN
C
      END SUBROUTINE YLV2WIG
C
C***********************************************************************