    module nystrom

    use, intrinsic :: iso_fortran_env, wp=>real64

    implicit none

    private

    ! derivative function prototype
    !abstract interface
    ! function deriv(n,t,x,xd) result(xdd)
    !  import :: wp
    !  integer,intent(in) :: n
    !  real(wp),intent(in) :: t
    !  real(wp),dimension(n),intent(in) :: x
    !  real(wp),dimension(n),intent(in) :: xd
    !  real(wp),dimension(n) :: xdd
    !  end function deriv
    !end interface
    public :: deriv
    public :: nystrom_lear_44
    public :: RKN


    contains

    !------------------------------------------------------------------------------------------------
    !
    ! Runge Kutta Nystrom Routine
    !
    !------------------------------------------------------------------------------------------------

    subroutine RKN(N, h, XStart, X, num, t0, dt, dX)
    use Base
    implicit none
    integer N, num, I, J, K, L
    REAL (KIND=dp) h,xstart,X(N),xf(N),xdf(N),t0,dt,dX(N)

    write(*,200)
    write(srB,200)
200 Format('  Algorithm :: Runge Kutta Nystrom Solver from Lear 1975 and Williams 2019')
    write(srB,100)N
100 Format('    Number of equations :: ',i5)

    if(N .eq. 1) then
        write(*,500)
500     format('  Count  Time     XF     XDF')
    elseif(N .eq. 2) then
        write(*,510)
510     format('  Count  Time     XF      XF     XDF     XDF')
    endif

    do 300 I = 1,num
        call nystrom_lear_44(N,deriv,t0,X,DX,dt,xf,xdf)
         if(N .eq. 1) then
        write(*,400)I,t0,(xf(k),k=1,N),(xdf(L),L=1,N)
        write(srB,400)I,t0,(xf(k),k=1,N),(xdf(L),L=1,N)
        write(srA,401)I,t0,(xf(k),k=1,N),(xdf(L),L=1,N)
         elseif(N .eq. 2) then 
             write(*,410)I,t0,(xf(k),k=1,N),(xdf(L),L=1,N)                  
             write(srB,410)I,t0,(xf(k),k=1,N),(xdf(L),L=1,N)      
             write(srA,411)I,t0,(xf(k),k=1,N),(xdf(L),L=1,N)
         endif
         
        do 310 J=1,N
            X(j) = XF(j)
            DX(j) = XDF(j)
310 end do
    
400 Format(I5, 3(2x,F6.3))
401 Format(I5, 3(',',F6.3))
410 Format(I5, 5(2x,F6.3))
411 Format(I5, 5(',',F6.3))
    t0 = t0 + dt
