Program Test_Fit
USE Numeric_Kinds_Module
USE DAFSigma_Module
IMPLICIT NONE
CHARACTER(LEN=1),   PARAMETER:: Jobz='V', UpLo='L'
CHARACTER(LEN=27), PARAMETER:: OutDir='E:\ParkerE\DAF_And_Numerov\'
CHARACTER(LEN=255) Potential_Name
INTEGER, PARAMETER:: N_Knot=500, Mval=32, NDAF=2, nd_DAF=38, ninterp=500
INTEGER, PARAMETER:: Nlocal=3, lwork=100000
INTEGER i, j, kmax, jmin, jmax, iop(2), K, info, S_work(lwork), IPIV(lwork)
REAL(KIND=WP_Kind), PARAMETER:: a=3.0, b=40.0
REAL(KIND=WP_Kind) X_KNot(N_Knot), F_Knot(N_Knot), F_Knotp(N_Knot), F_Knotpp(N_Knot), dx, tab(3), w(N_Knot), at(N_Knot), bt(N_Knot), ct(N_Knot), xbar(ninterp), Eig(N_Knot)
REAL(KIND=WP_Kind) X_Interp(Ninterp), xdif, expfac, fac, xarg, twopi, delta2, Potential, r
REAL(KIND=WP_Kind) Gx(Ninterp), Gxp(Ninterp), Gxpp(Ninterp)
REAL(KIND=WP_Kind) Gl(Ninterp), Glp(Ninterp), Glpp(Ninterp)
REAL(KIND=WP_Kind) Gs(Ninterp), Gsp(Ninterp), Gspp(Ninterp)
REAL(KIND=WP_Kind) Gd(Ninterp), Gdp(Ninterp), Gdpp(Ninterp)
REAL(KIND=WP_Kind) F(Ninterp), FP(Ninterp), FPP(Ninterp)
REAL(KIND=WP_Kind) DAFx(0:nd_daf,0:NDAF), hermit(0:Mval+NDAF+1), dm0(N_Knot,N_Knot), dm1(N_Knot,N_Knot), dm2(N_Knot,N_Knot), TSAVE(N_Knot,N_Knot), Amat(N_Knot,N_Knot), Bmat(N_Knot,N_Knot)
REAL(KIND=WP_Kind) TOE(0:nd_DAF), Val, dxinterp, dmat0(ninterp,N_Knot), dmat1(ninterp,N_Knot), dmat2(ninterp,N_Knot), vpot(N_Knot), Hpot(N_Knot), Mpot(N_Knot), MLRpot(N_Knot)
REAL(KIND=WP_Kind), PARAMETER::  UScale=1.0_WP_Kind, one=1.0_WP_Kind, two=2.0_WP_Kind, ten=10.0_WP_Kind, twelve=12.0_WP_Kind, usys2=6.0*1822.d0

twopi=8.d0*atan(1.d0)

OPEN(UNIT=4,File=TRIM(OutDIr)//"Test_DAF_And_Numerov_out.txt")
OPEN(UNIT=85,FILE=TRIM(OutDIr)//"DAF_out.txt")

!now we will solve the Schrodinger Equation

!Matrix Numerov
Potential_Name="Harmonic"
dx=(b-a)/(N_Knot-1)
DO i=1,N_Knot
   X_Knot(i)=a+(i-1)*dx
   HPot(i)=0.5d0*(x_Knot(i)-7.8775)**2
ENDDO
Call DAF_Matrix(N_Knot, dx, usys2, X_Knot, HPot, Eig, dm2, nd_daf, ndaf, DAFx, dm2, Potential_Name)
Call Numerov_Matrix(N_Knot, dx, usys2, X_Knot, HPot, Eig, dm2, dm0, dm1, Potential_Name)

!Morse Oscillators
Potential_Name="Morse"
DO i=1,N_knot
    r=X_Knot(i)
    CALL MorsePotential(r, Potential)
    MPot(i)=Potential
    dm2(i,i)=dm2(i,i)+Potential
ENDDO
Call DAF_Matrix(N_Knot, dx, usys2, X_Knot, MPot, Eig, dm2, nd_daf, ndaf, DAFx, dm2, Potential_Name)
Call Numerov_Matrix(N_Knot, dx, usys2, X_Knot, MPot, Eig, dm2, dm0, dm1, Potential_Name)

!Morse Long Range Potential (MLR)
Potential_Name="Morse Long Range"
OPEN(Unit=141,File=TRIM(OutDIr)//"Potentials.csv")
dm2=TSAVE
!CALL MLR_Potential(reau, Potential)
DO i=1,N_knot
    r=X_Knot(i)
    CALL MLR_Potential(r, Potential)
    vpot(i)=potential
    MLRPot(i)=potential
    WRITE(141,'(4(ES15.7,","))')X_Knot(i), MPot(i), MLRPot(i)
    dm2(i,i)=dm2(i,i)+Potential
ENDDO
Call DAF_Matrix(N_Knot, dx, usys2, X_Knot, MLRPot, Eig, dm2, nd_daf, ndaf, DAFx, dm2, Potential_Name)
Call Numerov_Matrix(N_Knot, dx, usys2, X_Knot, MLRPot, Eig, dm2, dm0, dm1, Potential_Name)

STOP "Execution Complete"
END