module integrator

  use parameters_module
 
  real (kind=8), dimension(:,:), allocatable, target:: xe
  real (kind=8), dimension(:,:), allocatable :: x, f
  real (kind=8), dimension(:), allocatable :: r22, r21, r2n
  real (kind=8), dimension(:), allocatable :: r1, r2, rn
  real (kind=8), dimension(:), allocatable :: a1, a2, a3

  real (kind=8) :: m1, m2, m3
  real (kind=8):: eps1, eps2

contains

!-----------------------------------------------------------------
!
subroutine init_rkvar(x0, mass1, mass2, epsilon1, epsilon2)
  implicit none
  
  real (kind=8), dimension(:,:), intent(in) :: x0
  real (kind=8), intent(in) :: mass1, mass2, epsilon1, epsilon2
  integer (kind=4) :: n

  n = size(x0, 2)

  allocate(x(6,n))
  allocate(f(6,n))
  allocate(xe(6,n))

  allocate(r22(n))
  allocate(r21(n))
  allocate(r2n(n))
  allocate(r1(n))
  allocate(r2(n))
  allocate(rn(n))
  allocate(a1(n))
  allocate(a2(n))
  allocate(a3(n))


  m1 = mass1
  m2 = mass2
  m3 = mass2
  eps1 = epsilon1*epsilon1
  eps2 = epsilon2*epsilon2
 

end subroutine init_rkvar

!-----------------------------------------------------------------
!
subroutine deallocate_rkvar

  deallocate(r22)
  deallocate(r21)
  deallocate(r2n)
  deallocate(r1)
  deallocate(r2)
  deallocate(rn)
  deallocate(a1)
  deallocate(a2)  
  deallocate(a3)
  deallocate(f)

  deallocate(x)
  deallocate(xe)

  return
end subroutine deallocate_rkvar


!-----------------------------------------------------------------
!
subroutine rk4(x0, h, xout)
  
  implicit none

  real (kind=8), dimension(:,:), intent(in) :: x0
  real (kind=8) :: h
  real (kind=8), dimension(:,:), intent(inout) :: xout

  integer (kind=4) ::  n

  n = size(x0,2)
  x = x0

  call diffeq(x)
  xe = x0 + h * f / 6.0d0
  x  = x0 + h * f / 2.0d0

  call diffeq(x)
  xe = xe + h * f / 3.0d0
  x  = x0 + h * f / 2.0d0

  call diffeq(x)
  xe = xe + h * f / 3.0d0
  x  = x0 + h * f

  call diffeq(x)
  xe = xe + h * f / 6.0d0

  xout = xe
  
return
end subroutine rk4
!
!
!----------------------------------------------------------------
subroutine diffeq(x)

  implicit none

  real (kind=8), dimension(:,:), intent(in) :: x  
  real (kind=8), dimension(6) ::xn
  integer (kind=4) :: n
  
  n = size(x,2)
  xn = x(:,n)

  r22 = (x(1,:)-xn(1))**2  +(x(2,:)-xn(2))**2 + (x(3,:)-xn(3))**2
  r21 = x(1,:)**2  + x(2,:)**2  + x(3,:)**2
  r2n = xn(1)**2 + xn(2)**2 + xn(3)**2
  
  r2 = sqrt(r22)
  r1 = sqrt(r21)
  rn = sqrt(r2n)

  a1 = -m1 / (r21 + eps1)
  a2 = -m2 / (r22 + eps2)
  a3 = -m3 / (r2n + eps2)

  ! this is a correction to prevent NaN errors in the vectorized
  ! function evalution at the location of the second mass
  r2(n) = 1.0d0

  ! calculate the RHS of the diffeq

  f(1,:) = x(4,:)
  f(2,:) = x(5,:)
  f(3,:) = x(6,:)

  f(4,:) = a1*x(1,:)/r1 + a2*(x(1,:)-xn(1))/r2 + a3*xn(1)/rn
  f(5,:) = a1*x(2,:)/r1 + a2*(x(2,:)-xn(2))/r2 + a3*xn(2)/rn
  f(6,:) = a1*x(3,:)/r1 + a2*(x(3,:)-xn(3))/r2 + a3*xn(3)/rn

  return
end subroutine diffeq

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



end module integrator