SUBROUTINE NumerovProp(nmax, mu, hstep, v, Emin, Emax, NEnergy, rvals, Psi, phase)
USE Numeric_Kinds_Module
USE Numbers_Module
USE FileUnits_Common_Module
USE FileUnits_OneDim_Module
USE QState_Module
IMPLICIT NONE
INTEGER n, nmax, KEnergy, NEnergy, j, npi
REAL(KIND=WP_Kind) E, mu, hstep, tn, qn, Scale, Emin, Emax, DeltaE, PlotPsi, SCLength
REAL(KIND=WP_Kind) z, zjn, zjnp, zyn, zynp, a, b, kmat, k, c, d, rmat, f, dphase, phaseguess
REAL(KIND=WP_Kind) Psi(0:nmax), v(0:nmax), rvals(0:nmax), phase(0:nmax)
REAL(KIND=WP_Kind), ALLOCATABLE:: u(:)
COMPLEX(KIND=WP_Kind) smat
COMPLEX(KIND=WP_Kind), PARAMETER:: i=(zero,one), zone=(one,zero)

IF(ALLOCATED(u))THEN
  IF(SIZE(Psi).ne.nmax+1)STOP "Num"
ELSE
   ALLOCATE(u(0:nmax))
ENDIF

OPEN(Unit=Numerov_Prop_Unit,File=TRIM(OutDIR)//"GraphicsOut/Numerov_Prop.csv")
OPEN(Unit=Continuum_Unit,File=TRIM(OutDIR)//"GraphicsOut/Continuum.csv")
OPEN(Unit=ContinuumBin_Unit,File=TRIM(OutDIR)//"BinOut/Continuum.bin",FORM="unformatted")
OPEN(UNIT=Psi_Unit,File=TRIM(OutDir)//"GraphicsOut/Psi.csv") 
OPEN(UNIT=RMat_Unit,File=TRIM(OutDir)//"GraphicsOut/RMat.csv") 
OPEN(UNIT=Phase_Unit,File=TRIM(OutDir)//"GraphicsOut/Phase.csv") 
OPEN(UNIT=DPhase_Unit,File=TRIM(OutDir)//"GraphicsOut/DPhase.csv") 
OPEN(UNIT=SMat_Unit,File=TRIM(OutDir)//"GraphicsOut/SMat.csv")

WRITE(Numerov_Prop_Unit,'(6A15)')"rvals,  ", "Psi,  ", "Sin,  "
WRITE(Continuum_Unit,'(A5,6A15)')"jval,  ", "E,  ", "K,  ", "Real(S),  ", "Imag(S),  ", "Phase,  ", "SCLength"
WRITE(Out_Unit,'(A5,6A15)')" jval  ", "E", "K", "Real(S)", "Imag(S)", "Phase", "SCLength"
WRITE(Psi_Unit,'(6A15)')"rvals,  ", "Psi,  "
WRITE(RMat_Unit,'(6A15)')"rvals,  ", "RMat,  ", "KMat,  "
WRITE(SMat_Unit,'(6A15)')"rvals,  ", "SReal,  ", "SImag,  "
WRITE(Phase_Unit,'(6A15)')"rvals,  ", "Phase,  "

DeltaE=(EMax-EMin)/NEnergy
DO KEnergy=0,NEnergy-1
   E=EMin+DeltaE*KEnergy
   k=SQRT(two*mu*E)
   WRITE(Screen_Unit,'(2(A,e15.7))')" E= ",E," k= ",k
   u=Zero
   Psi=Zero
   u(0)=Zero
   u(1)=One
   Psi(0)=0.d0
   WRITE(Psi_Unit,'(5(e15.8,","))')rvals(0),psi(0)
   phase(0)=0.d0
   DO n=1,nmax              ! Outward propagation using enhanced renormalized numerov method
      tn=(hstep**2/Twelve)*(Two*mu)*(v(n)-E)
      IF(tn.ge.Zero)THEN
         IF(SQRT(Twelve*tn)<700.d0)THEN
            qn=Two*COSH(SQRT(Twelve*tn))
         ELSE
            qn=One
            u(n-1)=Zero
            u(n)=one
         ENDIF
      ELSE
         qn=Two*COS(SQRT(-Twelve*tn))
      ENDIF
      IF(n<nmax)u(n+1)=qn*u(n)-u(n-1)
      psi(n)=u(n)*qn/(Two+Ten*tn)
      IF(u(n).gt.1.d+20)THEN    ! If wavefuntion gets too large rescale the wavefunction
         Scale=u(n)
         u=u/Scale
         Psi=Psi/Scale
      ENDIF
      WRITE(Psi_Unit,'(5(e15.8,","))')rvals(n),psi(n)
      k=SQRT(two*mu*E)
      z=k*rvals(n)
      CALL rbes (jval,z,zjn,zjnp,zyn,zynp)  ! Get asymptotic wave function at current rval
      a=zjn/SQRT(k)
      b=zyn/SQRT(k)
      z=k*rvals(n-1)
      CALL rbes (jval,z,zjn,zjnp,zyn,zynp)  ! Get asymptotic wave function at previous rval
      c=zjn/SQRT(k)
      d=zyn/SQRT(k)
      IF(n>1)THEN
         rmat=psi(n)/psi(n-1)
         kmat=(a-c*rmat)/(b-d*rmat)
      ELSE
         rmat=500.d0
         kmat=c/d
      ENDIF
      smat=(zone+i*kmat)/(zone-i*kmat)           ! Single channel S-Matrix
      phase(n)=ATAN2(IMAG(smat),REAL(smat))/two  ! Phase shift
      WRITE(RMat_Unit,'(5(e15.8,","))')rvals(n),rmat,kmat
      WRITE(SMat_Unit,'(5(e15.8,","))')rvals(n),smat
      SCLength=-Tan(Phase(n))/k                  ! Scattering length in atomic units
      WRITE(Phase_Unit,'(6(es15.7,","))')rvals(n), phase(n)
   ENDDO
   
   k=SQRT(two*mu*E)
   z=k*rvals(nmax-1)
   CALL rbes (jval,z,zjn,zjnp,zyn,zynp)
   a=zjn/SQRT(k)
   b=zyn/SQRT(k)
   z=k*rvals(nmax-2)
   CALL rbes (jval,z,zjn,zjnp,zyn,zynp)
   c=zjn/SQRT(k)
   d=zyn/SQRT(k)
   
   rmat=psi(nmax-1)/psi(nmax-2)
   kmat=(a-c*rmat)/(b-d*rmat)
   smat=(zone+i*kmat)/(zone-i*kmat)
   phase(nmax)=ATAN2(IMAG(smat),REAL(smat))/two
   SCLength=-Tan(Phase(nmax))/k
   WRITE(Continuum_Unit,'(i5,",",6(es15.7,","))')jval, E, kmat, smat, phase(nmax), SCLength
   WRITE(Out_Unit,'(i5,6(es15.7))')jval, E, kmat, smat, phase(nmax), SCLength
   
   f=(a*d-c*b)/(psi(nmax-1)*d-psi(nmax-2)*b)
   psi=psi*f
   DO n=0,nmax
      scale=psi(1)/Sin(k*rvals(1))
      PlotPsi=Psi(n)
      IF(ABS(PlotPsi)<1.d-30)PlotPsi=Zero
      WRITE(Numerov_Prop_Unit,'(ES15.7,",",ES15.7,",",ES15.7)')rvals(n), PlotPsi, Sin(k*rvals(n))/SQRT(k)
   ENDDO 
   
   IF(kenergy/=NEnergy-1)WRITE(Numerov_Prop_Unit,*)'&'
   IF(Kenergy==0)THEN
     WRITE(ContinuumBin_Unit)rvals
     WRITE(ContinuumBin_Unit)v
   ENDIF
   WRITE(ContinuumBin_Unit)jval, E, kmat, smat, phase(nmax)
   WRITE(ContinuumBin_Unit)psi
   
   WRITE(dphase_unit,*)"rvals,  ", "Phase1,  ", "Phase2,  ", "DPhase,  " 
   phaseguess=-k*rvals(0)
   phase(0)=phaseguess
   WRITE(dphase_unit,'(7(e15.7,","))')0.d0,0.d0,0.d0,hstep*(-k)
   WRITE(dphase_unit,'(7(e15.7,","))')rvals(0),phase(0)/pi,phaseguess/pi,hstep*(-k)   !,k**2-Two*mu*v(n)
   
   DO n=1,nmax
      k=SQRT(two*mu*E)
      z=k*rvals(n)
      CALL rbes (jval,z,zjn,zjnp,zyn,zynp)
      dphase=-Two*mu*(v(n-1)/k)*(cos(phase(n-1))*zjn-sin(phase(n-1))*zyn)**2 !Derivative of the phase shift with respect to r
      phaseguess=phase(n-1)+dphase*hstep
      IF(ABS(phaseguess-phase(n))/Pi>.9)THEN
         write(Screen_Unit,'(3e15.7)')(phaseguess-phase(n))/Pi,phase(n)/pi,phaseguess/pi
         npi=nint((phaseguess-phase(n))/Pi)
         DO j=n,nmax
            phase(j)=phase(j)+npi*pi
         ENDDO
      ENDIF
      WRITE(dphase_unit,'(7(e15.7,","))')rvals(n),phase(n)/pi,phaseguess/pi,hstep*dphase   !,k**2-Two*mu*v(n)
   ENDDO
ENDDO

CLOSE(Unit=Numerov_Prop_Unit)
CLOSE(Unit=Continuum_Unit)
CLOSE(Unit=ContinuumBin_Unit)
CLOSE(UNIT=Psi_Unit) 
CLOSE(UNIT=RMat_Unit) 
CLOSE(UNIT=SMat_Unit)
CLOSE(UNIT=Phase_Unit)
CLOSE(UNIT=DPhase_Unit)

RETURN
ENDSUBROUTINE NumerovProp