MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i, j, n_show, n_show_read, steps
  INTEGER, PARAMETER :: planet2_dat = 10 !manage UNIT numbers as a resource
  REAL(dp) :: G_times_M, earth_x, earth_x_read, pluto_x, pluto_x_read, dx, earth_y, pluto_y, dy, r, r2, r3, dr, dr2, dr3,          &
    earth_v_x, pluto_v_x, a_x, earth_v_y, earth_v_y_read, pluto_v_y, pluto_v_y_read, a_y, ratio1, ratio2
  REAL(dp), PARAMETER :: PI = ATAN2(0.0_dp, -1.0_dp), RATIO_PLUTO_SUN = 7.0D-09, RATIO_EARTH_SUN = 3.0D-06, dt = 0.1_dp !time-step in years
  CHARACTER(*), PARAMETER :: my_format = '(4f8.4)'
END MODULE Shared

PROGRAM planet2
  USE Shared
  IMPLICIT NONE
  WRITE(6, '(a23)', advance='no') "Earth x-position [AU]: "
  READ *, earth_x_read
  earth_x = earth_x_read !x-position of Earth; units of AU
  earth_y = 0.0_dp
  earth_v_x = 0.0_dp
  WRITE(6, '(a28)', advance='no') "Earth y-velocity [AU/year]: "
  READ *, earth_v_y_read
  earth_v_y = earth_v_y_read !y-velocity of Earth; units of AU/year
  G_times_M = 4.0_dp * PI * PI !G x mass of Sun
  ratio1 = RATIO_PLUTO_SUN * G_times_M !in units of solar masses; Pluto: 1.303E+22 kg
  WRITE(6, '(a23)', advance='no') "Pluto x-position [AU]: "
  READ *, pluto_x_read
  pluto_x = pluto_x_read !x-position of Earth; units of AU
  pluto_y = 0.0_dp
  pluto_v_x = 0.0_dp
  WRITE(6, '(a28)', advance='no') "Pluto y-velocity [AU/year]: "
  READ *, pluto_v_y_read
  pluto_v_y = pluto_v_y_read !y-velocity of Pluto; units of AU/year
  ratio2 = -RATIO_EARTH_SUN * G_times_M !in units of solar masses; Earth: 5.97237E+24 kg
  OPEN(planet2_dat, file='planet2.dat', access='sequential', status='unknown', action='write')
  WRITE(planet2_dat, my_format) earth_x, earth_y, pluto_x, pluto_y
  WRITE(6, '(a25)', advance='no') "Simulation time [years]: "
  READ *, n_show_read
  n_show = n_show_read !simulation time in years
  steps = INT(FLOAT(n_show) / dt) !number of steps corresponding to N_SHOW years
  orbit: DO i = 1, steps
    dx = pluto_x - earth_x
    dy = pluto_y - earth_y
    dr2 = dx * dx + dy * dy
    dr = SQRT(dr2)
    dr3 = dr2 * dr
    r2 = earth_x * earth_x + earth_y * earth_y !Earth
    r = SQRT(r2)
    r3 = r2 * r
    a_x = -G_times_M * earth_x / r3 + ratio1 * dx / dr3
    a_y = -G_times_M * earth_y / r3 + ratio1 * dy / dr3
    earth_v_x = earth_v_x + a_x * dt
    earth_v_y = earth_v_y + a_y * dt
    earth_x = earth_x + earth_v_x * dt
    earth_y = earth_y + earth_v_y * dt
    WRITE(planet2_dat, my_format) earth_x, earth_y, pluto_x, pluto_y !store
    r2 = pluto_x * pluto_x + pluto_y * pluto_y !Pluto
    r = SQRT(r2)
    r3 = r2 * r
    a_x = -G_times_M * pluto_x / r3 + ratio2 * dx / dr3
    a_y = -G_times_M * pluto_y / r3 + ratio2 * dy / dr3
    pluto_v_x = pluto_v_x + a_x * dt
    pluto_v_y = pluto_v_y + a_y * dt
    pluto_x = pluto_x + pluto_v_x * dt
    pluto_y = pluto_y + pluto_v_y * dt
    WRITE(planet2_dat, my_format) earth_x, earth_y, pluto_x, pluto_y !store
  END DO orbit
  CLOSE(planet2_dat)
END PROGRAM planet2