program rotating_particle
  implicit none

  integer, parameter  :: dp      = kind(0.0d0)
  real(dp), parameter :: pi      = acos(-1.0_dp)
  real(dp), parameter :: t_max   = 60 * pi
  real(dp), parameter :: dt      = pi / 10.0_dp
  integer, parameter  :: n_steps = nint(t_max/dt)
  real(dp), parameter :: x0(2)   = [0.0_dp, 1.0_dp]
  real(dp), parameter :: v0(2)   = [1.0_dp, 0.0_dp]

  integer             :: n, my_unit
  real(dp)            :: x(2), v(2), a(2)
  real(dp)            :: x_prev(2), v_prev(2)

  open(newunit=my_unit, file="forward_euler.txt")
  x = x0
  v = v0

  do n = 1, n_steps
     a = acceleration(x)
     x = x + v * dt
     v = v + a * dt
     write(my_unit, *) n, x, v, norm2(v), 0.5_dp*sum(v**2) - 1.0_dp/norm2(x)
  end do
  close(my_unit)

  open(newunit=my_unit, file="rk2.txt")
  x = x0
  v = v0

  do n = 1, n_steps
     x_prev = x
     v_prev = v

     a = acceleration(x)
     x = x + v * 0.5_dp * dt
     v = v + a * 0.5_dp * dt

     a = acceleration(x)
     x = x_prev + v * dt
     v = v_prev + a * dt
     write(my_unit, *) n, x, v, norm2(v), 0.5_dp*sum(v**2) - 1.0_dp/norm2(x)
  end do
  close(my_unit)

  open(newunit=my_unit, file="leap_frog.txt")
  x = x0
  v = v0
  do n = 1, n_steps
     x = x + v * 0.5_dp * dt
     v = v + acceleration(x) * dt
     x = x + v * 0.5_dp * dt
     write(my_unit, *) n, x, v, norm2(v), 0.5_dp*sum(v**2) - 1.0_dp/norm2(x)
  end do
  close(my_unit)

contains

  function acceleration(x) result(a)
    real(dp), intent(in) :: x(2)
    real(dp)             :: a(2)

    a = -x/norm2(x)**3
  end function acceleration

end program rotating_particle
