C N. R. BADNELL                   NRB  v6.0.1                   11/05/21
C
      PROGRAM TLCHNG
C
C-----------------------------------------------------------------------
C
C TEST-DRIVER/WRAPPER FOR L-CHANGING COLLISIONS (DELTA-N=0~DELTA-E).
C
C EVALUATES QM RESULTS OF VRINCEANU ET AL AP.J.747:56 (2012) AND,
C FOR DIPOLE TRANSITIONS, COMPARES WITH MODIFIED PENGELLY & SEATON
C RESULTS (PSM20) OF BADNELL ET AL MNRAS XXX, XXX (2021) WHICH 
C IMPROVES-UPON THE ORIGINAL MODIFIED PENGELLY SEATON (PSM17) METHOD
C INTRODUCED BY GUZMAN ET AL MNRAS 464, 312 (2017).
C
C INTERACTIVE VERSION ASKS USER FOR TEMPERATURE, DENSITY, N-RANGE,
C L-RANGE AND MULTIPOLE (=DELTA-L).
C
C HARDWIRED IS THE SET-UP OF PROJECTILE & TARGET MASSES & CHARGES.
C DEFAULT IS PROTONS ON HYDROGEN - OTHER EXAMPLES ARE GIVEN BELOW.
C
C OUTPUT IS MAXWELLIAN RATE COEFFICIENTS QM AND PSM20 (DIPOLE ONLY).
C THE PSM20 RESULTS ARE LABELLED PS AND PSX. THE FIRST (PS) IS A
C NUMERICAL MIRROR OF THE QM CALCULATION. PSX IS THE ANALYTIC RESULT.
C THUS (PSX-PS)/PSX PROVIDES A GUIDE TO THE QUADRATURE ERROR OF THE QM.
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
      LOGICAL bout,bout7,bsum,brslv
C
      PARAMETER (DZERO=0.0D0)
      PARAMETER (DONE=1.0D0)
      PARAMETER (DTEN=10.0D0)
      PARAMETER (D1P10=1.D10)
      PARAMETER (D1M20=1.D-20)
C
      PARAMETER (IZERO=0)
      PARAMETER (IONE=1)
C
      PARAMETER (TKAY=1.5789D5)         !KELVIN TO RYDBERGS
      PARAMETER (BOHRA=5.29177D-9)      !CM TO A.U.
      PARAMETER (EMASS=5.485798D-4)     !ELECTRON MASS A.M.U
C                                       !SAMPLE PROJECTILES
      PARAMETER (PROTON=  1.00727647)
      PARAMETER (DEUTERON=2.01355323)
      PARAMETER (ALPHAP=  4.00150618)
C                                       !SAMPLE TARGETS
      PARAMETER (HYDROGEN= 1.007825)
      PARAMETER (DEUTERIUM=2.014102)
      PARAMETER (TRITIUM=  3.016049)
      PARAMETER (HELIUM3=3.016029)
      PARAMETER (HELIUM4=4.002603)
      PARAMETER (HELIUM=HELIUM4)
      PARAMETER (CARBON=12.0)
C
      PARAMETER (T2NEBIG=1.D+20)        !MAX ALLOWED VALUE OF T**2/Ne
      PARAMETER (TOLT2NE=1.D-03)        !MIN ALLOWED VALUE OF T**2/Ne
c
      data ntest0/400/,ac0/5.e-4/       !for largest "safe" n for T/N
C
      PI=ACOS(-DONE)                    !F2003 PARAMETER
C
C***********************************************************************
C
      write(*,*)
      write(*,*)
     &       '*** L-CHANGING COLLISION RATE COEFFICIENTS (v6.0.1) ***'
      write(*,*)
      write(*,*)
     &       '*** MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***'
      write(*,*)
C
C I/O SET-UP
C    (SET UNIT NUMBER IOUT NEGATIVE TO SUPPRESS WRITES THERE)
C
      iout= 0                                  !=0 rate coeffs to screen
      bout=iout.ge.0
      if(iout.gt.0)open(iout,file='tlchng.out')           !=6 to outfile
c
      if(iout.gt.0)write(iout,
     &"(/'# *** L-CHANGING COLLISION RATE COEFFICIENTS (v6.0) ***'/)")
      if(iout.gt.0)write(iout,
     &"('# *** MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***'/)"
     &)
c
      iout7= 7                                  !.lt. to suppress writes
      bout7=iout7.gt.0
      if(bout7)open(iout7,file='tlchng.sum')!summd/unresolvd rate coeffs