300 end do


    return
    end

    !------------------------------------------------------------------------------------------------------------------------------
    !
    ! 4th order Nystrom integration routine
    ! Take a step from t0 to t0+dt (x0,xd0 -> xf,xdf)
    !
    !------------------------------------------------------------------------------------------------------------------------------

    subroutine nystrom_lear_44(n,df,t0,x0,xd0,dt,xf,xdf)

    ! 4th order Nystrom integration routine
    ! Take a step from t0 to t0+dt (x0,xd0 -> xf,xdf)

    use Base

    implicit none

    integer,intent(in) :: n !number of equations
    procedure (deriv) :: df !derivative routine
    REAL (KIND=dp),intent(in) :: t0 !initial time
    REAL (KIND=dp),dimension(n),intent(in) :: x0,xd0 !initial state
    REAL (KIND=dp),intent(in) :: dt !time step
    REAL (KIND=dp),dimension(n),intent(out) :: xf,xdf !final state

    REAL (KIND=dp) :: t
    REAL (KIND=dp),dimension(n) :: x,xd,k1,k2,k3,k4

    REAL (KIND=dp),parameter :: s5 = sqrt(5.0_dp)
    REAL (KIND=dp),parameter :: delta2 = (5.0_dp-s5)/10.0_dp
    REAL (KIND=dp),parameter :: delta3 = (5.0_dp+s5)/10.0_dp
    REAL (KIND=dp),parameter :: delta4 = 1.0_dp
    REAL (KIND=dp),parameter :: a1 = (3.0_dp-s5)/20.0_dp
    REAL (KIND=dp),parameter :: b1 = 0.0_dp
    REAL (KIND=dp),parameter :: b2 = (3.0_dp+s5)/20.0_dp
    REAL (KIND=dp),parameter :: c1 = (-1.0_dp+s5)/4.0_dp
    REAL (KIND=dp),parameter :: c2 = 0.0_dp
    REAL (KIND=dp),parameter :: c3 = (3.0_dp-s5)/4.0_dp
    REAL (KIND=dp),parameter :: alpha1 = 1.0_dp/12.0_dp
    REAL (KIND=dp),parameter :: alpha2 = (5.0_dp+s5)/24.0_dp
    REAL (KIND=dp),parameter :: alpha3 = (5.0_dp-s5)/24.0_dp
    REAL (KIND=dp),parameter :: alpha4 = 0.0_dp
    REAL (KIND=dp),parameter :: ad1 = (5.0_dp-s5)/10.0_dp
    REAL (KIND=dp),parameter :: bd1 = -(5.0_dp+3.0_dp*s5)/20.0_dp
    REAL (KIND=dp),parameter :: bd2 = (3.0_dp+s5)/4.0_dp
    REAL (KIND=dp),parameter :: cd1 = (-1.0_dp+5.0_dp*s5)/4.0_dp
    REAL (KIND=dp),parameter :: cd2 = -(5.0_dp+3.0_dp*s5)/4.0_dp
    REAL (KIND=dp),parameter :: cd3 = (5.0_dp-s5)/2.0_dp
    REAL (KIND=dp),parameter :: beta1 = 1.0_dp/12.0_dp
    REAL (KIND=dp),parameter :: beta2 = 5.0_dp/12.0_dp
    REAL (KIND=dp),parameter :: beta3 = 5.0_dp/12.0_dp
    REAL (KIND=dp),parameter :: beta4 = 1.0_dp/12.0_dp

    if (dt==0.0_dp) then
        xf = x0
        xdf = xd0
    else

        t = t0
        x = x0
        xd = xd0
        k1 = dt*df(n,t,x,xd)                       ! Df routine called

        t = t0+delta2*dt
        x = x0+delta2*dt*xd0+a1*dt*k1
        xd = xd0+ad1*k1
        k2 = dt*df(n,t,x,xd)                       ! Df routine called

        t = t0+delta3*dt
        x = x0+delta3*dt*xd0+dt*(b1*k1+b2*k2)
        xd = xd0+bd1*k1+bd2*k2
        k3 = dt*df(n,t,x,xd)                       ! Df routine called

        t = t0+dt
        x = x0+dt*xd0+dt*(c1*k1+c3*k3)
        xd = xd0+cd1*k1+cd2*k2+cd3*k3
        k4 = dt*df(n,t,x,xd)                       ! Df routine called

        xf = x0+dt*xd0+dt*(alpha1*k1+alpha2*k2+alpha3*k3)
        xdf = xd0+beta1*k1+beta2*k2+beta3*k3+beta4*k4

    end if

    end subroutine nystrom_lear_44



    function deriv(n,t,x,xd) result(xdd)
    import :: wp
    integer,intent(in) :: n
    real(wp),intent(in) :: t
    real(wp),dimension(n),intent(in) :: x
    real(wp),dimension(n),intent(in) :: xd
    real(wp),dimension(n) :: xdd

    real(wp) R

    R = SQRT(X(1)**2+X(2)**2)**3
    xdd(1) = -X(1)/R
    xdd(2) = -X(2)/R

    end function deriv


    end module nystrom