MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i, j, nshow, nshow_read, steps
  INTEGER, PARAMETER :: planet2_dat = 10 !manage UNIT numbers as a resource
  REAL(dp) :: GM, earth_x, earth_x_read, pluto_x, pluto_x_read, dx, earth_y, pluto_y, dy, r, r2, r3, dr, dr2, dr3, earth_vx,       &
    pluto_vx, ax, earth_vy, earth_vy_read, pluto_vy, pluto_vy_read, ay, 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, FOUR = 4.0, dt = 0.1 !time-step in years
  CHARACTER(*), PARAMETER :: my_format = '(4f8.4)'
END MODULE Shared
!1. dp?, D+00?
!2. i + 1?
!3. Stupid things with _plot?
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_vx = 0.0_dp
  WRITE(6, '(a28)', advance='no') "Earth y-velocity [AU/year]: "
  READ *, earth_vy_read
  earth_vy = earth_vy_read !y-velocity of Earth; units of AU/year
  GM = FOUR * PI * PI !mass of Sun x G
  ratio1 = RATIO_PLUTO_SUN * GM !in units of solar masses; Pluto: 1.303E22 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_vx = 0.0_dp
  WRITE(6, '(a28)', advance='no') "Pluto y-velocity [AU/year]: "
  READ *, pluto_vy_read
  pluto_vy = pluto_vy_read !y-velocity of Pluto; units of AU/year
  ratio2 = -RATIO_EARTH_SUN * GM !in units of solar masses; Earth: 5.97237E24 kg
  OPEN(planet2_dat, file='planet2.dat', access='sequential', status='unknown', action='write', asynchronous='yes')
  WRITE(planet2_dat, my_format) earth_x, earth_y, pluto_x, pluto_y
  WRITE(6, '(a25)', advance='no') "Simulation time [years]: "
  READ *, nshow_read
  nshow = nshow_read !simulation time in years
  steps = INT(FLOAT(nshow) / dt) !number of steps corresponding to NSHOW 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
    ax = -GM * earth_x / r3 + ratio1 * dx / dr3
    ay = -GM * earth_y / r3 + ratio1 * dy / dr3
    earth_vx = earth_vx + ax * dt
    earth_vy = earth_vy + ay * dt
    earth_x = earth_x + earth_vx * dt
    earth_y = earth_y + earth_vy * 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
    ax = -GM * pluto_x / r3 + ratio2 * dx / dr3
    ay = -GM * pluto_y / r3 + ratio2 * dy / dr3
    pluto_vx = pluto_vx + ax * dt
    pluto_vy = pluto_vy + ay * dt
    pluto_x = pluto_x + pluto_vx * dt
    pluto_y = pluto_y + pluto_vy * dt
    WRITE(planet2_dat, my_format) earth_x, earth_y, pluto_x, pluto_y !store
  END DO orbit
  CLOSE(planet2_dat)
END PROGRAM planet2