c
      if(bout7)write(iout7,
     &"(/'# *** L-CHANGING COLLISION RATE COEFFICIENTS (v6.0) ***'/)")
      if(bout7)write(iout7,
     &"('# *** MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***'/)
     &")
C
C***********************************************************************
C************************** BEGIN USER SET-UP **************************
C***********************************************************************
C
C SET-UP PROJECTILE AND TARGET CHARGES (ZP,ZT) AND MASSES (PMASS,TMASS):
C
      ZP=1                              !PROJECTILE CHARGE
      PMASS=PROTON                      !PROJECTILE MASS
C
C FOR IONIC TARGETS, ATTEMPT TO ACCOUNT FOR MASS OF MISSING ELECTRONS
C
      NZ0=1                             !TARGET NUCLEAR CHARGE
      NELC=1                            !NO. OF TARGET ELCTRNS (INC RYD)
C
      ZT=NZ0-NELC+1                     !TARGET CHARGE (AS SEEN BY RYD)
C
      TMASS=HYDROGEN                    !TARGET MASS (INC. RYD)
C
      TMASS=TMASS-(NZ0-NELC)*EMASS      !APPROX FOR IONIC TARGETS
C
C***********************************************************************
C*************************** END USER SET-UP ***************************
C***********************************************************************
C
C-----------------------------------------------------------------------
c
      write(*,"(/' TARGET MASS=',f14.6,' AMU'
     &        ,5x,'TARGET CHARGE(RYD)=',f4.1)")tmass,zt
      write(*,"(/' PROJECTILE MASS=',f10.6,' AMU'
     &        ,5x,'PROJECTILE CHARGE=',f5.1/)")pmass,zp
      if(iout.gt.0)write(iout,"(/'# TARGET MASS=',f14.6,' AMU'
     &                         ,5x,'TARGET CHARGE(RYD)=',f4.1)")tmass,zt
      if(iout.gt.0)write(iout,"(/'# PROJECTILE MASS=', f10.6,' AMU'
     &                         ,5x,'PROJECTILE CHARGE=',f5.1/)")pmass,zp
      if(bout7)write(iout7,"(/'# TARGET MASS=',f14.6,' AMU'
     &                      ,5x,'TARGET CHARGE(RYD)=',f4.1)")tmass,zt
      if(bout7)write(iout7,"(/'# PROJECTILE MASS=',f10.6,' AMU'
     &                      ,5x,'PROJECTILE CHARGE=',f5.1/)")pmass,zp
C
C-----------------------------------------------------------------------
C
C READ TEMPERATURE IN KELVIN AND DENSITY IN /CM^3 FOR ANY CUT-OFF
C (DENE.LT.0 CAN FLAG A RANGE OF DENSITIES, SEE BELOW.)
C
   1  WRITE(0,*)
      WRITE(0,*)"Enter T(K) and N(/cm^3):    (Negative to EXIT)"
C
      READ(*,*)TEE,DENE
C
      IF(TEE.LE.DZERO.OR.DENE.EQ.DZERO)STOP 'WE ARE DONE'
C
      IF(DENE.LT.DZERO)THEN
C
        WRITE(0,*)"Enter No. densities and Log increment (canbe .lt.0):"
C
        READ(*,*)NDEN,DENINC
C
        IF(NDEN.LE.0)STOP 'NEED NDEN.GT.0 FOR DENE.LT.0 FLAG'
        IF(DENINC.eq.DZERO)STOP 'NEED DENINC.NE.0 FOR DENE.LT.0 FLAG'
C
        DENE=-DENE
        if(deninc.gt.0)then
          denmx=log10(dene)+(nden-1)*deninc
          denmx=dten**denmx
        else
          denmx=dene
        endif
      ELSE
        NDEN=1
        DENINC=DZERO
        denmx=dene
      ENDIF
c
      t2ne=tee*tee/denmx
c
      if(t2ne.lt.tolt2ne)then
        write(*,"(/' *** Temperature/Density ratio too small',
     &  ' - set T**2/Ne =',1p,e9.2,' .gt.',e9.2)")t2ne,tolt2ne
        go to 1
      endif
c
      if(t2ne.gt.t2nebig)then
        write(*,"(/' *** SR.LCHNG: Temperature/Density ratio too large'
     &  ,' - set T**2/Ne =',1p,e9.2,' .lt.',e9.2)")t2ne,t2nebig
        go to 1
      endif
c
      write(*,*)
      write(*,"(/' T=',1p,e9.2,' K ',' N=',e9.2,' /cm**3')")
     &                   tee,abs(dene)
      debye=tee/(tkay*8*pi*abs(dene)*bohra**3)   !debye radius in a.u.
      debye=sqrt(debye)
      write(*,"(/' R_Debye= ',6x,1p,e8.2,' au',3x,1p,e8.2,' cm'/)")
     &                   debye,debye*bohra
C
C USER SPECIFIES RYDBERG L-RANGE AND N-RANGE TO LOOP-OVER AND DELTA-L
C (THE RANGES CAN BE IN EITHER DIRECTION, DESPITE MIN,MAX LABELS)
C
      NINC=1
      LINC=1
C
  10  WRITE(0,*)
      WRITE(0,*)
     &  "ENTER lmin, lmax, delta-l; nmin, nmax:   (5*0 to EXIT)"
C
      READ(*,*)L1MIN,L1MAX,LD0,NMIN,NMAX      !,LINC,NINC
c
      write(*,*)
C
      IF(NMIN+NMAX.EQ.0)GO TO 1
C
C DIRECTION OF L-LOOP IS GOVERNED BY L1MIN AND L1MAX, SINCE LINC IS
C NOT NORMALLY READ, BUT IN CASE READ, CHECK AND RESET IF NECESSARY.
C
      LINC=ABS(LINC)
      LINC=MAX(IONE,LINC)
C
C RESTRICT MULTIPOLES FOR PRODUCTION CODE
C
      LD=ABS(LD0)
      LD=MAX(LD,IONE)
      IF(LD.GT.3)THEN                          !FOR PRODUCTION CODE
        WRITE(0,*)'*** SET MULTIPOLE (LD) .LE.3 FOR PRODUCTION CODE'
        GO TO 10
      ENDIF
C
C-----------------------------------------------------------------------
c
c set some (outout) flags
c
      bsum=bout.and.abs(linc).eq.ld    !remove bout?
      bout7=bout7.and.bsum
      brslv=.false.
C
C SET-UP FOR UNRESOLVED PS CASE.
C
      LDX=LD0
      LDM=LD0
C
      IF(LD0.EQ.0)THEN           !SYNC RESOLVED QM WITH UNRESOLVED PS
C
        write(*,*)
        write(*,*)
     &  '*** MODIFIED PENGELLY & SEATON USES UNRESOLVED DIPOLE ***'
        if(iout.gt.0)write(iout,
     &"(/'# *** MODIFIED PENGELLY & SEATON USES UNRESOLVED DIPOLE ***')"
     & )
        if(bout7)write(iout7,
     &"(/'# *** MODIFIED PENGELLY & SEATON USES UNRESOLVED DIPOLE ***'/)
     &")
C
        IF(L1MIN.GT.L1MAX)THEN
C
C          L=L1MIN               !REVERSE SINCE SR.LCHNG DOES NOT KNOW
C          L1MIN=L1MAX           !TO "CHANGE" DIRECTION FOR QM AS LD0=0,
C          L1MAX=L               !UNLESS FLAG L1->-L1-1 (USER DANGEROUS)
c
c comment-out above and uncomment "crev" below to "test"
c backwards recursion - a little less accurate than forwards
c also need to uncomment "crev" flag l10=-l10-1 above call sr.lchng
c  *and* in sr.lchng uncomment l1=-l1-1
c
crev
          l1min=l1min+1
          ldm=-ld
crev
        ENDIF
C
        IF(L1MIN.LE.L1MAX)THEN
          L1MIN=L1MIN-1
          LDM=LD
        ENDIF
c
c see if we can generate resolved ps from unresolved ps by recurrence
c
        brslv=(ldm.gt.0.and.l1min.le.0
     &     .or.ldm.lt.0.and.l1min+1.ge.max(nmin,nmax))
        if(brslv)then
          ldx=ldm
          write(*,*)
          write(*,*)
     &      '*** RESOLVED MODIFIED PENGELLY & SEATON FROM RECURSION ***'
          write(*,*)
          if(iout.gt.0)write(iout,"(/
     &    '# *** RESOLVED MODIFIED PENGELLY & SEATON FROM RECURSION ***'
     &    /)")
        endif
C
      ENDIF
C
C ALLOW L TO RANGE IN EITHER DIRECTION (AND CHECK/ADJUST MIN VALUE)
C
      LI=L1MIN
      LF=L1MAX
      IF(LI.LE.LF)THEN
        LI=MAX(-LDX,LI,IZERO)
        LF0=LF
      ELSE
        LF=MAX(-LDX,LF,IZERO)
        LI0=LI
        LINC=-LINC
      ENDIF
C
C CHECK MAX N FOR PRODUCTION CODE
C
      NMN=MIN(NMIN,NMAX)
      IF(NMN.GT.1000)THEN                    !FOR PRODUCTION CODE
        WRITE(0,*)'*** NMIN.GT.1000 TOO LARGE FOR PRODUCTION CODE'
        GO TO 10
      ENDIF
      NMN=MAX(IZERO,NMN)
C
      NMX=MAX(NMIN,NMAX)
      IF(NMX.GT.1000)THEN                   !FOR PRODUCTION CODE
        WRITE(0,*)'*** REDUCING NMAX TO 1000 FOR PRODUCTION CODE'
        NMX=1000
      ENDIF
C
      NINC=MAX(IONE,NINC)
C
C ALLOW N TO RANGE IN EITHER DIRECTION
C
      IF(NMIN.LE.NMAX)THEN
        NI=NMN
        NF=NMX
      ELSE
        NI=NMX
        NF=NMN
        NINC=-NINC
      ENDIF
c
c try warn and flag largest recommended n for the smallest t2ne, where
c ps agrees with psx to within ~5%, i.e. we flag numerical failure.
c this assumes we know ac, the accuracy parameter which defines the mesh
c this is for pure debye cut-off only.
c tbd: accommodate lifetime/delta-e cut-off - more severe.
c
      ac=5.d-4        !***must sync with sr.lchng usage to be meaningful
      t2=sqrt(t2ne)
      ntest=5*ntest0*sqrt(t2)*(ac0/ac)**0.175d0
      if(nmx.gt.ntest)then
        write(*,
     &    "(/'(N.B. Debye only) Largest recommended N for this T/Ne is '
     &     ,i4/)")ntest
      endif
C
C-----------------------------------------------------------------------
C
C *** START DENSITY LOOP
C
C-----------------------------------------------------------------------
C
      call cpu_time(timei)
C
      DEN0=LOG10(DENE)
C
      DO NE=1,NDEN
C
        TE=DEN0+(NE-1)*DENINC
        DENE=DTEN**TE
c
c output headers
c
        if(bout)then
          debye=tee/(tkay*8*pi*dene*bohra**3)      !debye radius in a.u.
          debye=sqrt(debye)
          if(iout.gt.0.or.ne.gt.1)then
            write(iout,20010)tee,dene
            write(iout,20020)debye,debye*bohra
          endif
          if(iout7.gt.0)then
            write(iout7,20030)tee,dene
            write(iout7,20040)debye,debye*bohra
          endif
          if(ld.eq.1)then
            if(iout.eq.0)then
              write(iout,20050)
              write(iout,20060)
            else
              write(iout,20070)
              write(iout,20080)
            endif
            if(iout7.gt.0)then                                   !summed
              write(iout7,20090)
              write(iout7,20100)
            endif
          else                !quadrupole+, resolved onky - maybe summed
            if(iout.eq.0)then
              write(iout,20110)
            else
              write(iout,20120)
            endif
            if(iout7.gt.0)then
              write(iout7,20130)
            endif
          endif
          nold=-1
        endif
C
C-----------------------------------------------------------------------
C
C *** START N-LOOP
C
C-----------------------------------------------------------------------
C
        DO N=NI,NF,NINC
C
          N1=N-1
C
C-----------------------------------------------------------------------
C
C *** START L-LOOP
C
C-----------------------------------------------------------------------
C
          IF(LINC.GT.0)THEN
            LF=MIN(LF0,N1-LDX,N1)
          ELSE
            LI=MIN(LI0,N1-LDX,N1)
          ENDIF
C
          DO L1=LI,LF,LINC
C
            TLIFE=D1P10             !USER CAN SUPPLY RYD LIFETIME IN SEC
C
            DELTAE=D1M20            !USER CAN SUPPLY RYD DELTA-E IN RYDS
C
            IERR=0                  !INIT AS .LT. 0 JUST DE-ALLOCATES
C
            L10=L1
crev
            if(ld0.eq.0.and.ldm.lt.0)l10=-l10-1!flag backwards recursion
C
C
            CALL LCHNG(PMASS,TMASS,ZP,ZT,N,L10,LD0,TEE,DENE,TLIFE,DELTAE
     &                ,QQM,QPS,QPSX,IERR)
C
C
            IF(IERR.GT.0)THEN
              WRITE(0,*)'IERR=',IERR
              STOP '*** FAILURE IN SR.LCHNG'
            ENDIF
C
C-----------------------------------------------------------------------
c
c tbd: write info on R_delta-e/life(v_rms) if exists e.g. .lt. 2*R_Debye
c
c write-out QQM,QPS,QPSX rate coefficients (cm^3/s) in illustrative ways
c *** the user will want to manipulate and store them according to their
c own needs rather than slavishly following what is in sr.output.***
c
            if(bout)then
c
              call output(iout,iout7,bsum,brslv,nold,pmass,tmass,n,l1
     &             ,ld0,ld,ldm,linc,tee,debye,tlife,deltae,qqm,qps,qpsx)
c
              nold=n
            endif
C
C-----------------------------------------------------------------------
C
          ENDDO                         !END L-LOOP
C
C-----------------------------------------------------------------------
C
        ENDDO                           !END N-LOOP
C
C-----------------------------------------------------------------------
C
      ENDDO                             !END DENSITY LOOP
C
C-----------------------------------------------------------------------
C
C OPTIONAL TIDY-UP: DE-ALLOCATE ARRAYS USED & RETURN TO INITIAL STATE
C
      IERR=-1
C
      CALL LCHNG(PMASS,TMASS,ZP,ZT,N,L1,LD0,TEE,DENE,TLIFE,DELTAE
     &          ,QQM,QPS,QPSX,IERR)
C
      IF(IERR.GT.0)THEN
        WRITE(0,*)'IERR=',IERR
        STOP '*** FAILURE IN SR.LCHNG'
      ENDIF
C
C-----------------------------------------------------------------------
c
      if(iout.gt.0)call flush(iout)     !case user exits ungracefully
      if(iout7.gt.0)call flush(iout7)   !case user exits ungracefully
c
      call cpu_time(timef)
      times=timef-timei
      write(*,"(/1pe9.2,'  secs')")times
C
      GO TO 10                          !LETS DO IT ALL AGAIN
C
C-----------------------------------------------------------------------
C
20010 format(/'# T=',1p,e9.2,' K ',' N=',e9.2,' /cm**3')
20020 format(/'# R_Debye= ',1p,e9.2,' au',3x,1p,e8.2,' cm'/)
20030 format(/'# T=',1p,e9.2,' K ',' N=',e9.2,' /cm**3')
20040 format(/'# R_Debye= ',1p,e9.2,' au',3x,1p,e8.2,' cm'/)
20050 format(/20x,'RATE COEFFS (CM**3/S)')
20060 format(3x,'N',2x,'L1->L2',3x,'L',4x
     &       ,'QM',7x,'PS',7x,'PSX',4X,'(QM-PS)/QM',1X,'(PSX-PS)/PSX')
20070 format(/'#',19x,'RATE COEFFS (CM**3/S)')
20080 format('#',2x,'N',2x,'L1->L2',3x,'L',4x
     &       ,'QM',7x,'PS',7x,'PSX',4X,'(QM-PS)/QM',1X,'(PSX-PS)/PSX')
20090 format(/'#',15x,'RATE COEFFS (CM**3/S)')
20100 format('#',2x,'N',3x,'L',8x
     &       ,'QM',7x,'PS',7x,'PSX',4X,'(QM-PS)/QM',1X,'(PSX-PS)/PSX')
20110 format(/3x,'N',2x,'L1->L2',3x,'L',4x,'QM',4X
     &                  ,': RATE COEFFS (CM**3/S)')
20120 format(/'#',2x,'N',2x,'L1->L2',3x,'L',4x,'QM',4X
     &                  ,': RATE COEFFS (CM**3/S)')
20130 format(/'#',2x,'N',3x,'L',8x,'QM',4X
     &                  ,': RATE COEFFS (CM**3/S)')
C-----------------------------------------------------------------------
C
      END PROGRAM TLCHNG
C
C***********************************************************************
C
      subroutine output(iout,iout7,bsum,brslv,nold,pmass,tmass,n,l1
     &        ,ld0,ld,ldm,linc,tee,debye,tlife,deltae,qqm,qps,qpsx)
c
c-----------------------------------------------------------------------
c
c sr.output illustrates example output based-on numerical qqm and qps
c rate coefficients and the analytic qpsx.
c the remaining input is described/defined by the calling routine.
c
c e.g.
c
c 1/ the orginal qqm,qps,qpsx is written, it may be final-state resolved
c or, in the case of ps, unresolved by final-state.
c
c 2/ if resolved then l->l' can be summed-over final l'=l+/-1.
c
c 3/ if unresolved then resolved p&s data can be determined by recursion
c from l=0 or l=n-1, if the input l-data *starts* there for a given n.
c
c note, starting at l=0 and recurring-up is more accurate
c overall than starting at l=n-1 and recurring down.
c
c these last two options are for dipole data only.
c
c In addition, fractional differences are constructed between numerical
c qps and qqm as well as analystic qpsx and numerical qps.
c warnings can be flagged if the differences become too large, or even
c the run aborted.
c
c-----------------------------------------------------------------------
c
      implicit real*8(a-h,o-z)
      implicit integer*4 (i-n)
c
      logical bsum,brslv
c
      parameter (dzero=0.0d0)
      parameter (d1p10=1.d10)
      parameter (d1m20=1.d-20)
c
      parameter (tkay=1.5789d5)         !kelvin to rydbergs
      parameter (bohra=5.29177d-9)      !cm to a.u.
      parameter (bohrt=2.41889d-17)     !secs
      parameter (emass=5.485798d-4)     !electron mass a.m.u
c
      parameter (tolp =2.d-2)           !warn when (psx-ps)/psx .gt.
      parameter (tolpp=2.d-1)           !abort when (psx-ps)/psx .gt.
      parameter (tolpq=1.d10)           !warn when (qm-ps)/qm .gt.
c
      parameter (xl=0.72)               !p&s (30)
      parameter (xd=1.12)               !p&s (29)
      parameter (idtest=2)
c
c-----------------------------------------------------------------------
c
      rmass=pmass*tmass/(pmass+tmass)   !reduced mass
      rmass=rmass/emass
      vrms=sqrt(tee/(rmass*tkay))       !v_rms
c
      dlife=min(tlife,d1p10)
      dlife=xl*dlife*vrms/bohrt         !tlife in sec
c
      ddele=max(deltae,d1m20)
      ddele=2*xd*vrms/ddele             !deltae in ryd
c
      dtest=min(dlife,ddele)
c
      if(dtest.lt.idtest*debye)then     !write for info, at least
        if(dtest.lt.dlife)then
          write(iout,"(/1x,'R_delta-E(RMS)='
     &   ,1p,e9.2,' au',3x,1p,e8.2,' cm'/)")ddele,ddele*bohra
        else
          write(iout,"(/1x,'R_life(RMS)   ='
     &   ,1p,e9.2,' au',3x,1p,e8.2,' cm'/)")dlife,dlife*bohra
        endif
      endif
c
c initialze for summing-over, if resolved
c
      if(bsum.and.n.ne.nold)then
        qqm0=0
        if(brslv.or.ld0.ne.0)then
          qps0=0
          qpsx0=0
        endif
      endif
c
c basic resolved
c
      if(ld0.eq.0)then
c ps is unresolved but qm is still to be summed, so skip until later
      elseif(ld.eq.1)then                                  !all resolved
        diffq=dzero
        if(qqm.ne.dzero)diffq=(qqm-qps)/qqm
        diffp=dzero
        if(qpsx.ne.dzero)diffp=(qpsx-qps)/qpsx
        write(iout,"(3i4,5x,1p,3e9.2,2(2x,e9.2))")n,l1,l1+ld0
     &                                    ,qqm,qps,qpsx,diffq,diffp
c flag numerics
        if(abs(diffp).gt.tolp)
     &  write(iout,
     &  "('*** LIKELY NUMERICAL ERROR IN QM: (PSX-PS)/PSX ='
     &  ,1P,E9.2)")diffp
        if(abs(diffp).gt.tolpp)then
          write(iout,
     &    "('*** NUMERICAL BREAKDOWN: (PSX-PS)/PSX ='
     &    ,1P,E9.2)")diffp
          if(iout.gt.0)call flush(iout)
          STOP '*** NUMERICAL FAILURE, INCREASE T**2/Ne'
        endif
        if(abs(diffq).gt.tolpq)then
          write(iout,
     &    "('*** PSM INACCURACY: (QM-PS)/QM ='
     &    ,1P,E9.2)")diffq
        endif
      else
        write(iout,"(3i4,5x,1p,e9.2)")n,l1,l1+ld0,qqm
      endif
c
c sum resolved (may only need to sum qm) and write-out totals
c
      if(bsum)then                                         !sum resolved
c
        l11=2*l1+1                                            !symetrize
        qqm=qqm*l11
        qps=qps*l11
        qpsx=qpsx*l11
c
c qm is always resolved but ps may be explicitly unresolved. if so, then
c if we start at l1=0 or n-1 we can use recurrence to obtain resolved ps
c which is consistent with the unresolved.
c
        lp=l1
        if(ldm.lt.0)lp=lp-ld                                 !min(l1,l2)
        lw=lp
        if(linc.lt.0)lw=lw+ld
        lw1=2*lw+1
c
        if(qqm0.ne.dzero)then
          qqm_sum=(qqm+qqm0)/lw1
          if(ld0.ne.0)then                                 !resolved sum
            qps_sum=(qps+qps0)/lw1
            qpsx_sum=(qpsx+qpsx0)/lw1
          else                                               !unresolved
            if(l11.ne.lw1)stop 'l11.ne.lw1'
            qps_sum=qps/lw1
            qpsx_sum=qpsx/lw1
            if(brslv)then                  !can apply recurrence
              qps=qps-qps0        !to get consistent resolved ps
              qpsx=qpsx-qpsx0
            endif
          endif
        else
          go to 250
        endif
c
  150   continue
c
        iout0=(iout+abs(ld0))*iout7
        if(ld.eq.1)then                                          !dipole
          if(brslv)then
            diffq=dzero
            if(qqm.ne.dzero)diffq=(qqm-qps)/qqm
            diffp=dzero
            if(qpsx.ne.dzero)diffp=(qpsx-qps)/qpsx
            write(iout,"(3i4,5x,1p,3e9.2,2(2x,e9.2))")
     &          n,l1,l1+ldm,qqm/lw1,qps/lw1,qpsx/lw1,diffq,diffp
          endif
          diffq=dzero
          if(qqm_sum.ne.dzero)diffq=(qqm_sum-qps_sum)/qqm_sum
          diffp=dzero
          if(qpsx_sum.ne.dzero)diffp=(qpsx_sum-qps_sum)/qpsx_sum
          if(iout0.eq.0)then
            write(iout0,"(i4,8x,i4,1x,1p,3e9.2,2(2x,e9.2))")n,lw
     &       ,qqm_sum,qps_sum,qpsx_sum,diffq,diffp
          endif
          if(iout7.gt.0)then
            write(iout7,"(i4,i4,4x,1x,1p,3e9.2,2(2x,e9.2))")n,lw
     &       ,qqm_sum,qps_sum,qpsx_sum,diffq,diffp
          endif
c flag numerics
          if(ld0.eq.0)then
            if(abs(diffp).gt.tolp)
     &      write(iout,
     &      "('*** LIKELY NUMERICAL ERROR IN QM: (PSX-PS)/PSX ='
     &      ,1P,E9.2)")diffp
            if(abs(diffp).gt.tolpp)then
              write(iout,
     &        "('*** NUMERICAL BREAKDOWN: (PSX-PS)/PSX ='
     &        ,1P,E9.2)")diffp
              if(iout.gt.0)call flush(iout)
              STOP '*** NUMERICAL FAILURE, INCREASE T**2/Ne'
            endif
            if(abs(diffq).gt.tolpq)then
              write(iout,
     &        "('*** PSM INACCURACY: (QM-PS)/QM ='
     &        ,1P,E9.2)")diffq
            endif
          endif
c
        else                                                !quadrupole+
          if(brslv)then
            write(iout,"(3i4,5x,1p,e9.2)")n,l1,l1+ld0,qqm
          endif
          if(iout0.eq.0)then
            write(iout0,"(i4,8x,i4,1x,1p,e9.2)")n,lw,qqm_sum
          endif
          if(iout7.gt.0)then
            write(iout7,"(i4,i4,4x,1x,1p,e9.2)")n,lw,qqm_sum
          endif
        endif
c
  250   continue
c
        if(lp+ld.eq.n-1)then                  !pick-up last/first single
          lp=lp+ld
          lw=lp
          lw0=lw1
          lw1=2*lw+1
          qqm_sum=qqm/lw1
          if(ld0.ne.0.or.ldm.lt.0.or.brslv)then
            qps_sum=qps/lw1
            qpsx_sum=qpsx/lw1
            if(brslv)lw1=lw0
            go to 150
          else                         !ld0=0 has one more l-loop for ps
            qqm0=qqm                                 !as qm.eq.zero then
          endif
        elseif(lp.eq.0)then                   !pick-up first/last single
          lp=-1
          lw=0
          lw1=1
          qqm_sum=qqm
          if(ld0.ne.0.or.ldm.gt.0.or.brslv)then
            qps_sum=qps
            qpsx_sum=qpsx
            if(brslv)lw1=3
            go to 150
          else                         !ld0=0 has one more l-loop for ps
            qqm0=qqm                                 !as qm.eq.zero then
          endif
        else
          qqm0=qqm
          if(brslv.or.ld0.ne.0)then
            qps0=qps
            qpsx0=qpsx
          endif
        endif
c
      endif
c
      return
c
      end subroutine output
C
C***********************************************************************
C
      SUBROUTINE LCHNG(PMASS,TMASS,ZP,ZT,N,L10,LD0,TEE0,DENE0
     &                ,TLIFE0,DELTAE0,QQM,QPS,QPSX,IERR)
C
C-----------------------------------------------------------------------
C
C  SR.LCHNG CALCULATES MAXWELLIAN RATE COEFFICIENTS FOR L-CHANGING
C    COLLISIONS, I.E. THOSE WHERE THE TRANSITION ENERGY DELTA-E IS
C    APPROX ZERO. E.G. N L1->N L2, WHERE L2=L1+LD0 AND LD=ABS(LD0)=1
C    FOR DIPOLE. LD=0 FLAGS UNRESOLVED (DIPOLE) PS, QM RESOLVED STILL.
C    RESULTS FOR L1->L1+LD ARE OBTAINED FROM L1+LD->L1 AND THE
C    APPLICATION OF THE RECIPROCITY RELATION, IN THIS PRODUCTION CODE.
C
C    RESULTS ARE PROVIDED FOR THE
C    QUANTUM MECHANICAL (QM) APPROXIMATION DETAILED BY
C    VRINCEANU ET AL AP.J. 747:56 (2012)
C    AND (DIPOLE ONLY) FOR THE
C    MODIFIED PENGELLY SEATON (PSM20) APPROXIMATION DETAILED BY
C    BADNELL ET AL MNRAS XXX , XXX (2021), WHICH IMPROVES-UPON THE
C    ORIGINAL MODIFIED PENGELLY SEATON (PSM17) INTRODUCED BY
C    GUZMAN ET AL MNRAS 464, 312 (2017)
C    IN THIS CASE BOTH NUMERICAL AND ANALYTIC RESULTS ARE GIVEN.
C
C  INPUT (COMPULSORY):
C
C    PMASS: IS THE PROJECTILE MASS IN AMU. E.G. =1.0072764 FOR A PROTON.
C
C    TMASS: IS THE TARGET MASS IN AMU, CONTAINING A RYDBERG ELECTRON
C          .E.G. =1.007825 FOR THE HYDROGEN ATOM.
C
C    ZP: IS THE CHARGE OF THE PROJECTILE E.G. PROTON=1, ALPHA PARTICLE=2
C
C    ZT: IS THE CHARGE OF THE TAGRET **AS SEEN BY THE RYDBERG ELECTRON**
C        SO ZT=1 FOR NEUTRAL HYDROGEN, AND ZT=2 FOR HE+.
C
C    N:  IS THE PRINCIPAL QUANTUM NUMBER OF THE RYDBERG ELECTRON,
C        RESTRICTED TO BE .LE. 1000 FOR THE PRODUCTION CODE.
C
C    L10:IS THE INITIAL ORBITAL ANGULAR MOMENTUM OF THE RYDBERG ELECTRON
C
C    LD0:DENOTES THE 2**LD-POLE ELECTRIC MULTIPOLE TRANSITION OPERATOR.
C        NORMALLY, LD=ABS(LD0)=1 FOR DIPOLE. BUT LD=2 OR 3 ARE ALSO
C        ALLOWED. LD=0 FLAGS UNRESOLVED DIPOLE PS, BUT QM SUMMED.
C
C    TEE:  IS THE ELECTRON TEMPRATURE IN KELVIN USED TO CALCULATE THE
C          MAXWELLIAN RATE COEFFICIENTS.
C
C    DENE: IS THE *ELECTRON* NUMBER DENSITY (*NOT* HYDROGEN DENSITY)
C          IN /CM^3 USED TO DETERMINE THE DEBYE CUT-OFF FOR THE
C          LARGE IMPACT PARAMETER.
C
C    IERR: .LT.0 FLAGS DE-ALLOCATION OF PREVIOUSLY ALLOCATED ARRAYS AND
C              RETURNS THE ROUTINE TO ITS ORIGINAL UNINITIALIZED STATE.
C          .GE.0 IS OTHERWISE IGNORED.
C
C  INPUT (OPTIONAL):
C
C    TLIFE:  IS THE LIFETIME, IN SECONDS, OF THE RYDBERG NL-STATE.
C
C    DELTAE: IS THE DE/EXCITATION ENERGY, IN RYD, OF THE L-CHANGING
C            TRANSITION, .LE. ZERO IS IGNORED.
C
C    TLIFE AND/OR DELTAE PROVIDE AN ALTERNATIVE EFFECTIVE CUT-OFF,
C    FOLLOWING PENGELLY & SEATON MNRAS 127, 165 (1964).
C    THE MINIMUM CUT-OFF VALUE IS USED.
C
C    N.B. THE INPUT (EXCLUDING IERR) IS NOT REDEFINED.
C
C  OUTPUT
C
C    QQM: IS THE QUANTUM MECHANICAL (QM) RATE COEFFICIENT FOR THE
C         TRANSITION N L1->L1-LD IN CM^3/S (NUMERICAL).
C
C    QPS: IS THE MODIFIED PENGELLY-SEATON (PSM20) RATE COEFFICIENT
C         FOR THE TRANSITION N L1->L1-LD IN CM^3/S (NUMERICAL).
C
C    QPSX:IS THE MODIFIED PENGELLY-SEATON (PSM20) RATE COEFFICIENT
C         FOR THE TRANSITION N L1->L1-LD IN CM^3/S (ANALYTIC).
C
C    IERR: .GT. ZERO FLAGS FAILURE, USUALLY DUE TO USER INPUT.
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8(A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
      LOGICAL BMOD,BPLTOT,BREV,BMIN,BTWO,BTWOH,BCUT,BFIRST
     &       ,bexp,bout,bout8,bout9,bex
C
C DO *NOT* CHANGE (INCREASE) MXNVAL - THE PRODUCTION CODE IS LIMITED TO
C A MAX N OF 1000.
C
      PARAMETER (MXNVAL=1000)           !MAX PRINCIPAL QUANTUM NO. ALLOW
C
C SET THE USUAL CONSTANTS
C
      PARAMETER (IZERO=0)
      PARAMETER (IONE=1)
C      PARAMETER (ITWO=2)                !CONFLICTS WITH BTWO VARIABLE!!
C
      PARAMETER (DZERO=0.0D0)
      PARAMETER (DHALF=0.5D0)
      PARAMETER (D3QRT=0.75D0)
      PARAMETER (D3HALF=1.5D0)
      PARAMETER (DONE=1.0D0)
      PARAMETER (DTWO=2.0D0)
      PARAMETER (DTHREE=3.0D0)
C
      PARAMETER (D1P10=1.D10)
      PARAMETER (D1M20=1.D-20)
C
      PARAMETER (EULER=0.57721566490153286D0)
C
C SET SOME I/O CONVERSION FACTORS
C
      PARAMETER (EMASS=5.485798D-4)     !ELECTRON MASS A.M.U
      PARAMETER (TKAY=1.5789D5)         !KELVIN TO RYDBERGS
      PARAMETER (BOHRA=5.29177D-9)      !CM TO A.U.
      PARAMETER (BOHRT=2.41889D-17)     !SECS
C
      PARAMETER (CQ=BOHRA**3/BOHRT)     !a_0^3/tau_0
C      PARAMETER (CR=2.6775D9)           !alpha^3/3*hbar
C
      PARAMETER (T2NEBIG=1.D+20)        !MAX ALLOWED VALUE OF T**2/Ne
      PARAMETER (TOLT2NE=1.D-15)        !MIN ALLOWED VALUE OF T**2/Ne
C
C FOR SR.GAMAF (FACTORIALS)
C
      PARAMETER (SCALE=1.0D6)           !SYNC WITH SR.GAMAF
      PARAMETER (MXDFS=2*MXNVAL)        !REQUIRE 2*NMAX
C
      ALLOCATABLE :: GAM(:),JGAM(:)
      ALLOCATABLE :: TQM(:),TPS(:),ENERG(:)
      ALLOCATABLE :: AL(:)
      ALLOCATABLE :: W6J(:)
C
      SAVE
C
      DATA BFIRST/.TRUE./,ierrout/0/,iun6/-1/,debyo/dzero/
      data lsum0/999/                            !0 to restrict QM L-SUM
C
C-----------------------------------------------------------------------
C
      IF(IERR.LT.0)THEN                           !ALLOW USER TO TIDY-UP
        IF(ALLOCATED(GAM))THEN
          DEALLOCATE (GAM,JGAM,STAT=JERR)
          IF(JERR.NE.0)THEN
            if(ierrout.ge.0)WRITE(ierrout,*)
     &                  '*** SR.LCHNG: DE-ALLOCATION FAILS FOR SR.GAMAF'
            IERR=22
            RETURN
          ENDIF
        ENDIF
        BFIRST=.TRUE.
        IERR=0
        RETURN
      ENDIF
C
C-----------------------------------------------------------------------
C
      IF(BFIRST)THEN
C
        BFIRST=.FALSE.
C
C-----------------------------------------------------------------------
C
C SET SOME CONSTANTS ETC
C
        PI=ACOS(-DONE)
        RTPI=SQRT(PI)
        DPHI=-PI
C
C-----------------------------------------------------------------------
C
C INITIALIZE FACTORIAL ARRAY
C
        IF(ALLOCATED(GAM))THEN
          if(ierrout.ge.0)WRITE(ierrout,*)
     &                 '*** SR.LCHNG: ERROR, SR.GAMAF ALREADY ALLOCATED'
          IERR=20
          RETURN
        ENDIF
C
        ALLOCATE(GAM(MXDFS),JGAM(MXDFS),STAT=JERR)
C
        IF(JERR.NE.0)THEN
          if(ierrout.ge.0)WRITE(ierrout,*)
     &                     '*** SR.LCHNG: ALLOCATION FAILS FOR SR.GAMAF'
          IERR=21
          RETURN
        ENDIF
C
        CALL GAMAF(GAM,JGAM,MXDFS)
C
C***********************************************************************
C************  USER JOE SHOULD NOT CHANGE ANYTHING BELOW HERE **********
C***********************************************************************
C
C-----------------------------------------------------------------------
C
C I/O SET-UP
C    (SET UNIT NUMBER NEGATIVE TO SUPPRESS WRITES THERE)
C
      iout=-4                            !=0 xsctn/rate coeffs to screen
      bout=iout.ge.0         !more details/debug than tlchng.out in main
      if(iout.eq.6)inquire(file='tlchng.out',number=iun6) !=6 to outfile
      if(iout.gt.0.and.iun6.ne.6)open(iout,file='lchng.out')
c
      iout8=-8                                  !.lt. to suppress writes
      bout8=iout8.gt.0
      if(bout8)open(iout8,file='lchng.probs')       !probabilities vs IP
c
      iout9=-9                                  !.lt. to suppress writes
      bout9=iout9.gt.0
      if(bout9)open(iout9,file='lchng.xscns')       !cross sections vs E
C
C-----------------------------------------------------------------------
C
C
C AC IS AN ACCURACY PARAMETER. TIME SCALES AS THE SQUARE OF LOG(AC).
C
C     AC=5.D-3 IS A GOOD VALUE FOR PRODUCTION, SAVE HIGN-N HIGH DENSITY.
C     AC=5.D-4 COVERS MOST CASES AND IS A TEST OF PRODUCTION PRECISION.
C     AC=5.D-5 IS FOR EXTREME CASES WHERE RATE COEFF. IS VERY SMALL
C          E.G. VERY LOW T (<100K) AND HIGHLY NON-DEGENERATE (SMALL IP).
C
      AC=5.D-4
C
C-----------------------------------------------------------------------
C
C ALLOW FOR VARIATION FROM USUAL PENGELLY & SEATON (1964):
C
C BMOD=.T. THEN PHALF=1/2 PMOD=2/3 GIVES MODIFIED; FIXED B-RATIO=1/2.
C          CAN ALSO SWITCH SMALL R BEHAVIOR FROM ORIGINAL R/R_C TO R/R_1
C
C BMOD=.F. THEN PHALF=1/2 PMOD=1 GIVES USUAL P&S; L-DEPENDENT B-RATIO.
C
      BMOD=.TRUE.                       !            ***DO NOT CHANGE***
C
      IF(BMOD)THEN                      !***NOTE CHANGE IN DEF OF PHALF!
C
        PHALF=DHALF                     !R1 DEFINED BY P=PONE=PHALF*B.R.
        IPOW=1                          !IF NON-INTEGER DEFINE AS REAL*8
        PMOD=DTWO/(2+IPOW)              !THEN P=PONE*(R/R1)**IPOW :R<R1
        imin=2                          !use PSM17 r/rcut when rcut<r1
        IMIN=1                          !USE PSM20 R/R1 WHEN RCUT<R1
C
        if(bout)then
          if(imin.eq.1)then
            write(iout,*)'*** SR.LCHNG: '
     &          ,'MODIFIED PENGELLY & SEATON (PSM20) IN USE (R/R_1) ***'
          else
            write(iout,*)'*** SR.LCHNG: '
     &          ,'MODIFIED PENGELLY & SEATON (PSM17) IN USE (R/R_C) ***'
          endif
        endif
c
      ELSE
c
c vos12 define p&s p_ij=1/2 for r<r1, then "need" phalf=done
        PHALF=DHALF                     !R1 DEF. BY P=PONE=PHALF*B-RATIO
        IPOW=0                          !ONLY USED IF BMOD=.TRUE.
        PMOD=DONE                       !THEN P=PONE FOR R<R1
        IMIN=0
        if(bout)write(iout,*)'*** SR.LCHNG: '
     &          ,'CLASSIC PENGELLY & SEATON (PS64) IN USE!!! ***'
        if(ierrout.ge.0)write(ierrout,*)'*** SR.LCHNG: '
     &          ,'CLASSIC PENGELLY & SEATON (PS64) IN USE!!! ***'
c
      ENDIF
C
      if(imin.eq.1.and.ipow.gt.1.and.bout)write(iout,*)'*** SR.LCHNG: '
     &       ,'ANALYTIC MODIFIED PENGELLY & SEATON IPOW>1 NOT CODED ***'
C
C-----------------------------------------------------------------------
C
C BMIN=.T./.F. FLAGS +/- CONTRIB FROM E=0-E_MIN FOR P&S
C
      bmin=.False.                       !STANDARD P&S NEGLECTS 0-E_EMIN
c
      BMIN=bmin.or.BMOD                                 !PSM INCLUDES IT
C
C-----------------------------------------------------------------------
c
c bexp=.f. analytic PS uses original P&S approx of exponential integral
c bexp=.t. analytic PS uses exact exponential integral E_min-infty
c          numerical is always exact, of course.
c
      bexp=.False.                            !ORIGINAL P&S APPROXES E_1
c
      BEXP=bexp.or.BMOD                             !PSM DOES NOT APPROX
C
      if(bout.and..not.bexp)write(iout,*)'*** SR.LCHNG: PS64 '
     & ,'Exponential integral approx. in use (analytic) - numerical '
     & ,'result is exact still'
C
C-----------------------------------------------------------------------
C
C BPLTOT=.F. GENERATES FINAL STATE RESOLVED DATA WITH APPROPRIATE
C            BRANCHING RATIOS AND MATCHING POINTS.
C
C BPLTOT=.T. GENERATES TOTAL UNRESOLVED BY FINAL STATE DATA SO UNIT
C            BRANCHING AND MATCHING POINT BASED ON TOTAL PROBABILTY.
C            NOTE, THIS DOES *NOT* EQUAL L-SUMMED MODIFIED P&S, WHICH
C            USES DIFFERENT MATCHING POINTS FOR L-1,L+1; NOR ORGINAL
C            P&S DETERMINED BY RECIPROCITY.
C
      BPLTOT=LD0.EQ.0
C
      if(bout.and.bpltot)write(iout,*)
     x'*** NOTE: PENGELLY & SEATON USES UNRESOLVED DIPOLE TERM ***'
C
C-----------------------------------------------------------------------
C
C BREV=.T. CALCULATE L->L+1 FROM L+1->L AND RECIPROCITY
C     =.F. CALCULATE L->L+1 EXPLICITLY *AND* WARN - NOT FOR PRODUCTION
C
      BREV=.TRUE.                       !            ***DO NOT CHANGE***
C
      IF(.NOT.BREV.AND.LD0.GT.0)THEN
        if(bout)write(iout,*)
     &    '*** SR.LCHNG: PSM WARNING, FOR L -> L+1 YOU SHOULD CALCULATE'
     &           ,' L+1 -> L AND APPLY DETAILED BALANCE ***'
      ENDIF
C
C-----------------------------------------------------------------------
C
C ENERGY DEPENDENT CUT-OFF WILL EXCEED DEBYE FOR LARGE ENOUGH ENERGY.
C
C BTWO=.TRUE.COMBINES USE OF BOTH, SWITCHING OVER AT THE MATCHING ENERGY
C            SEE GUZMAN ET AL (2017) OR BADNELL ET AL (2021) FOR DETAILS
C            THIS IS APPLIED TO ALL NUMERICAL AND ANALYTIC PSM METHODS.
C            (IF R(E=V_RMS).GT.IDTEST*R(DEBYE) IT IS RE-SET .FALSE.)
C
C BTWO=.FALSE. USES DEBYE OR ENERGY DEPENDENT AT ALL ENERGIES,
C              DEPENDING ON WHICH IS THE SMALLER AT V_RMS.
C
      btwo=.False.                      !STANDARD P&S USES ONE OR OTHER
C
      BTWO=BTWO.OR.BMOD                 !PSM (&QM) USES COMBINED 
C
      if(.not.btwo.and.bout)write(iout,*)'*** ATTENTION: CUT-OFF IS '
     &   ,'DEBYE *OR* LIFETIME/DELTA-E **ONLY**, *NOT* COMBINED ***'
C
      IF(BTWO)THEN
        IDTEST=2         !USE PURE DEBYE WHEN R(VRMS).GT.IDTEST*R(DEBYE)
      ELSE
        IDTEST=1                        !AS NO COMBINED CASE POSSIBLE
      ENDIF
C
      BTWOH=BTWO                        !HOLD TO ALLOW BTWO RE-SET .F.
C
C SET XL & XD FOR POSSIBLE LIFETIME & DELTA-E USAGE
C
      x0=d3qrt                          !use 3*V_RMS/4
      X0=DONE                           !P&S
C
      X=0.72                            !P&S (30)
      XL=X*X0
C
      X=1.12                            !P&S (29)
      XD=X*X0
C
C-----------------------------------------------------------------------
C
C BCUT=.T. APPLIES LARGE I.P. CUT-OFF TO ALL DELTA-L TRANSITIONS.
C
C     =.F. APPLIES TO DIPOLE ONLY; (EFFECTIVE) COLLISION STRENGTHS
C          FOR DELTA-L.GT.1 ARE THEN THEIR CONSTANT LIMITING VALUE.
C          "ONLY" REQUIRED FOR RCUT.LT.~0.1 CM, I.E.
C          FOR H, VERY LOW-T/HIGH DENSITY;
C          FOR HE, AND/OR HIGHLY NON-DEGENERATE.
C
C RCUT0 .GT. 0 IS ANY USER-SUPLLIED TEST CUT-OFF (IN A.U.)
C       .LT. 0 (ANY) ALLOWS REDEFINITION BY DEBYE/LIFETIME/DELTA-E
C
      BCUT=.TRUE.
C
      RCUT0=-DONE                        !ALLOWS CUT-OFF TO BE REDEFINED
C
      BCUT=RCUT0.NE.DZERO.AND.(BCUT.OR.ABS(LD0).LE.1)  !SET I.P. CUT-OFF
C
      IF(.NOT.BCUT)THEN
        IF(ABS(LD0).LE.1)THEN
          if(ierrout.ge.0)WRITE(ierrout,*)
     &  '*** SR.LCHNG: MUST FLAG AN I.P. CUT-OFF FOR DIPOLE TRANSITIONS'
          IERR=12
          RETURN
        ENDIF
        if(bout)WRITE(iout,*)
     &    '*** SR.LCHNG: ANY I.P. CUT-OFF IGNORED FOR DELTA-L=',ABS(LD0)
      ENDIF
C
C-----------------------------------------------------------------------
C
C***********************************************************************
C************  USER JOE SHOULD NOT CHANGE ANYTHING ABOVE HERE **********
C************  AND CERTAINLY SHOULD NOT CHANGE ANYTHING BELOW **********
C***********************************************************************
C
C
      ENDIF                             !END INITIALIZATION
C
C-----------------------------------------------------------------------
C
C INITIALIZE OUTPUT INCASE ABORT
C
      QQM=DZERO
      QPS=DZERO
      QPSX=DZERO
C
C-----------------------------------------------------------------------
C
C CHECK USER INPUT
C
      IERR=0
C
      IF(PMASS.LE.DZERO)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                  '*** SR.LCHNG: PROJECTILE MASS MUST BE POSITIVE'
        IERR=101
        RETURN
      ENDIF
C
      IF(NINT(PMASS).LT.IONE)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                 '*** SR.LCHNG: PROJECTILE MASS MUST BE IN A.M.U.'
        IERR=102
        RETURN
      ENDIF
C
      IF(TMASS.LE.DZERO)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                      '*** SR.LCHNG: TARGET MASS MUST BE POSITIVE'
        IERR=103
        RETURN
      ENDIF
C
      IF(NINT(TMASS).LT.IONE)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                     '*** SR.LCHNG: TARGET MASS MUST BE IN A.M.U.'
        IERR=104
        RETURN
      ENDIF
C
      IF(ZP.LE.DZERO)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                '*** SR.LCHNG: PROJECTILE CHARGE MUST BE POSITIVE'
        IERR=201
        RETURN
      ENDIF
C
      IF(ZT.LE.DZERO)THEN              !ZT IS CHARGE AS SEEN BY RYDBERG
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                '*** SR.LCHNG: TARGET CHARGE MUST BE NON-NEGATIVE'
        IERR=202
        RETURN
      ENDIF
C
      IF(N.LT.0)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
        WRITE(0,*)'*** SR.LCHNG: N MUST BE NON-NEGATIVE'
        IERR=301
        RETURN
      ENDIF
C
      IF(N.GT.MXNVAL)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &           '*** SR.LCHNG: N (>1000) TOO LARGE FOR PRODUCTION CODE'
        IERR=302
        RETURN
      ENDIF
C
      L1=L10
crev
      if(l10.lt.0)l1=-l1-1               !re-set l1 (backwards recursion
C
      IF(L1.LT.0)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                            '*** SR.LCHNG: L MUST BE NON-NEGATIVE'
        IERR=401
        RETURN
      ENDIF
C
      IF(L1.GE.N)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                             '*** SR.LCHNG: L MUST BE LESS THAN N'
        IERR=402
        RETURN
      ENDIF
C
      IF(ABS(LD0).GT.3)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &      '*** SR.LCHNG: MULTIPOLE (>3) TOO LARGE FOR PRODUCTION CODE'
        IERR=403
        RETURN
      ENDIF
C
      IF(TEE0.LE.DZERO)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &             '*** SR.LCHNG: TEMPERATURE (KELVIN) MUST BE POSITIVE'
        IERR=501
        RETURN
      ENDIF
C
      IF(DENE0.LE.DZERO)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                   '*** SR.LCHNG: DENSITY (/CM3) MUST BE POSITIVE'
        IERR=502
        RETURN
      ENDIF
C
C TESTS ON USER TLIFE0 & DELTAE0?
C
C...
C
C-----------------------------------------------------------------------
C
C INPUT O.K. SO NOW CONVERT UNITS ETC
C
      ZP2=ZP*ZP                         !PROJECTILE CHARGE SQUARED
      ZT2=ZT*ZT                         !TARGET CHARGE SQUARED
C
      RMASS=PMASS*TMASS/(PMASS+TMASS)   !REDUCED MASS
      RMASS=RMASS/EMASS
C
      TEE=TEE0
      DENE=DENE0
c
      t2ne=tee*tee/dene
c
      if(t2ne.lt.tolt2ne)then
        write(0,"(/' *** SR.LCHNG: Temperature/Density ratio too small'
     &  ,' - set T**2/Ne =',1p,e9.2,' .gt.',e9.2)")t2ne,tolt2ne
        ierr=1
        return
      endif
c
      if(t2ne.gt.t2nebig)then
        write(0,"(/' *** SR.LCHNG: Temperature/Density ratio too large'
     &  ,' - set T**2/Ne =',1p,e9.2,' .lt.',e9.2)")t2ne,t2nebig
        ierr=1
        return
      endif
C
      TEE=TEE/TKAY                      !KELVIN TO RYD
C
      DENE=DENE*BOHRA**3                !CM TO A.U.
C
      VRMS2=TEE/RMASS
      VRMS=SQRT(VRMS2)                  !V_RMS
C
      CQQ=CQ*2*RTPI/VRMS                !FOR PROBABILITY TO RATE COEFF.
C
      TCM=2*RMASS*RMASS/CQQ             !AS SEATON INCLUDES RMASS^(-3/2)
C      TCM=TCM/(AN*AN)                 !FOR UPSILON/N^4
C
C-----------------------------------------------------------------------
C
C DEBYE CUT-OFF FOR USER T/NE (PRINCIPALLY FOR DIPOLE TRANSITIONS)
C
C-----------------------------------------------------------------------
C
      DEBYE2=TEE/(8*PI*DENE)
      DEBYE=SQRT(DEBYE2)                !DEBYE RADIUS IN A.U.
c
      if(debye.ne.debyo)then
        if(bout)write(iout,"(/' R_Debye= ',6x,1p,e9.2,' au',3x,1p
     &                     ,e8.2,' cm')")debye,debye*bohra
        if(bout8)write(iout8,"(/'# R_Debye= ',6x,1p,e9.2,' au',3x,1p
     &                     ,e8.2,' cm')")debye,debye*bohra
        if(abs(ld0).le.1.and.
     &     bout9)write(iout9,"(/'# R_Debye= ',6x,1p,e9.2,' au',3x,1p
     &                     ,e8.2,' cm')")debye,debye*bohra
      endif
c
      if(bout.and.debye.ne.debyo)then
        if(abs(ld0).le.1)then
          write(iout,"(/3x,'N',3x,'L',2x,'V*N',3x,'QM'
     x                  ,7x,'PS',7x,'PSX',5X,': SIGMA/N**4  (CM**2)')")
        elseif(abs(ld0).eq.2)then
          write(iout,"(/3x,'N',3x,'L',2x,'V*N',3x,'QM'
     x                 ,24X,': SIGMA/N**4  (CM**2)')")
        else
          write(iout,"(/3x,'N',3x,'L',2x,'V*N',3x,'QM'
     x                 ,6x,'--',7x,'---',6X,': SIGMA/N**4  (CM**2)')")
        endif
        write(iout,"(44x,'RATE COEFFS (CM**3/S)')")
      endif
c
      debyo=debye
C
C-----------------------------------------------------------------------
C
C SELECT CUT-OFF:
C
      BTWO=BTWOH                        !RE-INSTATE
C
      ITWO=1                            !INITIALIZE FOR ENERGY IND DEBYE
C
C   RCUT0.GT.0 IS ANY USER GLOBAL (IND. NL) ENERGY IND. I.P. CUT-OFF
C   FLAG RCUT0 NEGATIVE TO USE MIN OF DEBYE/DELTAE/LIFETIME
C
      IF(RCUT0.LT.DZERO)THEN
C
        RCUT0=-DEBYE                    !INITIALIZE FOR DEBYE
C
C-----------------------------------------------------------------------
C
C LIFETIME CUT-OFF BASED ON USER SUPPLIED TLIFE0 FOR THIS NL
C
C-----------------------------------------------------------------------
C
        TLIFE=MIN(TLIFE0,D1P10)
        TLIFE=TLIFE/BOHRT               !SEC TO A.U.
        DLIFE=XL*TLIFE*VRMS
C
C-----------------------------------------------------------------------
C
C ENERGY CUT-OFF BASED ON USER SUPPLIED DELTAE0 IN RYDS FOR NL->L+LD0
C
C-----------------------------------------------------------------------
C
        DELTAE=MAX(DELTAE0,D1M20)
        beta=5*deltae*n*n
        beta=5*beta*beta
        if(tee.lt.beta)then
          if(ierrout.ge.0)WRITE(ierrout,*)
     &          '*** SR.LCHNG: Delta-e too large for kT, increase kT to'
     &         ,nint(beta*tkay),' K'
          IERR=13
          GO TO 999
        endif
C
        DDELE=2*XD*VRMS/DELTAE          !DELTAE IN RYD
C
C-----------------------------------------------------------------------
C
C USE MIN OF ENERGY IND. DEBYE AND ENERGY DEP. LIFETIME/DELTA-E.
C
C-----------------------------------------------------------------------
C
        DTEST=MIN(DLIFE,DDELE)
C
        IF(DTEST*DTEST.LT.-IDTEST*RCUT0*DTEST)THEN
          if(bout)then
            if(dtest.lt.dlife)then
              write(iout,"(/' R_delta-E(RMS)=',1p,e9.2,' au',3x,1p,e8.2
     &                   ,' cm')")ddele,ddele*bohra
            else
              write(iout,"(/' R_life(RMS)   =',1p,e9.2,' au',3x,1p,e8.2
     &                   ,' cm')")dlife,dlife*bohra
            endif
          endif
          RCUT0=-DTEST                  !HOLD FLAG
          ITWO=2                        !USE E-DEP RCUT (BTWO=.T.OR.F.)
        ELSE
          BTWO=.FALSE.                  !AS PURE DEBYE
        ENDIF
C
      ELSE
C
        BTWO=.FALSE.
C
      ENDIF
C
C-----------------------------------------------------------------------
c
c define an effective deltae based-on rcut (for poss. future use.)
c      if(rcut0.ne.dzero)deltae=2*xd*vrms/abs(rcut0)
C
C-----------------------------------------------------------------------
C
C WE HAVE FINAL RCUT CUT-OFF
C
      RCUT=ABS(RCUT0)
      RCBAR=RCUT*RCUT                   !SQUARED
c
      rmean=(3*n*n-l1*(l1+1))/(2*zt)    !mean radius
      if(rcut*rcut.lt.rcut*rmean)then
        if(rcut.lt.debye*0.999D0)then   !warn
          if(ierrout.ge.0)
     &    WRITE(ierrout,"('*** SR.LCHNG: R_cut(RMS)=',1p,e9.2,' au',3x,
     &       '.LT. Mean radius=',1p,e9.2,' au'/)")rcut,rmean
        else                            !abort
          if(ierrout.ge.0)
     &    WRITE(ierrout,"('*** SR.LCHNG: R_Debye=',1p,e9.2,' au',3x,
     &       '.LT. Mean radius=',1p,e9.2,' au'/)")rcut,rmean
          IERR=14
          GO TO 999
        endif
      endif
C
C DEFINE MATCHING ENERGY AND ALPHA FOR COMBINED LIFE/DELTA & DEBYE:
C RCUT(DEBYE)=RCUT(LIFETME)
C
      IF(BTWO)THEN
        EEC=TEE*DEBYE2/RCBAR
c        write(0,*)'eec=',eec
        ALPHAC=3*ZP*N*SQRT(RMASS/EEC)/(2*ZT*DEBYE)
      ELSE                                              !JUST INITIALIZE
        IF(ITWO.EQ.1)THEN                                         !DEBYE
          EEC=DZERO
        ELSE                                               !LIFE/DELTA-E
          EEC=D1P10*TEE
        ENDIF
      ENDIF
C
C-----------------------------------------------------------------------
C
C N SET-UP
C
C-----------------------------------------------------------------------
C
      IF(N.GT.650)THEN                  !NEED TO RE-SCALE ULTRAS
        cmp0=1.d200
        if(bout)write(iout,*)
     &     '*** SR.LCHNG: n=',n,' re-scaling ultraspherical polynomials'
        IF(AC.GT.5.D-3)THEN
          if(ierrout.ge.0)WRITE(ierrout,*)
     &        '*** SR.LCHNG: ACCURACY PARAMETER TOO LARGE, TRY AC=1.D-3'
          IERR=10
          RETURN
        ENDIF
        IF(AC.GT.5.D-4)THEN
          if(bout)WRITE(iout,*)
     &                 '*** SR.LCHNG: POSSIBLE INACCURACY, TRY AC=1.D-4'
        ENDIF
      else
        cmp0=done
      ENDIF
C
      nv=n/100
      nv=nv*nv
      nv=nv/6
      tv=nv+1
      tv=log(tv)
C
      AN=N*N
      PIAN=PI*AN
      N1=N-1
      TJ=N1
      TJ=TJ/2
      N2=N-2
C
      ALLOCATE(AL(N1),STAT=JERR)
C
      IF(JERR.NE.0)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                           '*** SR.LCHNG: ALLOCATION FAILS FOR AL'
        IERR=30
        GO TO 999
      ENDIF
C
C-----------------------------------------------------------------------
C
C L SET-UP
C
C-----------------------------------------------------------------------
C
      LD=ABS(LD0)
      LD=MAX(LD,IONE)
C
      DO L=N1,LD,-1
        AL(L)=DZERO
      ENDDO
C
      IF(BREV.AND.LD0.GT.0)THEN         !IF L->L+1 THEN
        LDM=-LD                         !CALCULATE L+1->L
        L1=L1+LD                        !AND APPLY RECIPROCITY
      ELSE
        LDM=SIGN(LD,LD0)                !SO L->L+1 CASE LD0=0
        if(l10.lt.0)ldm=-ld             !test l->l-1
      ENDIF
C
      L11=2*L1+1
      W1=L11                          !INITIAL L STAT WEIGHT (QM)
      TL1=L1
      L2=L1+LDM
      L22=2*L2+1
      W2=L22                          !FINAL L STAT WEIGHT (QM)
      W12=W1/W2
      TL2=L2
      L12=L1+L2
      LGT=MAX(L1,L2)
      LLT=MIN(L1,L2)
C
      lsum=max(lsum0,5+ld*lgt/4)
      LX=MIN(L12,N1,lsum)
C
      NDIMW=LX
      ALLOCATE(W6J(-1:NDIMW+1),STAT=JERR)
C
      IF(JERR.NE.0)THEN
        if(ierrout.ge.0)WRITE(ierrout,*)
     &                          '*** SR.LCHNG: ALLOCATION FAILS FOR W6J'
        IERR=40
        GO TO 999
      ENDIF
C
C-----------------------------------------------------------------------
C
C INITIALIZE FOR QM & PS:
C
C-----------------------------------------------------------------------
C
      IF(LD.EQ.1)THEN                 !SET-UP FOR PENGELLY&SEATON
        SLINIJ=9*AN*LGT*(AN-LGT*LGT)/(2*ZT2)                   !RESOLVED
        SLINI=9*AN*(AN-L1*(L1+1)-1)/(2*ZT2)                  !UNRESOLVED
        IF(BPLTOT)THEN !SET HISTORIC SLIN FOR P&S, COULD JUST SET PPS0
          W=DONE
          SLIN=SLINI
          BR=DONE
        ELSE
          W=L11
          SLIN=SLINIJ
          IF(BMOD)THEN
            BR=DHALF                  !GUZMAN ET AL PSM
c            br=slinij/slini           !summers, storey etc
c            br=br/w
          ELSE
            BR=SLINIJ/SLINI
            BR=BR/W
          ENDIF
        ENDIF
        PPS0=2*ZP2*SLIN/(3*W)         !includes AS spin factor 2
        RONE=PPS0/BR
        RONE=RONE/PHALF               !=(R1*V)**2
        PONE=PPS0/RONE                !=PHALF*BRANCHING RATIO
        R1BAR=RONE/VRMS2              !SQUARED
        EMIN=R1BAR/RCBAR              !R1=RCUT
        if(itwo.eq.2)emin=sqrt(emin)  !energy-dependent emin
        EMIN=TEE*EMIN
        EMIND=TEE*R1BAR/DEBYE2        !DEBYE EMIN
      ENDIF
C
C-----------------------------------------------------------------------
C                                       !SET-UP FOR QUANTUM MECHANICAL
C
      CALL WIG6JR(W6J,TL1,TL2,TJ,TJ,TJ,JMIN,JMAX,IH,NDIMW)    !RECURSION
C
      if(jmin.lt.0)then                                         !failure
        if(jmin.eq.-2.and.ld0.eq.0)then               !for unresolved ps
          w6j=dzero                                             !o.k.
        else                                          !should not happen
          ierr=2
          go to 999
        endif
      endif
C
      DO L=LX,LD,-1
C
        LP=L+1
        NL=N-L
        F=GAM(LP)*GAM(LP)*GAM(NL)/GAM(N+LP)
        M1=2*JGAM(LP)+JGAM(NL)-JGAM(N+LP)
        F=(L+LP)*F*SCALE**M1
        AL(L)=F*W6J(L)*W6J(L)
C
      ENDDO
C
C-----------------------------------------------------------------------
C
C INITIAL SET-UP FOR ALPHA MESH (COMBINES VEL. AND I.P. "LOOPS")
C
      la=min(ld,5)
      if(ac.lt.1.d-10)ac=1.d-2          !reset to production
      iac=-nint(log10(ac))
      iac=max(iac,1)                    !don't allow too crude
      iamax=100*iac**2                  !no. alpha steps
      if(ld.eq.1) then                  !need P&S cut-off
        a0=1.0d0                        !min alpha (max i.p.)
        a1=720.d0*iac**3                !max alpha (min i.p.)
      else                              !to convergence?
        a0=0.5d0*la                     !min alpha
c        a1=14400.d0*la                  !max alpha (min i.p.)
        a1=720.d0*iac**3                !max alpha (min i.p.)
      endif
      iamx0=iamax
c      write(0,*)iac,iamx0,nint(a1)
C
C-----------------------------------------------------------------------
C
C CROSS SECTIONS ARE CALCULATED AT A SINGLE N-SCALED REPRESENTATIVE
C VELOCITY VREP (A.U.) FOR ALL NECESSARY IMPACT PARAMETERS (B).
C CROSS SECTIONS AT ALL OTHER VELOCITIES ARE DETERMININED FROM THIS
C GENERATED SET OF ALPHA=1/(B*VREP).
C
      VREP=5.0D0                        !0.25 FOR VRINCEANU ET AL 2012
C
 100  v0=vrep                           !min vel a.u.
      v1=v0                             !max vel a.u.
C
      IVMAX=1                           !NO. OF V/E STEPS
      VMIN=LOG(V0)                      !MIN VEL/ENERG
      VMAX=LOG(V1)                      !MAX VEL/ENERG
      DV=VMAX-VMIN
      DV=DV/IVMAX
C      IVMAX=IVMAX+1                    !UNCOMMENT FOR .GT. 1 ENERGIES
c
      if(bout8)write(iout8,"(/' n=',i4,3x,'l=',i4)")n,l1
C
      V=tv+VMIN-DV
C
      DO IV=1,IVMAX                     !START ANY VELOCITY/ENERGY-LOOP
C
        V=V+DV
C
        VS=EXP(V)
        VS2=VS*VS
        VEE2=VS2/AN                     !UNSCALED SQUARE
c         write(0,*)'vrms,vrep=',vrms,vs/n
        E=RMASS*VEE2
        ALPHA0=3*ZP/(2*VS*ZT)
c
c tbd: test that vrep is large enough so that sigma has
c      reached its uncut asymptotic form there.
c
        if(ld.gt.1.and.-deltae/e.gt.0)write(0,*)deltae/e
C
C-----------------------------------------------------------------------
C
C SET-UP ALPHA
C
        IAMAX=IAMX0
        IF(BCUT)THEN
          T=RCUT
          if(itwo.eq.2)then                  !e-dep cut-off
            t=t*vs/(n*vrms)             !i.e. rcut(e)
            iflagd=0
            if(btwo.and.t.gt.debye)then
              iflagd=1
c              if(bout)write(iout,
c     &               "(/'Reducing rcut(e) to Debye:'1p,2e10.2)")t,debye
              t=debye
            endif
          endif
          AMIN=AN*ALPHA0/T              !MIN ALPHA
c          if(btwo)write(0,*)'amin,alphac:',amin,alphac
        ELSE
          AMIN=A0/AN                    !MIN ALPHA
        ENDIF
c
        if(ld.eq.1)then
          bmax=alpha0/amin
          t=sqrt(rone/vee2)/an
c          write(0,*)t,bmax
          if(t.gt.bmax)then     !NUMERICAL PS *AND* QM WILL BE WRONG!!!
           vrep=vs*sqrt(1.2*t/bmax)
           if(bout)then
             WRITE(iout,*)
     &               '*** SR.LCHNG: R_1 .gt. b_max: ',nint(t),nint(bmax)
             WRITE(iout,*)
     &               '*** SR.LCHNG: Increasing velocity VREP to ',vrep
           endif
           go to 100                    !might need to watch for runaway
          endif
        endif
C
        AMAX=A1/AN                      !MAX ALPHA
        AMIN=LOG(AMIN)
        AMAX=LOG(AMAX)
        DA=AMAX-AMIN
        DA=DA/IAMAX
        IAMAX=IAMAX+1
C
        ALLOCATE(TQM(0:IAMAX),TPS(0:IAMAX),ENERG(0:IAMAX),STAT=JERR)
C
        IF(JERR.NE.0)THEN
          WRITE(0,*)'*** SR.LCHNG: ALLOCATION FAILS FOR IAMAX'
          IERR=50
          GO TO 999
        ENDIF
C
        AS=EXP(AMIN)
C
        if(btwo)then                    !allows for energy dep cut-off
          if(as.gt.alphac)then          !lifetime/delta-e
            if(iflagd.eq.1)stop 'e-dep/debye mix-up'
            iaea=1
          else                          !debye
            if(iflagd.eq.0)stop 'debye/e-dep mix-up'
            iaea=2
          endif
        else                            !constant, debye or v_rms
          iaea=2/itwo
        endif
C
        TE=E*AS**iaea
C
        SQM=DZERO
        PA=DZERO
        IF(LD.EQ.1)THEN
          SPS=DZERO
          PC=DZERO
        ENDIF
c
        if(bout8)then
          if(ld.eq.1)then
            write(iout8,*)'#   b     pqm    pps'
          else
            write(iout8,*)'#   b     pqm'
          endif
        endif
C
C-----------------------------------------------------------------------
C
C *** START ALPHA-LOOP
C
C-----------------------------------------------------------------------
C
c        iflag=0
        A=AMIN-DA
        IREV0=IAMAX+1
        IREVE=1
C
        DO I=1,IAMAX
C
          AS0=AS
          A=A+DA
          AS=EXP(A)
          TAS3=AS**3
          TAS0=AS-AS0
C
          ALPHA=AS
          ALPHA2=ALPHA*ALPHA
          ALPHA1=DONE+ALPHA2
          T=DPHI*SQRT(ALPHA1)
          CKHI=DONE+ALPHA2*COS(T)
          CKHI=CKHI/ALPHA1
          SKHI=DONE-CKHI*CKHI           !SQUARED
c        write(0,*)i,alpha,ckhi,skhi
C
          IREV=IREV0-I
c
c allow for combined cut-offs R(Debye) & R(E), switch over at alphac
          if(btwo)then                  !allows for energy dep cut-off
            if(as.gt.alphac)then        !lifetime/delta-e
              if(iaea.eq.2)then         !first point so switch from old
                iaea=1
                e=3*zp*n*vrms*rmass/(2*zt*rcut*alpha)  !recalculate e
                te=e*alpha              !as in general do not hit alphac
c                write(0,*)'alpha match=',alpha
              endif
            else                        !debye
              iaea=2
            endif
c          else                          !fixed, debye or v_rms
          endif
c
          ENERG(IREV)=TE/ALPHA**iaea
          IF(ENERG(IREV).GE.EMIN)IREVE=IREV    !TO CUT-OFF PS AT E=E_MIN
C
          B=ALPHA0/ALPHA                !I.P.
C
C PQM
          IM=N+N
          IP=0
C
          C0M=DZERO
          CMP=done/cmp0                 !maybe rescaling
C
          IF(N1.LE.L12.AND.AL(N1).NE.DZERO)THEN
            PQM=AL(N1)*(4*SKHI)**N1
          ELSE
            PQM=DZERO
          ENDIF
C
          DO L=N2,LD,-1                 !NOW FOR A NEAT RECURRENCE
C
            IM=IM-2
            IP=IP+1
C
            C00=CMP
            CMP=IM*(CKHI*C00-C0M)/IP
            C0M=IM*(C00-CKHI*C0M)/(IP+IM-1)
C
            IF(L.LE.L12.AND.AL(L).NE.DZERO)
     &        PQM=PQM+cmp0*CMP*AL(L)*cmp0*CMP*(4*SKHI)**L
C
          ENDDO
C
          PA0=PA
          PA=PQM/TAS3
          SQM=SQM+(PA0+PA)*TAS0         !QUICK TRAPEZOIDAL
          TQM(IREV)=SQM
C
C PPS
          IF(LD.EQ.1)THEN
            BN=B*AN
            T=BN*BN*VEE2
            IF(T.GT.RONE)THEN           !N.B. R1<RCUT
              PPS=PPS0/T
            ELSE
              PPS=PONE                  !P=PHALF*B-RATIO for R<R1
              IF(BMOD)THEN              !P=PHALF*1/2*(R/R1)**IPOW : R<R1
                T=(T/RONE)**IPOW
                PPS=PPS*SQRT(T)
              ENDIF
            ENDIF
            PC0=PC
            PC=PPS/TAS3
            SPS=SPS+(PC0+PC)*TAS0       !QUICK TRAPEZOIDAL
            TPS(IREV)=SPS
          ENDIF
C
C ALL
c
c --- output probabilities to file
c
          if(bout8)then
c            b=b*vs                     !aligns different vs-curves
c            b=b*an                     !aligns different n-curves
            if(ld.eq.1)then
c              write(iout8,*)b,b*pqm*l22,b*pps
              write(iout8,*)b,pqm*l22,pps           !,alpha,e12
            else
              write(iout8,*)b,pqm*l22               !,alpha,e12
            endif
          endif
C
        ENDDO                           !END ALPHA-LOOP
c
        if(bout8)call flush(iout8)      !case user exits ungracefully
C
C-----------------------------------------------------------------------
C
        T=B/ALPHA0                      !CONTRIB [0:B], & ADJUST NORM
        T=T*T
        t=2*t/3                         !take p(r)=p(b)*r/b for all
c        write(0,*)b*an,pqm,pqm*t,sqm
        SQM=SQM+T*PQM
        TQM(0)=SQM
        if(.not.bmod)t=3*t/2            !take p=const for P&S
        SPS=SPS+T*PPS
        TPS(0)=SPS
        ENERG(0)=DZERO
C
C-----------------------------------------------------------------------
C
C CROSS SECTIONS & COLLISION STRENGTHS:
C
C n.b. sigma~n^4 for fixed-v, here we use fixed-vs, so really ~n^6.
C
C-----------------------------------------------------------------------
C
        T0=ALPHA0*N                     !V=VS/N
        T0=T0/AN                        !REMOVE N^4 FROM SIGMA
        T0=T0*T0
C
        TEM=2*vs2*RMASS*RMASS           !2*E*RMASS*N^2
        VN4=VS2*AN*AN
C
        SQM=SQM*L22*T0                  !SIGMA/N^4   units pi*a0**2
        OQM=SQM*L11*TEM                 !OMEGA
C
        IF(LD.EQ.1)THEN                 !DIPOLE
C
          SPS=SPS*T0                    !SIGMA/N^4   units pi*a0**2
          OPS=SPS*L11*TEM               !OMEGA
C
          R1=RONE/VEE2                  !SQUARED STILL
          R2=RCBAR                      !SQUARED STILL
          if(itwo.eq.2)r2=r2*e/tee      !as rcbar energy dependent
          T1=MIN(R1,R2)                 !FOR DEBYE
C
          SPSX=PONE*T1*(PMOD+LOG(R2/T1))/AN**3             !EXACT P&S
          if(imin.eq.1)spsx=spsx*sqrt(t1/r1)
          OPSX=SPSX*L11*TEM             !OMEGA
c
          if(bout)write(iout,"(2i4,f5.2,1p,5e9.2)")
     &                       n,l1,vs,sqm,sps,spsx
c          if(bout)write(iout,"(2i4,f5.2,1p,5e9.2)")
c     &                       n,l1,vs,oqm,ops,opsx
C
        ELSEIF(LD.EQ.2)THEN             !QUADRUPOLE
C
          if(bout)write(iout,"(2i4,f5.2,1p,5e9.2)")
     &                       n,l1,vs,sqm
c          if(bout)write(iout,"(2i4,f5.2,1p,5e9.2)")
c     &                       n,l1,vs,oqm
c
        else                            !OCTUPOLE ETC.
c
          if(bout)write(iout,"(2i4,f5.2,1p,5e9.2)")
     &                       n,l1,vs,sqm
c          if(bout)write(iout,"(2i4,f5.2,1p,5e9.2)")
c     &                       n,l1,vs,oqm
        ENDIF
C
      ENDDO                             !END VEL/ENERGY-LOOP
C
C-----------------------------------------------------------------------
C
C RATE COEFFICIENTS
C
C-----------------------------------------------------------------------
C
      if(bout9)write(iout9,"(/' n=',i4,3x,'l=',i4)")n,l1
C
      IF(BCUT)THEN                      !WE NEED TO CONVOLUTE QM (& PS)
C
        IREV=0
        TCQ=CQQ/(2*TEE)                 !CONVERT TO RATE COEFF
C
        QQM=0                           !INIT QM
        OQ2=0
        TQSUM=TQM(IREV)
C        IREV=IREVE                     !CUT-OFF QM AT E_MIN AS WELL
        TQM(IREV)=0
        TQ=L22*T0*VN4                   !CONVERT TO E*SIGMA
C
        IF(LD.EQ.1)THEN                 !***DIPOLE
C
          QPS=0                         !INIT P&S
          OP2=0
          TPSUM=TPS(IREV)
C          IREV=IREVE                   !CUT-OFF PS AT E_MIN AS WELL
          TPS(IREV)=0
          TP=T0*VN4                     !CONVERT TO E*SIGMA
C
          if(bout9)write(iout9,*)
     &     '# i,     e,     sqm,     sps,     spsx'
c
c         write(0,*)itwo,emin
          DO I=IREV+1,IAMAX             !CONTRIB FROM ~0 TO E=M*V^2
            TQM(I)=TQSUM-TQM(I)         !REVERSE SUM SINCE WE SUMMED
            TPS(I)=TPSUM-TPS(I)         !FROM ALPHA MIN TO ALPHA MAX
            IF(ENERG(I).LT.EMIN)THEN
              IF(.NOT.BMIN)THEN         !STANDARD P&S NEGLECT E<E_MIN
                if(i+1.eq.ireve)then        !energ(i)<e_min<energ(ireve)
                  p1=tps(i)
                  p2=tpsum-tps(ireve)
                  h=energ(ireve)-energ(i)
                  hm=(energ(ireve)-emin)
                  pm=((emin-energ(i))*p2+hm*p1)/h         !assume linear
                  y=energ(ireve)/tee
                  p2=p2*exp(-y)
                  y=emin/tee
                  pm=pm*exp(-y)
                  qps=qps+(p2+pm)*hm      !contrib from emin to e(ireve)
                  qps=qps-p2*h                       !avoid double count
                endif
                TPS(I)=DZERO
              ELSEIF(IMIN.EQ.2)THEN     !CORRECT R/R_1->R/R_c, R<R_c
                t=(EMIN/ENERG(I))**IPOW
                if(itwo.eq.1)t=SQRT(t)
                TPS(I)=TPS(I)*t
              ENDIF
            ENDIF
c
c --- output cross sections to file
c
            if(bout9)then
              sqm9=rmass*tq*tqm(i)/energ(i)
              sps9=rmass*tp*tps(i)/energ(i)
              if(ld.eq.1.and.(energ(i).gt.emin.or.bmin))then
                r1=rone*rmass/energ(i)                    !squared still
                r2=rcbar                                  !squared still
                if(itwo.eq.2)r2=r2*energ(i)/tee
                t1=min(r1,r2)                             !for debye
                spsx9=pone*t1*(pmod+log(r2/t1))           !exact p&s
                if(imin.eq.1)spsx9=spsx9*sqrt(t1/r1)
              else
                spsx9=0
              endif
c
c small E scaling
              if(imin.eq.-1)then
                if(energ(i).gt.eec)then                    !debye
                  t=sqrt(energ(i))
                else
                  t=energ(i)*energ(i)
                endif
              else
                if(energ(i).gt.eec)then                    !debye
                  t=done
                else
                  t=energ(i)
                endif
              endif
c
c large E scaling
c              t=-log(energ(i))
c
c global factors
              tt=vn4               !*exp(-energ(i)/tee)
c x-axis
              e=energ(i)
c              e=e/tee
c              if(energ(i).gt.eec)e=sqrt(e)
c              b=alpha0*e
c              e=b
c
c --- write cross sections to file
c
              write(iout9,*)i,e
     &          ,sqm9,sps9,spsx9
c     &          ,tqm(i)*tq/tt,tps(i)*tp/tt
c     &          ,sqm9/t,sps9/t,spsx9/t
c     &          ,sps9/spsx9,emin/energ(i),r1,r2
            endif
c
            H=ENERG(I)-ENERG(I-1)
            Y=ENERG(I)/TEE
            EY=EXP(-Y)
            OQ1=OQ2
            OP1=OP2
            OQ2=TQM(I)*EY
            OP2=TPS(I)*EY
            QQM=QQM+(OQ1+OQ2)*H
            QPS=QPS+(OP1+OP2)*H
          ENDDO
C
          if(bout9)call flush(iout9)       !case user exits ungracefully
c
          QQM=TCQ*QQM*TQ
          QPS=TCQ*QPS*TP
C
C-----------------------------------------------------------------------
C
C ASYPTOTIC ENERGY CONTRIBUTION
C
C-----------------------------------------------------------------------
C
          EJI=ENERG(IAMAX)
c
          if(eji.lt.emin)stop'asymptotic energy .lt. emin...'
C
          JTWO=ITWO
          T2=ITWO*LOG(EJI/EMIN)
          IF(BTWO)THEN
            IF(EJI.GE.EEC)THEN
              JTWO=1
              T2=LOG(EJI/EMIND)
            ENDIF
          ENDIF
C
          T=TQ*TQM(IAMAX)
          PQM0=PPS0                             !ASSUME (AS PPS0 EXACT)
c          pqm0=t/(pmodq+t2)                     !to define
          PMODQ=T/PQM0-T2                       !EFFECTIVE PMODQ
C
C SIGMA(E)=INTEGRAL OVER PROBS ONLY COMPLETE FOR E.LE.EJI. SO, NOW USE
C ANALYTIC PS TO ADD-IN REMAINING PROBS TO QM & PS,
C FOR E=EJI TO EJF, AND ADD THE CONTRIBUTION FROM THESE SIGMA(E) TO THE
C CORRESPONDING RATES.
C
          EJF=5*TEE                     !SAY
          H=0.04
          NS0=LOG(EJF/EJI)/H
          NS0=NS0+1
          NS0=2*((NS0+1)/2)+1           !NEED ODD
C
          IF(NS0.LT.3)THEN
            E2=EJI
          ELSE
C
            X1=LOG(EJI)
            X=X1-H
            IFAC=1
            OQ=TQSUM*TQ                 !CONVERT TO E*SIGMA
            OP=TPSUM*TP                 !CONVERT TO E*SIGMA
            TQSUM=DZERO
            TPSUM=DZERO
            E2=EJI
C
            DO NS=1,NS0
              E1=E2
              X=X+H
              E2=EXP(X)
              Y=E2/TEE
              if(e2.ge.eec)jtwo=1
              TE=LOG(E2/E1)*jtwo
              T=PQM0*TE                          !NORMALLY SET PQM0=PPS0
              OQ=OQ+T                            !SO TOP-UP QM WITH PS
              T=PPS0*TE                          !AS PS USES EXACT SLIN
              OP=OP+T                            !TOP-UP PS
              T=IFAC*EXP(-Y)*E2
              TQSUM=TQSUM+T*OQ
              TPSUM=TPSUM+T*OP
              IF(IFAC.LT.3)THEN
                IFAC=4
              ELSE
                IFAC=2
              ENDIF
            ENDDO
C
            T=T/IFAC
            TQSUM=TQSUM-T*OQ
            TPSUM=TPSUM-T*OP
            T=CQQ*H/(3*TEE)             !CONVERT TO RATE COEFF
            QTOPQ=T*TQSUM
            QTOPP=T*TPSUM
c         write(0,*)'TOP-UP FROM-TO:',eji,ejf,' To '
c     &             ,qqm,qps,'add top-ups:',qtopq,qtopp
            QQM=QQM+QTOPQ
            QPS=QPS+QTOPP
C
          ENDIF
C
C NOW ADD-IN CONTRIBUTION TO (NUMERICAL) RATES FROM SIGMA(E) FOR E=EJF
C TO INFINITY USING ANALYTIC PS SIGMA(E).
C
          EJF=E2
C
          EMIN=EMIN/TEE                            ! /KT NOW
c        write(0,*)'Emin/kT=',emin
          EMAX=EJF/TEE                             !DITTO
          TEXP=EXP(-EMAX)
C
          QTOPEP=(PMOD+ITWO*LOG(EMAX/EMIN))*TEXP
     &                               +ITWO*ENINT(IONE,.FALSE.,EMAX)
C
          IF(BTWO)THEN
            EEC=EEC/TEE
            IF(EEC.GT.EMAX)THEN
              QTOPEP=QTOPEP-ENINT(IONE,.FALSE.,EEC)
            ELSE
              QTOPEP=QTOPEP-ENINT(IONE,.FALSE.,EMAX)
            ENDIF
            EMIND=EMIND/TEE
          ENDIF
C
          QTOPEQ=QTOPEP+(PMODQ-PMOD)*TEXP
          QTOPEQ=CQQ*PQM0*QTOPEQ
          QTOPEP=CQQ*PPS0*QTOPEP
c       write(0,*)'TOP-UP FROM (TO INF)=',ejf,'; To '
c     &           ,qqm,qps,'add top-up:',qtopeq,qtopep
          QQM=QQM+QTOPEQ                !TOP-UP QM WITH PS NORMALLY
          QPS=QPS+QTOPEP                !AS PS USES EXACT SLIN
C
        ELSE                            !*** NON-DIPOLE WITH RCUT
C
          DO I=IREV+1,IAMAX             !CONTRIB FROM ~0 TO E=M*V^2
            TQM(I)=TQSUM-TQM(I)         !REVERSE SUM SINCE WE SUMMED
c
            if(bout9)then
              sqm9=rmass*tq*tqm(i)/energ(i)
              write(iout9,*)i,energ(i),sqm9     !pi*a0**2
            endif
c
            H=ENERG(I)-ENERG(I-1)
            Y=ENERG(I)/TEE
            EY=EXP(-Y)
            OQ1=OQ2
            OQ2=TQM(I)*EY
            QQM=QQM+(OQ1+OQ2)*H
          ENDDO
C
          if(bout9)call flush(iout9)    !case user exits ungracefully
c
          QQM=TCQ*QQM*TQ
C
          QQMA=CQQ*SQM*VN4              !ASYMPTOTIC VALUE
          EMAX=ENERG(IAMAX)/TEE
          QTOPE=QQMA*EXP(-EMAX)
c       write(0,*)'TOP-UP FROM (TO INF)=',energ(iamax),'; To '
c     &           ,qqm,'add top-up:',qtope
          QQM=QQM+QTOPE
c       write(0,*)'RATIO QM TO ASYMP:',qqm/qqma
        ENDIF
C
      ELSE                              !*** NON-DIPOLE ***
        QQM=CQQ*SQM*VN4                 !ASSUMING CAN IGNORE RCUT
      ENDIF
C
      UQM=QQM*TCM*L11                   !UPSILON
C
      IF(LDM*LD0.LT.IZERO)THEN         !APPLY RECIPROCITY FOR PRODUCTION
        W12=W1/W2
        QQM=QQM*W12
      ENDIF
C
C-----------------------------------------------------------------------
C
C NOW EVALUATE ANALYTIC PS RATE COEFFICIENTS
C
C-----------------------------------------------------------------------
C
      IF(LD.EQ.1)THEN                                            !DIPOLE
C                                       
C CONTRIB E_MIN->\infty
C
        E1=ENINT(IONE,.FALSE.,EMIN)     !USE EXACT FOR E_MIN NOT SMALL
C
        TX=EXP(-EMIN)
        TT=PMOD*TX+ITWO*E1              !DEBYE OR LIFE/DELTA-E, ALL PS
C
        IF(BTWO)THEN                    !COMBINED
C          EEC=DEBYE2/RCBAR              !E/KT:RCUT(DEBYE)=RCUT(LIFETME)
          IF(EEC.GT.EMIN)THEN           !ALL LIFE/DELTA-E
            TT=TT-ENINT(IONE,.FALSE.,EEC)
          ELSE                          !TRUNCATE BY DEBYE (UNLIKELY)
C            EMIND=R1BAR/DEBYE2          !DEBYE EMIN/KT
            TX=EXP(-EMIND)
            TT=PMOD*TX+ENINT(IONE,.FALSE.,EMIND)
          ENDIF
        ENDIF
C
C CONTRIB 0->E_MIN
C
        IF(BMIN)THEN
C
          IF(IMIN.EQ.1)THEN
C
C PSM20
            IF(BTWO)THEN                               !COMBINED
C
              IF(EEC.LT.EMIN)THEN          !TRUNCATE BY DEBYE (UNLIKELY)
c
c                write(0,*)"Note, matching energy  E_c .lt. E_min"
c     &                                                 ,eec,emin
                TC=EXP(-EEC)
                T=6*(DONE-TC*(DONE+(DONE+EEC/4)*EEC))/EMIN**3
                RTED=SQRT(EMIND)
                RTEC=SQRT(EEC)
                T0=D3QRT*RTPI*(ERF(RTED)-ERF(RTEC))
                T0=T0/RTED-TX*(D3HALF+EMIND)
                T0=T0/EMIND+T
C
              ELSE                                     !ALL LIFE/DELTA-E
C
                T0=(DONE-TX*(DONE+EMIN))/EMIN
                T0=DTWO*T0/EMIN-TX
                T0=DTHREE*T0/EMIN-TX
c
c                write(0,*)"Note, matching energy  E_c .gt. E_min"
c     &                                                  ,eec,emin
              ENDIF
C
            ELSE                                       !ONE OR OTHER
C
              IF(ITWO.EQ.1)THEN                        !DEBYE
C
                RTEM=SQRT(EMIN)
                T0=D3QRT*RTPI*ERF(RTEM)/RTEM-TX*(D3HALF+EMIN)
                T0=T0/EMIN
C
              ELSE                                     !LIFE/DELTA-E
C
                T0=(DONE-TX*(DONE+EMIN))/EMIN
                T0=DTWO*T0/EMIN-TX
                T0=DTHREE*T0/EMIN-TX
C
              ENDIF
C
            ENDIF
C
          ELSE
C
C PSM17
C
            IF(BTWO)THEN                               !COMBINED
C
              IF(EEC.LT.EMIN)THEN          !TRUNCATE BY DEBYE (UNLIKELY)
c
c                write(0,*)"Note, matching energy  E_c .lt. E_min"
c     &                                                 ,eec,emin
                TC=EXP(-EEC)
                T=(DTWO-TC*(DTWO+EEC))/EMIN
                T0=T/EMIN-TX*(DONE+DONE/EMIND)
C
              ELSE                                     !ALL LIFE/DELTA-E
C
                T0=(DONE-TX*(DONE+EMIN))/EMIN
                T0=DTWO*T0/EMIN-TX
c
c                write(0,*)"Note, matching energy  E_c .gt. E_min"
c     &                                                  ,eec,emin
              ENDIF
C
            ELSE                                       !ONE OR OTHER
C
              T0=(DONE-TX*(DONE+EMIN))/EMIN            !DEBYE
              IF(ITWO.EQ.2)T0=DTWO*T0/EMIN-TX          !LIFE/DELTA-E
C
            ENDIF
C
          ENDIF
C
          T0=PMOD*T0
          TT=TT+T0
C
        ELSE
C
C STANDARD P&S OMITS 0-E_MIN, SYNC WITH NUMERICAL ABOVE
c
c also original PS approximates of exponential integral
c
          if(.not.bexp)tt=pmod-itwo*(euler+log(emin))
C
        ENDIF
C
        QPSX=CQQ*PPS0*TT
        UPSX=QPSX*TCM*L11
        UPS=QPS*TCM*L11
C
        IF(LDM*LD0.LT.IZERO)THEN       !APPLY RECIPROCITY FOR PRODUCTION
          QPSX=QPSX*W12
          QPS=QPS*W12
        ENDIF
c
        if(bout)write(iout,"(13x,1p,5e9.2)")qqm,qps,qpsx
c        if(bout)write(iout,"(13x,1p,5e9.2)")uqm,ups,upsx
C
      ELSEIF(LD.EQ.2)THEN                                    !QUADRUPOLE
C
        if(bout)write(iout,"(13x,1p,4e9.2)")qqm
c        if(bout)write(iout,"(13x,1p,4e9.2)")uqm
      else                                                    !OCTUPOLE+
        if(bout)write(iout,"(13x,1p,3e9.2)")qqm
c        if(bout)write(iout,"(13x,1p,3e9.2)")uqm
      ENDIF
c
      if(bout)call flush(iout)            !case user exits ungracefully
C
C-----------------------------------------------------------------------
C-----------------------------------------------------------------------
C
 999  CONTINUE
C
      IF(ALLOCATED(TQM))THEN
        DEALLOCATE(TQM,TPS,ENERG,STAT=JERR)
        IF(JERR.NE.0)THEN
          WRITE(0,*)'*** SR.LCHNG:DE-ALLOCATION FAILS FOR TQM,TPS,ENERG'
          IERR=51
        ENDIF
      ENDIF
C
      IF(ALLOCATED(W6J))THEN
        DEALLOCATE(W6J,STAT=JERR)
        IF(JERR.NE.0)THEN
          WRITE(0,*)'*** SR.LCHNG: DE-ALLOCATION FAILS FOR W6J'
          IERR=41
        ENDIF
      ENDIF
C
      IF(ALLOCATED(AL))THEN
        DEALLOCATE(AL,STAT=JERR)
        IF(JERR.NE.0)THEN
          WRITE(0,*)'*** SR.LCHNG: DE-ALLOCATION FAILS FOR AL'
          IERR=31
        ENDIF
      ENDIF
C
      RETURN
C
      END SUBROUTINE LCHNG
C
C***********************************************************************
C
      REAL*8 FUNCTION ENINT(N,BHALF,X)
C
C-----------------------------------------------------------------------
C
C EVALUATES EXPONENTIAL INTEGRAL
C N R BADNELL: ADAPTED FROM NUMERICAL RECIPES AND ABRAMOWITZ & STEGUN
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
      LOGICAL BHALF
C
      PARAMETER (DZERO=0.0D0)
      PARAMETER (DONE=1.0D0)
      PARAMETER (DTWO=2.0D0)
      PARAMETER (BIG=1.D30)
      PARAMETER (TOLE=1.D-5)
      PARAMETER (MAXIT=100)
      PARAMETER (EULER=0.5772156649D0)
C
C GENERATE INTEGER (BHALF.EQ.FALSE.) EXPONENTIAL INTEGRAL ORDER N
C GENERATE HALF-INTEGER (BHALF.EQ.TRUE.) EXPONENTIAL INTEGRAL ORDER N+.5
C
      IF(N.LT.0.OR.X.LT.DZERO.OR.X.EQ.DZERO.AND.(N.EQ.0.OR.N.EQ.1))THEN
        WRITE(0,*)'ILLEGAL ARGUMENT IN FN.ENINT: N,X=',N,X
        ENINT=DZERO
        RETURN
      ENDIF
C
      NM=N-1
C
      H=DZERO
      IF(BHALF)H=DONE/DTWO
C
      IF(N.EQ.0.AND..NOT.BHALF)THEN
        EXPINT=EXP(-X)/X
      ELSEIF(X.EQ.DZERO)THEN
        EXPINT=DONE/(NM+H)
      ELSEIF(X.GT.DONE)THEN              !CONTINUED FRACTION
        BX=X+N+H
        CX=BIG
        DX=DONE/BX
        HX=DX
        DO I=1,MAXIT
          AX=-I*(NM+H+I)
          BX=BX+DTWO
          DX=DONE/(AX*DX+BX)
          CX=BX+AX/CX
          DEL=CX*DX
          HX=HX*DEL
          IF(ABS(DEL-DONE).LT.TOLE)THEN
             EXPINT=HX*EXP(-X)
             GO TO 100
          ENDIF
        ENDDO
c        write(0,*)'FN.ENINT continued fraction not converged'
        EXPINT=-DONE
      ELSE                               !SERIES
        FACT=DONE
        IF(BHALF)THEN
          PI=ACOS(-DONE)
          FACH=-SQRT(PI*X)/(N-H)
        ENDIF
        IF(NM.LE.0)THEN
          IF(BHALF)THEN
            IF(NM.EQ.0)THEN
              EXPINT=FACH+DTWO
            ELSE
              EXPINT=SQRT(PI/X)-DTWO
            ENDIF
          ELSE
            EXPINT=-LOG(X)-EULER
          ENDIF
        ELSE
          EXPINT=DONE/(NM+H)
        ENDIF
        DO I=1,MAXIT
          FACT=-FACT*X/I
          IF(BHALF)FACH=-FACH*X/(I-H)
          IF(I.NE.NM)THEN
            DEL=-FACT/(I-NM-H)
          ELSE
            IF(BHALF)THEN
              DEL=FACH+DTWO*FACT
            ELSE              !DIGAMMA
              PSI=-EULER
              DO II=1,NM
                PSI=PSI+DONE/II
              ENDDO
              DEL=FACT*(-LOG(X)+PSI)
            ENDIF
          ENDIF
          EXPINT=EXPINT+DEL
          IF(ABS(DEL).LT.ABS(EXPINT)*TOLE)GO TO 100
        ENDDO
c       write(0,*)'ENINT series not converged'
        EXPINT=-DONE
      ENDIF
C
  100 ENINT=EXPINT
C
      RETURN
C
      END FUNCTION ENINT
C
C***********************************************************************
C
      SUBROUTINE GAMAF(GAM,JGAM,NMAX)
C
C-----------------------------------------------------------------------
C
C ALAN BURGESS, D.A.M.T.P. CAMBRIDGE.
C
C CALCULATES GAMMA(N), N=1,2...NMAX.
C OUTPUT IS IN GAM(N), JGAM(N), WHERE GAMMA(N)=GAM(N)*SC1**JGAM(N).
C
C LAST CHANGED BY NRB ON 22 JUL 21 TO REMOVE COMMON BLOCK USAGE
C
C-----------------------------------------------------------------------
C
      IMPLICIT REAL*8 (A-H,O-Z)
      IMPLICIT INTEGER*4 (I-N)
C
      DIMENSION GAM(NMAX),JGAM(NMAX)
C
      PARAMETER (DONE=1.0D0)
      PARAMETER (SC1=1.0D6)
      PARAMETER (SC2=1.0D-6)
C
      GAM(1)=DONE
      JGAM(1)=0
C
      DO N=2,NMAX
        N1=N-1
        X=N1
        X1=X*GAM(N1)
        J1=JGAM(N1)
  100   IF(X1.GT.SC1)THEN
          X1=SC2*X1
          J1=J1+1
          GO TO 100
        ENDIF
        GAM(N)=X1
        JGAM(N)=J1
      ENDDO
C
      RETURN
C
      END SUBROUTINE GAMAF
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 BADNELL AT AL MNRAS (2021)
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/.false./,debug2/.false./
      data iwritd/0/     !=0 to screen; set < 0 for non-interactive use
C
      E(A1)=SQRT(
     &           (A1*A1-(A2-A3)*(A2-A3))*((A2+A3+1)*(A2+A3+1)-A1*A1)*
     &           (A1*A1-(B2-B3)*(B2-B3))*((B2+B3+1)*(B2+B3+1)-A1*A1)
     &           )
c     &           (A1-A2+A3)*(A1+A2-A3)*(A2+A3+1-A1)*(A2+A3+1+A1)*
c     &           (A1-B2+B3)*(A1+B2-B3)*(B2+B3+1-A1)*(B2+B3+1+A1)
c     &           )
C
      F(A1)=(A1+A1+1)*(
     &           A1*(A1+1)*(-A1*(A1+1)+A2*(A2+1)+A3*(A3+1)-2*B1*(B1+1))+
     &           B2*(B2+1)*(A1*(A1+1)+A2*(A2+1)-A3*(A3+1))+B3*(B3+1)*
     &          (A1*(A1+1)-A2*(A2+1)+A3*(A3+1))
     &                 )
c     &           A1*(A1+1)*(-A1*(A1+1)+A2*(A2+1)+A3*(A3+1)+
c     &           B2*(B2+1)+B3*(B3+1)-2*B1*(B1+1))+
c     &          (A2*(A2+1)-A3*(A3+1))*(B2*(B2+1)-B3*(B3+1))
c     &                 )
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(
     &   NINT(A2+B1+B3-DEPS2).NE.NINT(A2+B1+B3+DEPS2)
     &                       .OR.
     &   NINT(A3+B1+B2-DEPS2).NE.NINT(A3+B1+B2+DEPS2)
     &                                               )THEN
        if(debug2)write(0,*)'*** sr.wig6jr: non-integer triangle sum...'
        JMIN=-1
        RETURN
      ENDIF
      IF(
     &   B3.GT.A2+B1+DEPS2.OR.B3.LT.ABS(A2-B1)-DEPS2
     &                    .OR.
     &   B2.GT.A3+B1+DEPS2.OR.B2.LT.ABS(A3-B1)-DEPS2
     &                                              )THEN
        if(debug2)write(0,*)'*** 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)
     &  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(0,*)'*** 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)-
     &  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)
     &  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
      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
     &               .AND.JMAX-JMIN+1.GT.2
     &                                )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(0,*)'*** 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 'jmin=0, jmid1>0...'
        JMID1=1
        JMD1M=JMID1-1
        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)
     &  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(0,*)'*** 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)
     &  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(0,*)'nodes expected/found=',nodes,nod
c      write(0,*)'wmax=',wmax                               !order unity
c
      if(debug2.and.jmid1.gt.jmid2)
     &  write(0,*)'jmid1,2=',jmid1,jmid2,jmin,jmax
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***********************************************************************