MODULE Shared
  USE, INTRINSIC :: ISO_Fortran_env, li=>INT64 !modern LONG INTEGER
  USE, INTRINSIC :: ISO_Fortran_env, dp=>REAL64 !modern DOUBLE PRECISION
  IMPLICIT NONE
  INTEGER :: i, i_max, n
  INTEGER(li) :: Factorial
  INTEGER, PARAMETER :: theta_dat = 10, omega_dat = 11, energy_dat = 12 !manage UNIT numbers as a resource
  REAL(dp) :: f_osc, theta_read, omega_read, L, L_read, g_over_L, elliptic, factor, period, dt, dt_read
  REAL(dp), DIMENSION(:), ALLOCATABLE, ASYNCHRONOUS :: theta, omega, a, E, t
  REAL(dp), PARAMETER :: g = 9.80665, PI = ATAN2(0.0_dp, -1.0_dp)
  REAL(dp), PARAMETER :: deg_to_rad = PI / 180.0_dp, rad_to_deg = 180.0_dp / PI
  CHARACTER(*), PARAMETER :: my_format = '(2f20.15)' !.15 is dp
  EXTERNAL Factorial
END MODULE Shared

PROGRAM amp
  USE Shared
  IMPLICIT NONE
  PRINT *, " "
  WRITE(6, '(a34)', advance='no') " Initial angle of the mass [deg]: "
  READ *, theta_read
  theta_read = theta_read * deg_to_rad !convert to radians
  n = 0
  elliptic = 0.0_dp
  sum_terms: DO
    IF (n > 5) EXIT !avoid long integers
    factor = GAMMA(FLOAT(2 * n + 1)) / (2**(2 * n)) / (GAMMA(FLOAT(n + 1))**2) !GAMMA(n + 1) == Factorial(n)
    elliptic = elliptic + (factor**2) * (SIN(theta_read / 2.0_dp)**(2 * n))
    n = n + 1
  END DO sum_terms
  WRITE(6, '(a29)', advance='no') " Length of the pendulum [m]: "
  READ *, L_read
  L = L_read
  g_over_L = g / L
  WRITE(6, '(a49)', advance='no') " Initial angular velocity of the mass [deg/sec]: "
  READ *, omega_read
  WRITE(6, '(a32)', advance='no') " Time-step of the system [sec]: "
  READ *, dt_read
  dt = dt_read
  PRINT *, " "
  period = 2.0_dp * PI * SQRT(1.0_dp / g_over_L) * elliptic
  WRITE(6, '(a25, f8.6, a4)') " Period of the pendulum: ", period, " sec"
  PRINT *, " "
  i_max = INT((2.0_dp * period) / dt)
  ALLOCATE(theta(i_max + 1), omega(i_max + 1), a(i_max + 1), t(i_max + 1), E(i_max + 1)) !keep "i_max + 1"
  theta(1) = theta_read !already converted
  omega(1) = omega_read * deg_to_rad !convert to rad/sec
  t(1) = 0.0_dp
  E(1) = -0.5_dp * omega(1)**2 + g_over_L * COS(theta(1))
  i = 1
  newton: DO
    IF (i > i_max) EXIT
    f_osc = -g_over_L * SIN(theta(i))
    a(i) = f_osc !N2L in "reduced units" of m = 1
    theta(i + 1) = theta(i) + omega(i) * dt + 0.5_dp * a(i) * dt * dt
    f_osc = -g_over_L * SIN(theta(i + 1))
    a(i + 1) = f_osc !N2L in "reduced units" of m = 1
    omega(i + 1) = omega(i) + 0.5_dp * (a(i + 1) + a(i)) * dt
    E(i + 1) = -0.5_dp * omega(i + 1) * omega(i + 1) + g_over_L * COS(theta(i + 1))
    t(i + 1) = t(i) + dt
    i = i + 1
  END DO newton
  DEALLOCATE(a)
  OPEN(theta_dat, file='theta.dat', access='sequential', status='unknown', action='write', asynchronous='yes')
  WRITE(theta_dat, my_format, asynchronous='yes') (t(i), theta(i) * rad_to_deg, i=1, i_max + 1) !convert to degrees
  CLOSE(theta_dat)
  DEALLOCATE(theta)
  OPEN(omega_dat, file='omega.dat', access='sequential', status='unknown', action='write', asynchronous='yes')
  WRITE(omega_dat, my_format, asynchronous='yes') (t(i), omega(i) * rad_to_deg, i=1, i_max + 1) !convert to deg/sec
  CLOSE(omega_dat)
  DEALLOCATE(omega)
  OPEN(energy_dat, file='energy.dat', access='sequential', status='unknown', action='write', asynchronous='yes')
  WRITE(energy_dat, my_format, asynchronous='yes') (t(i), E(i), i=1, i_max + 1)
  CLOSE(energy_dat)
  DEALLOCATE(t, E)
END PROGRAM amp

INTEGER(li) FUNCTION Factorial(nf)
  USE, INTRINSIC :: ISO_Fortran_env, li=>INT64 !modern LONG INTEGER
  IMPLICIT NONE
  INTEGER :: nf, i
  IF (nf < 0) ERROR STOP " Factorial is singular for negative integers."
  Factorial = 1
  product: DO i = 2, nf
    Factorial = Factorial * i
  END DO product
END FUNCTION Factorial