    module nystrom

    use, intrinsic :: iso_fortran_env, wp=>real64

    implicit none

    private

    public :: RKN, RKNB, RKNC


    contains
    
    !------------------------------------------------------------------------------------------------------------------------------
    !
    !
    !
    !------------------------------------------------------------------------------------------------------------------------------

    subroutine RKNC(N, h, XStart, X, numb, t0, dt, dX, ddx,dxR,ModelA)
    use Base
    use Bridge
    implicit none
    
    TYPE (Model), TARGET :: ModelA
    
    integer N, numb, I, J, K, L
    REAL (KIND=dp) h,xstart,X(N),xf(N),xdf(N),t0,dt,dX(N),xddf(N),lx(numb),lxon(numb), load(numb),loadQ(numb), PLoad, QLoad,omega,KT(numb),K1,K2,lp, rdot(numb),adt, tau, ddx(numb),dxR(numb)

    write(*,200)
    write(srB,200)
    200 Format(//'-------------------------------------------------------------------------------------------',/,&
                 '             Algorithm :: Runge Kutta Nystrom Solver from Lear 1975 and Williams 2019')

    write(*,100)N
    write(srB,100)N
100 Format(//'             Number of equations :: ',i5,/)
    
    lp = 2.345
    K1 = 1000000.0
    K2 = 0.2*K1
    adt = 0.005
    ModelA%chi = ModelA%Pload/ModelA%G
    tau = pi*ModelA%c*t0/ModelA%length
    dxR = ZERO
    
    write(*,1101)ModelA%chi
1101 Format('             Chi Factor Load Ratios :: ',F7.3,//)    

    if(N .eq. 1) then
        write(*,500)
        write(sRB,500)
500     format('  Count  Time     XF     XDF')
    elseif(N .eq. 2) then
        write(*,510)
        write(sRB,510)
510     format('             Count  Time Velocity Distance Bridge   PLoad      Q load      K(x)    rdot   XF      XF     XDF     XDF',/,&
               '                                            eta                                                                     ')
    elseif(N .eq. 3) then
        write(*,520)
        write(sRB,520)
520     format('             Count  Time Tau T   Distance Bridge   Load 1      Load 2      K(x)    rdot          XF           XF          DY         XDF          XDF         DDY',/,&
               '                                            eta                                                                                                                   ',/)
    endif

    do 300 I = 1,numb
        lx(i) = ModelA%c * t0
        lxon(i) = lx(i)/ModelA%length
        
        if(lxon(i) .le. (ONE + 0.001)) then 
        load(i) = ONE * ModelA%PLoad
        loadQ(i) = ONE*ModelA%QLoad*sin(ModelA%omega*t0)
        KT(i) = K1 + K2*cos(TWO*pi*lx(i)/lp)
        rdot(i) = PointFIVE * adt*(ONE - cos(TWO*pi*lx(i)/lp))
        else
            load(i) = ZERO
            loadQ(i) = ZERO
            KT(i) = ZERO
            rdot(i) = ZERO
        endif
        
        !   subsitute tau for t0 as we have dimensionless form
        call nystrom_lear_44C(N,derivC,tau,X,DX,dt,xf,xdf,xddf, ModelA)
        ddx(i) = xddf(3)
        dxR(i) = xdf(3)
        
        
        
        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,ModelA%c,lx(i),lxon(i),load(i),loadQ(i),KT(i),rdot(i),(xf(k),k=1,N),(xdf(L),L=1,N)
            write(srB,410)I,t0,ModelA%c,lx(i),lxon(i),load(i),loadQ(i),KT(i),rdot(i),(xf(k),k=1,N),(xdf(L),L=1,N)
            write(srA,411)I,t0,ModelA%c,lx(i),lxon(i),load(i),loadQ(i),KT(i),rdot(i),(xf(k),k=1,N),(xdf(L),L=1,N)
        elseif(N .eq. 3) then
            write(*,412)I,t0,tau,lx(i),lxon(i),load(i),loadQ(i),KT(i),rdot(i),(xf(k),k=1,N),(xdf(L),L=1,N),ddx(i)
            write(*,*)dxR(i)
            write(srB,412)I,t0,tau,lx(i),lxon(i),load(i),loadQ(i),KT(i),rdot(i),(xf(k),k=1,N),(xdf(L),L=1,N),ddx(i)
            write(srA,414)I,t0,tau,lx(i),lxon(i),load(i),loadQ(i),KT(i),rdot(i),(xf(k),k=1,N),(xdf(L),L=1,N),ddx(i)
        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(12x,I4, 4(2x,F6.3),3(2x,F9.0),5(2x,F7.3))
411 Format(12x,I4, 4(',',F6.3),3(',',F9.0),5(',',F6.3))
412 Format(12x,I4, 4(2x,F6.3),3(2x,F9.0),3(2x,E10.3),2x,E10.3,3(2x,E10.3),2x,E10.3)
414 Format(12x,I4, 4(',',F6.3),3(',',F9.0),8(',',E10.3))
    t0 = t0 + dt
    
    tau = pi*ModelA%c*t0/ModelA%length
300 end do


    return
    end subroutine RKNC


    function derivC(n,t,x,xd, alpha, chi) 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),intent(in) :: alpha
    real(wp),intent(in) :: chi
    real(wp),dimension(n) :: xdd

    real(wp) R, A1, A2, A3,rt
    
    A1 = 1.0 / (alpha*alpha)
    A2 = 1.0/0.05
    A3 = 1.0 / alpha

    !R = SQRT(X(1)**2+X(2)**2)**3
    rt = 0.0
    xdd(1) = (A1*(X(1) - X(2))) +(A3*(XD(1) - XD(2)))
    xdd(2) = (1.0*A1*(1.05/(2.0*chi*0.05))) - (alpha*A2*(x(1) - x(2))) - (xd(1) - Xd(2))
    R = 10.0*(1.0 + 0.3*cos(20.0*t))*(x(1) - (x(3)*sin(t))-rt)
    if(R .lt. 0.0) then
        R = 1.0
       ! write(*,*)t,x(3),R,cos(20.0*t),x(1),(x(1) - (x(3)*sin(t))-rt)
        !write(*,*) 'R is zero'
    endif
    
    xdd(3) = A1*R*sin(t) - (A1* x(3)) - (A3* xd(3))

    end function derivC



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

    subroutine nystrom_lear_44C(n,df,t0,x0,xd0,dt,xf,xdf,xddf, ModelA)

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

    use Base
    use Bridge
    implicit none
    
    TYPE (Model), TARGET :: ModelA

    integer,intent(in) :: n !number of equations
    procedure (derivC) :: 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,xddf !final state

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

    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,ModelA%alpha,ModelA%chi)                       ! 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,ModelA%alpha,ModelA%chi)                       ! Df routine called
        dum = df(n,t,x,xd,ModelA%alpha,ModelA%chi)

        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,ModelA%alpha,ModelA%chi)                      ! 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,ModelA%alpha,ModelA%chi)                       ! Df routine called

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

    end if

    end subroutine nystrom_lear_44C


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

    subroutine RKNB(N, h, m, t0, dt, L, vt, H20, NT, Freq,NR, X,DX, XStart, ModelA)
    use Base
    use Bridge

    implicit none

    TYPE (Model), TARGET :: ModelA

    TYPE (EquationA), TARGET :: Equation1

    integer N
    integer I
    Integer m
    INTEGER NT
    Integer NR
    Integer K
    Integer J
    REAL (KIND=dp) h
    REAL (KIND=dp) X(N)
    REAL (KIND=dp) DX(N)
    REAL (KIND=dp) z(2)

    REAL (KIND=dp) xstart,xf(N),xdf(N), xdd(N)
    REAL (KIND=dp) t0
    REAL (KIND=dp) dt
    REAL (KIND=dp) time
    REAL (KIND=dp) L
    REAL (KIND=dp) vt                   !   Velocity of the truck
    REAL (KIND=dp) LAxle1
    REAL (KIND=dp) LAxle2
    REAL (KIND=dp) H20(NT)
    REAL (KIND=dp) eps(mn1), BGH

    REAL (KIND=dp) Ri(m,mn1)
    REAL(KIND=dp)Freq(NR,mn1)
    logical Axle1On
    logical Axle2On

    !-----------------------------------------------------------------------------------------------
    ! Start points
    z(1) = 0.0_DP
    z(2) = 0.0_DP
    BGH = 100000000000.0

    Write(*,100)N,h
    Write(sRB,100)N,h
