!========================================================
module          rk4_force
 implicit none
 private
 save
!--------------------------------------------------------
 integer, parameter :: NEQ = 4
!--------------------------------------------------------
 public             :: NEQ, k1,k2,k3,k4
 public             :: f  , energy
!--------------------------------------------------------
 integer, parameter :: dp  = kind(0.0_8)
 real(dp)           :: k1,k2,k3,k4
 real(dp)           :: r
!--------------------------------------------------------
!========================================================
contains
!========================================================
 function       f(t,X)               result(dXdt)
  real(dp) :: t,X(:)
  real(dp) :: dXdt(size(X))
  real(dp) :: x1,x2,v1,v2,a1,a2
!-----------------------
  x1      = X(1) ;  v1      = X(3) 
  x2      = X(2) ;  v2      = X(4)
!-----------------------
! Force field:
  r       = (x1*x1+x2*x2)**(-1.5_dp)
  a1      =  k1*x1*r
  a2      =  k1*x2*r
!-----------------------
  dXdt(1) = v1   ;  dXdt(3) = a1
  dXdt(2) = v2   ;  dXdt(4) = a2
!-----------------------
 end function   f
!--------------------------------------------------------
 function       energy(t,X)
  real(dp) :: energy
  real(dp) :: t,X(:)
  real(dp) :: x1,x2,v1,v2
  real(dp) :: r
!-----------------------
  x1      = X(1) ;  v1      = X(3) 
  x2      = X(2) ;  v2      = X(4)
!-----------------------
  r       = sqrt(x1*x1+x2*x2)
  energy  = 0.5_dp * (v1*v1+v2*v2) + k1/r
!-----------------------
 end function   energy
!--------------------------------------------------------
end module      rk4_force
!========================================================

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