program kepler
  implicit none
  integer, parameter     :: neq = 4
  real, dimension(0:neq) :: trv
  real :: r_aphelion,eccentricity,dt,accuracy,a,T,vy0,G_m1_plus_m2,pi
  !
  logical :: step_using_y
  common /ystep/ step_using_y
  common /gcons/ pi, G_m1_plus_m2
  !
  step_using_y = .false.
  pi = 4 * atan(1.0);
  G_m1_plus_m2 = 4 * pi * pi  !  Units: G*M/(4*pi**2) = 1 -> T**2=a**3
  !
  print*, " Kepler orbit using fixed and then adaptive Runge-Kutta"

  print*, " Enter aphelion distance in AU: "
  read *, r_aphelion

  print*, " Enter eccentricity: "
  read *, eccentricity

  a = r_aphelion / (1 + eccentricity);
  T = a**1.5  ! 3rd Kepler law: T**2=a**3
  vy0 = sqrt(G_m1_plus_m2 * (2 / r_aphelion - 1 / a))

  print*, " Semimajor axis a = ", a, " AU "
  print*, " Period T = ",  T, " yr "
  print*, " v_y(0) = ", vy0, " AU/y "

  print*, " Enter step dt: "
  read *, dt

  print*, " Enter desired accuracy for adaptive integration: "
  read *, accuracy

  ! fixed step integration
  trv(0) = 0; trv(1) = r_aphelion; trv(2) = 0; trv(3) = 0; trv(4) = vy0;
  call integrate(trv, dt, T, accuracy, .false.);

  ! adaptive step integration
  trv(0) = 0; trv(1) = r_aphelion; trv(2) = 0; trv(3) = 0; trv(4) = vy0;
  call integrate(trv, dt, T, accuracy, .true.)

end program kepler


!  Derivative vector for Newton's law of gravitation
subroutine equations(trv,flow)
  integer, parameter   :: n=4
  real, dimension(0:n) :: trv,flow
  real                 :: t,x,y,vx,vy,r,ax,ay
  logical              :: step_using_y
  common /ystep/ step_using_y
  common /gcons/ pi, G_m1_plus_m2
  t  = trv(0)
  x  = trv(1)
  y  = trv(2)
  vx = trv(3)
  vy = trv(4)
  r  = sqrt(x*x + y*y)
  ax = - G_m1_plus_m2 * x / (r*r*r)
  ay = - G_m1_plus_m2 * y / (r*r*r)
  flow(0) = 1
  flow(1) = vx
  flow(2) = vy
  flow(3) = ax
  flow(4) = ay
  if (step_using_y) then   ! change independent variable from t to y
     do i = 0,n
        flow(i) = flow(i)/vy
     end do
  end if
  return
end subroutine equations


subroutine integrate(trv, dt, t_max, accuracy, adaptive)
  integer, parameter :: n = 4
  real, dimension(0:n) :: trv
  real         :: dt, dt_mmin, dt_max, t, t_max, accuracy, y_save
  logical      :: adaptive, step_using_y
  character*17 :: file_name
  character*10 :: f_or_a
  integer      :: step,iout
  common /ystep/ step_using_y
  common /gcons/ pi, G_m1_plus_m2
  external :: equations
  !
  file_name = "kepler_fixed.data"
  f_or_a="fixed"
  if (adaptive) then
     file_name = "kepler_adapt.data"
     f_or_a  = "adaptive"
  end if
  !
  open(newunit=iout,file=file_name,status="unknown");
  !
  print*, "Initial conditions:"
  call print_data(trv);
  print*, "Integrating with ", trim(f_or_a), " step size ..."
  step = 0;
  dt_min = dt
  dt_max = dt
  !  
  do while (.true.)
     call write_data(iout, trv)
     y_save = trv(2);
     step_using_y = .false.
     if (adaptive) then
        call RK4_adaptive_step(neq, trv, dt, equations, accuracy)
        dt_min = min(dt, dt_min);
        dt_max = max(dt, dt_max);
     else
        call RK4_step(neq, trv, dt, equations);
     end if
     t = trv(0); x = trv(1); y = trv(2); vx = trv(3); vy = trv(4)
     if (x > 0 .and. y * y_save < 0) then
        step_using_y = .true.
        call RK4_step(n, trv, -y, equations)
        call write_data(iout, trv)
        exit
     end if
     if (t > 10 * t_max) then
        write(*,*) "t too big, quitting ..."
        exit
     end if
     step = step + 1;
  end do
  !
  close(iout)
  !
  write(*,"(a,i3)") "Number of "//trim(f_or_a)//" steps = ", step
  if (adaptive) then
     write(*,*) "Minimum dt = ",dt_min
     write(*,*) "Maximum dt = ",dt_max
  end if
  call print_data(trv);
  write(*,*) "Trajectory data in: ", trim(file_name)
  !
end subroutine integrate


subroutine print_data(values)
  real, dimension(0:4) :: values
  character*4, dimension(0:4) :: names
  character*5, dimension(0:4) :: units
  names = ["t=  ", "x=  ", "y=  ", "v_x=", "v_y="]
  units = [" yr  ", " AU  ", " AU  ", "AU/yr", "AU/yr"]
  write(*,"(5(1x,a,f7.3,1x,a))") &
       (trim(names(i)), values(i), trim(units(i)),i=0,4)
end subroutine print_data


subroutine write_data(iout, values)
  real, dimension(0:4) :: values
  integer :: iout
  write(iout,"(5f8.4)") (values(i),i=0,4)
end subroutine write_data


