!========================================================
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
!--------------------------------------------------------
!========================================================
contains
!========================================================
 function       f(t,X)               result(dXdt)
  real(dp) :: t,X(:)
  real(dp) :: dXdt(sizeof(X))
  real(dp) :: x1,x2,v1,v2,a1,a2
!-----------------------
  x1      = X(1) ;  v1      = X(3) 
  x2      = X(2) ;  v2      = X(4)
!-----------------------
! Force field:
  a1      =  -k2*v1
  a2      =  -k2*v2-k1
!-----------------------
  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
!-----------------------
  x1      = X(1) ;  v1      = X(3) 
  x2      = X(2) ;  v2      = X(4)
!-----------------------
  energy  = 0.5_dp * (v1*v1+v2*v2) + k1*x2
!-----------------------
 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/>.
!  -----------------------------------------------------------------------