100 Format(/,'----------------------------------------------------------',//,'   Runge Kutta Nystrom Solution for Bridges!',//,'   Number of equations :: ',I4,/,'   Tolerance (m) :: ',F10.5)

    if(N .ne. FIVE) then
        stop 'There are five equations refer to Fryba(1999) Vibrations of Solids and Structures undermoving loads'
    else
        write(*,350)
        write(sRB,350)
350     Format(/'   There are five equations refer to ::',/,'   Fryba, L., (1999) Vibrations of Solids and Structures undermoving loads, Telford: London UK',/)
    endif
    do 200 i=1,N
        if(i .eq. 1) then
            Write(*,300)
            Write(sRB,300)
300         Format('   Equation 1 :: Fryba Equation 36 ')
            call FirstEquation(ModelA, Equation1)
        else if(i .eq. 2) then
            Write(*,310)
            Write(SRB,310)
310         Format('   Equation 2 :: Fryba Equation 37 ')
        else if(i .eq. 3) then
            Write(*,320)
            Write(sRB,320)
320         Format('   Equation 3 :: Fryba Equation 38 ')
        else if(i .eq. 4) then
            Write(*,330)
            Write(sRB,330)
330         Format('   Equation 4 :: Fryba Equation 39 ')
        else if(i .eq. 5) then
            Write(*,340)
            Write(sRB,340)
340         Format('   Equation 5 :: Fryba Equation 52 ',//)
        endif
200 end do


    if(N .eq. 5) then
        write(*,500)
        500     format('  Count  Time                 Equation 1                                    Equation 2                                    Equation 3',/,&
                       '                     Phi         dPhi/dt        ddPhi/dtt        y3             dy3           ddy3         y1             dy1           ddy1',/)
    elseif(N .eq. 2) then
        write(*,510)
510     format('  Count  Time     XF      XF     XDF     XDF')
    endif
    time = t0
    LAxle1 = t0
    LAxle2 = LAxle1 - H20(1)

    do 600 i=1,m
        Ri(i,1) = ZERO
        Axle1On = .false.
        Axle2On = .false.
        if(LAxle1 .gt. ZERO .AND. LAxle1 .lt. L )then
            AXLE1On = .true.
            Ri(i,1) = ((ONE/FortyEight)*(pi**FOUR))*(H20(2)/gr)*((Freq(1,1)/Freq(1,1))**TWO)
        endif
        if(LAxle2 .gt. ZERO .AND. LAxle2 .lt. L )then
            AXLE2On = .true.
        endif
        if(LAxle1 .gt. L)then
            !AXLE1On = .false.
            Ri(i,1) = ZERO
        endif
        if(LAxle2 .gt. L)then
            ! AXLE2On = .false.
        endif

        if( Laxle1/L .le. ONE) then
            eps(1) = Laxle1/L
        else
            eps(1) = ZERO
        endif

        if( Laxle2/L .ge. ZERO .AND. Laxle2/L .le. ONE ) then
            eps(2) = Laxle2/L
        else
            eps(2) = ZERO
        endif


        call nystrom_lear_44B(N,bridgeF,t0,X,DX,dt,xf,xdf,z, Equation1, ModelA, xdd)

        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. 5) then
            write(*,420)I,time,xf(1)/BGH,xdf(1)/BGH,xdd(1)/BGH,xf(2)/BGH,xdf(2)/BGH,xdd(2)/BGH,xf(3)/BGH,xdf(3)/BGH,xdd(3)/BGH,xf(4)/BGH,xdf(4)/BGH,xdd(4)/BGH
            write(srB,410)I,time,(xf(k),k=1,N),(xdf(L),L=1,N)
            write(srA,411)I,time,xf(1)/BGH,xdf(1)/BGH,xdd(1)/BGH,xf(2)/BGH,xdf(2)/BGH,xdd(2)/BGH,xf(3)/BGH,xdf(3)/BGH,xdd(3)/BGH,xf(4)/BGH,xdf(4)/BGH,xdd(4)/BGH
        endif

        do 1310 J=1,N
            X(j) = XF(j)
            DX(j) = XDF(j)
1310 end do


    ! write(*,700)i,time,LAxle1,AXLE1On,LAxle2,AXLE2On,eps(1),eps(2),Ri(i,1)
    time = time + dt
    Laxle1 = Laxle1 + vt*dt
    Laxle2 = Laxle2 + vt*dt


600 end do

    return
700 Format(i4, 2x,f5.2,2x, f7.2,2x,L2,2x,f7.2,2x,L2,F7.3,2x,F7.3,F12.3)

400 Format(I5, 6(2x,F6.3))
401 Format(I5, 6(',',F6.3))
410 Format(I5, 11(2x,F8.3))
411 Format(I5, 13(',',E12.5))
420 Format(I5, 2x, F6.3,' |',3(2x,E12.5),' |',3(2x,E12.5),' | ',3(2x,E12.5),' |',3(2x,E12.5),' |')

    end subroutine RKNB

    !---------------------------------------------------------------------------------------------------------------------------------------------
    !
    !
    !
    !---------------------------------------------------------------------------------------------------------------------------------------------

    subroutine FirstEquation(ModelA, Equation1)
    use Bridge
    use Base

    implicit none

    TYPE (Model), TARGET :: ModelA

    TYPE (EquationA), TARGET :: Equation1



    Equation1%lambda = ModelA%H20(6)/(((ModelA%H20(2) + ModelA%H20(3))/gr)*ModelA%H20(1)**TWO)
    Equation1%alpha = ModelA%c/(TWO*ModelA%Freq(1,1)*ModelA%length)

    !
    !   lamdba from page 135 Fryba
    !
    write(*,100)Equation1%lambda,ModelA%Freq(1,1),Equation1%Alpha
    
100 format(//,'    Equation One Set Up ',//,'   Lambda :: ',F10.3,'   Normal Range [ 0.1 - 0.2]'/,'   Natural Frequency :: 'F10.3,'  Normal Range < 30 Hz.'/,'   Alpha  :: ',F10.3,'   Normal Range [0.02 - 0.1]'/)


    return
    end subroutine FirstEquation



    !-----------------------------------------------------------------------------------------------------------------------------------------------------------------------
    !
    !
    !
    !-----------------------------------------------------------------------------------------------------------------------------------------------------------------------

    function bridgeF(n,t,x,xd,z, Equation1, ModelA) result(xdd)

    use Bridge
    use Base
    !implicit none

    !import :: wp

    TYPE (Model), TARGET :: ModelA

    TYPE (EquationA), TARGET :: Equation1

    integer i
    integer,intent(in) :: n
    real(kind=dp),intent(in) :: t
    real(kind=dp),dimension(n),intent(in) :: x
    real(kind=dp),dimension(n),intent(in) :: xd
    real(kind=dp),dimension(2),intent(in) :: z
    real(kind=dp),dimension(n) :: xdd
    real(kind=dp) d(2), f(2),f3(2), zF(2),dzF(2), Q(2)

    real(kind=dp) R(2)

    d(1) = ModelA%H20(7)/ModelA%H20(1)
    d(2) = ModelA%H20(8)/ModelA%H20(1)
    f(1) = (ONE/(TWO*pi))*sqrt(ModelA%KC(1)/((ModelA%H20(2)+ModelA%H20(3))/gr))
    f(2) = (ONE/(TWO*pi))*sqrt(ModelA%KC(2)/((ModelA%H20(2)+ModelA%H20(3))/gr))
    f3(1) = (ONE/(TWO*pi))*sqrt(ModelA%KC(3)/((ModelA%H20(2)+ModelA%H20(3))/gr))
    f3(2) = (ONE/(TWO*pi))*sqrt(ModelA%KC(4)/((ModelA%H20(2)+ModelA%H20(3))/gr))
    Q(1) = 0.5
    Q(2) = 0.5
    R(1) = 10000.0
    R(2) = 10000.0
    
    do 50 i = 1,2
        zF(i) = ONE * x(2) - (-1**i)*d(i)*x(1) - x(i + 2)                        ! these are Fryba's equations 9.46 but the second one is a guess
        dzF(i) = 0.001 ! xd(2) - (-1**i)*d(i)*xd(1) - xd(i + 2)
50  end do
    
    xdd = ZERO
       
    do 100 i=1,2
        !equation 9.36
        ! Loop through front and rear axles front is 1 rear is 2
        ! Assume 9.29 is 0.5 
        xdd(1) = xdd(1) + ((((0.0_dp-1.0_dp)**i)*d(i))*(f3(i)/ModelA%Freq(1,1))*((pi*pi*(f3(i)/ModelA%Freq(1,1))*zF(i))+(Equation1%Alpha*PointFIVE*dzF(i))))
        
100 end do
    xdd(1) =   (ONE/(Equation1%lambda*Equation1%Alpha))*xdd(1)
   ! write(*,*)'here C :: ', xdd(1),zF(1),zF(2)
    xdd(2) = ZERO
    do 200 i = 1,2
        xdd(2) = xdd(2) + ((f3(i)/ModelA%Freq(1,1))*(pi*pi*(f3(i)/ModelA%Freq(1,1))*zf(i)+ 0.5*Equation1%Alpha*dzF(i)))
200 end do        
    
    xdd(2) = ZERO - (ONE/(0.9*Equation1%Alpha*Equation1%Alpha))*xdd(2)
    
    xdd(3) = (48.0/(pi*pi*0.1*0.05*Equation1%Alpha*Equation1%Alpha)) * ((TWO * Q(1) + ((pi*pi*pi*pi/48.0)*0.5*(f3(1)/ModelA%Freq(1,1))*(f3(1)/ModelA%Freq(1,1))*zf(1) - R(1))))
    xdd(4) = (48.0/(pi*pi*0.1*0.05*Equation1%Alpha*Equation1%Alpha)) * ((TWO * Q(2) + ((pi*pi*pi*pi/48.0)*0.5*(f3(2)/ModelA%Freq(1,1))*(f3(2)/ModelA%Freq(1,1))*zf(2) - R(2))))
    xdd(5) = 0.0

    end function BridgeF


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

    subroutine nystrom_lear_44B(n,df,t0,x0,xd0,dt,xf,xdf,z, Equation1, ModelA, xdd)

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

    implicit none

    TYPE (Model), TARGET :: ModelA

    TYPE (EquationA), TARGET :: Equation1
    integer,intent(in) :: n !number of equations
    procedure (bridgeF) :: 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),dimension(2),intent(in) :: z !variables in 9.36
    REAL (KIND=dp),intent(in) :: dt !time step
    REAL (KIND=dp),dimension(n),intent(out) :: xf,xdf,xdd !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,z, Equation1, ModelA)                       ! 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,z, Equation1, ModelA)                       ! 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,z, Equation1, ModelA)                       ! Df routine called

        t = t0+dt
        x = x0+dt*xd0+dt*(c1*k1+c3*k3)
        xd = xd0+cd1*k1+cd2*k2+cd3*k3
        xdd = df(n,t,x,xd,z, Equation1, ModelA)
        k4 = dt*xdd                                                      ! 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_44B

    !------------------------------------------------------------------------------------------------------------------------------
    !
    !
    !
    !------------------------------------------------------------------------------------------------------------------------------

    subroutine RKN(N, h, XStart, X, numb, t0, dt, dX)
    use Base
    implicit none
    integer N, numb, 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(*,100)N
    write(srB,100)N
100 Format('    Number of equations :: ',i5)

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

    do 300 I = 1,numb
        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


    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



    !------------------------------------------------------------------------------------------------------------------------------
    !
    ! 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


    end module nystrom