!========================================================
program         Runge_Kutta_4
 use rk4_force
 implicit none
 integer , parameter                 :: dp = kind(0.0_8)
 real(dp), allocatable               :: t(:),x(:,:)
 real(dp), dimension(NEQ)            :: x0
 real(dp)                            :: ti, tf
 integer                             :: Nt
 integer                             :: i , unit
!--------------------------------------------------------
!Input:
 print *,'# Runge-Kutta 4th Order Method for ODEs'
 print *,'# NEQ= ',NEQ
 print *,'# Enter k1, k2:'
 read  *,         k1, k2
 print *,'# ',    k1, k2
 print *,'# Enter Nt, ti, tf, x0: '
 read  *,         Nt, ti, tf, x0
 print *,'# ',    Nt, ti, tf, x0
!--------------------------------------------------------
 ALLOCATE(t(Nt),x(Nt,NEQ))
!--------------------------------------------------------
 call rk
!--------------------------------------------------------
 open(newunit=unit,file='rk4.dat')
 do i=1,Nt
  write(unit,'(2000G28.16)')                            &
       t(i),x(i,:),energy(t(i),x(i,:))
 enddo
 close (unit)
!========================================================
contains
!========================================================
!--------------------------------------------------------
 subroutine     rk
  real(dp) ::   dt, ts, xs(NEQ)
  integer  ::   i

  dt       = (tf-ti)/(Nt-1)
  t(1)     = ti
  x(1,:)   = x0
  ts       = ti
  xs       = x0

  do i     = 2, Nt
   call rkstep(ts,xs,dt)
   t(i)    =  ti  + (i-1) * dt
   x(i,:)  = xs
  end do
 end subroutine rk
!--------------------------------------------------------
 subroutine     rkstep(t,x,dt)
  real(dp), dimension(NEQ) :: x
  real(dp), dimension(NEQ) :: k1,k2,k3,k4
  real(dp)                 :: t ,dt
  real(dp)                 :: h ,h2,h6
!--------------------------------------------------------
  h  =  dt                 !h =dt, integration step
  h2 =  0.5_dp * h         !h2=h/2
  h6 =h/6.0_dp             !h6=h/6
!--------------------------------------------------------
  k1 = f(t   ,x      )
  k2 = f(t+h2,x+h2*k1)
  k3 = f(t+h2,x+h2*k2)
  k4 = f(t+h ,x+h *k3)
!--------------------------------------------------------
  t =t+h
  x =x +h6*(k1+2.0_dp*(k2+k3)+k4)
!--------------------------------------------------------
 end subroutine rkstep
!--------------------------------------------------------
end program     Runge_Kutta_4
!========================================================

!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2021)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------