C N. R. BADNELL                                                 02/12/21
C
      PROGRAM TWIG3JRM
C
C-----------------------------------------------------------------------
C TEST
C EVALUATION OF WIGNER 3-J SYMBOL ARRAY BY RECURRENCE IN M (SR.WIG3JRM)
C
C EXAMPLES:
C
C SCHULTEN & GORDON  (CPC 11, 269 (1976))
C 240  120  140   -20
C  16   15   13     2
C  32   30   26     4
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
C
      PARAMETER (IWRITE=77)             !=0 TO SCREEN, >0 TO WIG3JRM.out
C
      PARAMETER (DEPS2=1.D-5)
C
      ALLOCATABLE :: W3J(:)
C
C
      IF(IWRITE.GT.0)OPEN(IWRITE,FILE='WIG3JRM.out')   !STATUS='UNKNOWN'
C
C ASK USER TO ENTER THE *FOUR* PREDEFINED ARGUMENTS a,b,c;d  (as f=-d-e)
C NO PROBLEMS IF THEY ENTER FIVE OR SIX HERE AS JUST IGNORED
C
      WRITE(*,*)
     X"DETERMINES WIGNER 3-J SYMBOL W3J(e) FOR ALL e FOR GIVEN a,b,c;d"
C
   1  WRITE(*,*)
      WRITE(*,*)"ENTER THE 4 INTEGERS (a,b,c;d) AT *TWICE* "
     X         ,"ACTUAL VALUE (all.lt.0 exits):"
C
      READ(*,*)IA,IB,IC,ID
C
      IF(IA.LT.0.OR.IB.LT.0.OR.IC.LT.0)THEN
        WRITE(*,*)'*** EXITING AS IA,IB,IC NOT ALL .GE. ZERO'
        IF(IWRITE.GT.0)WRITE(*,*)'*** W3J(e) OUTPUT IN WIG3JRM.out'
        STOP
      ENDIF
c
ct       call cpu_time(timei)
C
      A1=IA
      A1=A1/2
      A2=IB
      A2=A2/2
      A3=IC
      A3=A3/2
      B1=ID
      B1=B1/2
C
      T=MAX(-A2,-A3-B1)
      NL=NINT(T-DEPS2)
      T=MIN(A2,A3-B1)
      NU=NINT(T-DEPS2)
      NDIMW=MAX(-NL,NU)
C
      ALLOCATE (W3J(-1-NDIMW:NDIMW+1))
C
ct      do i=1,10000
      CALL WIG3JRM(W3J,A1,A2,A3,B1,JMIN,JMAX,IH,NDIMW)
ct      enddo
c
ct       call cpu_time(timef)
ct       times=timef-timei
ct       write(*,"(' recursion time=',1pe9.1,' msecs')")times*1000
C
      IF(JMAX.LT.JMIN)THEN
        IF(JMIN.EQ.-3)THEN                            !SHOULD NOT HAPPEN
          WRITE(*,*)'*** DIMENSION FAILURE, INCREASE NDIMW TO ',-JMAX
        ELSE
          WRITE(*,*)'*** SELECTION RULE FAILURE, ALL 3-J EQUAL ZERO!'
          IF(JMIN.EQ.-1)WRITE(*,*)'*** non-integer triangle sum'
          IF(JMIN.EQ.-2)WRITE(*,*)'*** triangle rule failure'
        ENDIF
        DEALLOCATE (W3J)
        GO TO 1
      ENDIF
C
      IF(IWRITE.GE.0)
     X   WRITE(IWRITE,"(/4X,'2E',5X,'W3J(A,B,C; D,E,F)',4X
     X         ,'2A 2B 2C; 2D 2F =',3I6,1X,I6,2X,'-2D-2E')")IA,IB,IC,ID
C
      DO JW=JMIN,JMAX
        IE=JW+JW+IH
        IF(IWRITE.GE.0)WRITE(IWRITE,"(I6,1P2D25.16)")IE,W3J(JW)
      ENDDO
      IF(IWRITE.GE.0)WRITE(IWRITE,*)
c
      if(iwrite.gt.0)call flush(IWRITE)
C
      DEALLOCATE (W3J)
C
      GO TO 1
C
      END PROGRAM TWIG3JRM
