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 :: theta, omega, a, E, t
  REAL(dp), PARAMETER :: g = 9.80665_dp, 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 = '(2f21.15)' !.15 is dp
  EXTERNAL Factorial
END MODULE Shared

PROGRAM amp
  USE Shared
  IMPLICIT NONE
  PRINT *, " "
  WRITE(6, '(a35)', advance='no') "  Initial angle of the mass [deg]: "
  READ *, theta_read !degrees
  IF (theta_read > 90.0_dp .OR. theta_read < -90.0_dp) STOP
  theta_read = theta_read * deg_to_rad !convert to radians
  WRITE(6, '(a30)', advance='no') "  Length of the pendulum [m]: "
  READ *, L_read
  L = L_read
  IF (L > 2.0_dp .OR. L < 0.2_dp) STOP
  g_over_L = g / L
  WRITE(6, '(a50)', advance='no') "  Initial angular velocity of the mass [deg/sec]: "
  READ *, omega_read !deg/sec
  IF (omega_read > 17.0_dp .OR. omega_read < -17.0_dp) STOP
  WRITE(6, '(a33)', advance='no') "  Time-step of the system [sec]: "
  READ *, dt_read
  dt = dt_read
  IF (dt > 0.100_dp .OR. dt < 0.001_dp) STOP
  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) !theta in radians
    n = n + 1
  END DO sum_terms
  period = 2.0_dp * PI * SQRT(1.0_dp / g_over_L) * elliptic
  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
  omega(1) = omega(1) * deg_to_rad; !convert to rad/sec
  a(1) = -g_over_L * SIN(theta(1))
  t(1) = 0.0_dp
  E(1) = -0.5_dp * omega(1) * omega(1) + g_over_L * COS(theta(1)) !equation of motion
  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)) !equation of motion
    t(i + 1) = t(i) + dt
    i = i + 1
  END DO newton
  DEALLOCATE(a)
  PRINT *, " "
  WRITE(6, '(a26, f9.7, a4)') "  Period of the pendulum: ", period, " sec"
  PRINT *, " "
  theta = theta * rad_to_deg !convert array to degrees
  OPEN(theta_dat, file='theta.dat', access='sequential', status='unknown', action='write')
  WRITE(theta_dat, my_format) (t(i), theta(i), i=1, i_max + 1)
  CLOSE(theta_dat)
  DEALLOCATE(theta)
  omega = omega * rad_to_deg !convert array to deg/sec
  OPEN(omega_dat, file='omega.dat', access='sequential', status='unknown', action='write')
  WRITE(omega_dat, my_format) (t(i), omega(i), i=1, i_max + 1)
  CLOSE(omega_dat)
  DEALLOCATE(omega)
  OPEN(energy_dat, file='energy.dat', access='sequential', status='unknown', action='write')
  WRITE(energy_dat, my_format) (t(i), E(i), i=1, i_max + 1)
  CLOSE(energy_dat)
  DEALLOCATE(t, E)
END PROGRAM amp

INTEGER(li) FUNCTION Factorial(nf) !not used
  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