C
C***********************************************************************
C
      SUBROUTINE WIG3JRM(W3J,A1,A2,A3,B1,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 B2 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  A1,A2,A3,B1: 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 DIMENSION OF THE ARRAY W3J(-1-NDIMW:NDIMW+1) (INTEGER).
C         REQUIRED IS AT LEAST MAX(JMAX,|JMIN|).
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 B2. (REAL*8).
C
C  JMIN.GT.JMAX FLAGS AN ERROR, THEN JMIN=
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 WIG3JRM 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 (DFOUR=4.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: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
      C(B2)=SQRT(
     X           (A2-B2+1)*(A2+B2)*(A3-B2-B1+1)*(A3+B2+B1)
     X           )
C
      D(B2)=A2*(A2+1)+A3*(A3+1)-A1*(A1+1)-(B2+B2)*(B2+B1)
C
      X(B2)=C(B2+1)
      Z(B2)=C(B2)
      Y(B2)=D(B2)
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=-4
C
      IF(
     X   NINT(A1+A2+A3-DEPS2).NE.NINT(A1+A2+A3+DEPS2)
     X                       .OR.
     X   NINT(A1+ABS(B1)-DEPS2).NE.NINT(A1+ABS(B1)+DEPS2)
     X                                               )THEN
       if(debug2)write(*,*)'*** sr.wig3jrm: non-integer triangle sum...'
        JMIN=-1
        RETURN
      ENDIF
      IF(
     X   ABS(B1).GT.A1+DEPS2
     X                    .OR.
     X   A1.GT.A2+A3+DEPS2.OR.A2.GT.A1+A3+DEPS2.OR.A3.GT.A1+A2+DEPS2
     X                                              )THEN
        if(debug2)write(*,*)'*** sr.wig3jrm: triangle rule failure...'
        JMIN=-2
        RETURN
      ENDIF
C
C QUANTUM LIMITS
C
      J31M=NINT(-A3-B1-DEPS2)
      J31P=NINT(A3-B1-DEPS2)
      J2M=NINT(-A2-DEPS2)
      J2P=NINT(A2-DEPS2)
      JMIN=MAX(J2M,J31M)
      JMAX=MIN(J2P,J31P)
      IF(NINT(A2-DEPS2).NE.NINT(A2+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
        JMAX=JMAX+1
        if(debug2)write(*,*)'*** SR.WIG3JRM: INCREASE NDIMW TO ',jmax
        jmax=-jmax
        RETURN
      ENDIF
      IF(JMIN.LT.-NDIMW)THEN
        JMAX=-JMIN
        JMIN=-3
        if(debug2)write(*,*)'*** SR.WIG3JRM: INCREASE NDIMW TO ',jmax
        jmax=-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
      NODES=NINT(A2+A3-A1)
C
C CLASSICAL LIMITS (TRIANGLE A=0)
C
      C1=A1+DHALF
      C1=C1*C1
      C2=A2+DHALF
      C2=C2*C2
      C3=A3+DHALF
      C3=C3*C3
      C123=C1+C2-C3
      D1=B1*B1
      E=(C1-D1)*(DFOUR*C1*C2-C123*C123)
      E=SQRT(E)
      E1=-B1*C123
      E2=C1+C1
      C1MIN=(E1-E)/E2
      C1MAX=(E1+E)/E2
      if(debug1)
     x  write(iwritd,"(' classical range:',2f9.1)")c1min,c1max
C
C MATCHINGS
C                                   !INT fallback -> JMID1=0 when JMIN=0
      JMID1=INT(C1MIN+sign(DEPS2,c1min))
      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.wig3jrm: jmid1,2=',jmid1,jmid2,jmin,jmax
        stop '*** sr.wig3jrm: 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.WIG3JRM: INCREASE NDIMM TO ',nmatch
        RETURN
      ENDIF
C
C BEGIN MAIN RECURSIONS
C
C FORWARD
C
      IF(ABS(JMIN)+IH.EQ.0)THEN
C
        Y0=Y(DZERO)
c
c can have Y(0)=0 but X(0).ne.0 (as W3J(1)=0) so use 3-term all the eway
C
        IF(ABS(Y0).LT.DEPS1)THEN
c
c          if(jmid1.ne.jmin)then               !for info
c            write(mw6,*)'*** sr.wig3jrm: jmin=0, jmid1,2=',jmid1,jmid2
cc            write(mw6,*)' a1,a2,a3,b1=',A1,A2,A3,B1
c          endif
C
          JMID1=1
          JMD1M=JMID1-1
C
          W3J(JMIN)=1
          W3J(JMIN+1)=-W3J(JMIN)*Y0/X(DZERO)       !CAN TEST FOR Y0.NE.0
          go to 100
c
c or could go backward all the way to zero, if 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
c
C       ELSE 2-TERM IS FINE STILL & SAFER FOR A TURNING POINT .GT. JMIN
C
        ENDIF
C
      ENDIF
C
      IF(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-B1)
      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
cp        if(debug2)write(*,*)'***sr.wig3jrm: unable to determine phase'
cp        jmin=-5
cp        return
        isign=isign*jsign
      ENDIF
C
C ABSOLUTE NORM
C
      SUM=0
      DO J=JMIN,JMAX
        SUM=SUM+W3J(J)*W3J(J)
      ENDDO
      SUM=SUM*(A1+A1+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,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 WIG3JRM
C
C***********************************************************************