C LEAPINT.F: program to integrate hamiltonian system using leapfrog. program leapint parameter(MAXPNT = 100) ! maximum number of points integer n, mstep, nout, nstep real*8 x(MAXPNT), v(MAXPNT), tnow, dt C first, set up initial conditions n = 1 ! set number of points x(1) = 1.0 ! set initial position v(1) = 0.0 ! set initial velocity tnow = 0.0 ! set initial time C next, set integration parameters mstep = 256 ! number of steps to take nout = 4 ! steps between outputs dt = 1.0 / 32.0 ! timestep for integration C now, loop performing integration do 10 nstep = 0, mstep-1 ! loop mstep times in all if (mod(nstep, nout) .eq. 0) ! if time to output state \$ call printstate(x, v, n, tnow) ! then call output routine call leapstep(x, v, n, dt) ! take integration step tnow = tnow + dt ! and update value of time 10 continue if (mod(mstep, nout) .eq. 0) ! if last output wanted \$ call printstate(x, v, n, tnow) ! then output last step end C LEAPSTEP: take one step using the leap-from integrator, formulated C as a mapping from t to t + dt. WARNING: this integrator is not C accurate unless the timestep dt is fixed from one call to another. subroutine leapstep(x, v, n, dt) real*8 x(*) ! positions of all points real*8 v(*) ! velocities of all points integer n ! number of points real*8 dt ! timestep for integration parameter(MAXPNT = 100) integer i real*8 a(MAXPNT) call accel(a, x, n) ! call acceleration code do 10 i = 1, n ! loop over all points... v(i) = v(i) + 0.5 * dt * a(i) ! advance vel by half-step 10 continue do 20 i = 1, n ! loop over points again... x(i) = x(i) + dt * v(i) ! advance pos by full-step 20 continue call accel(a, x, n) ! call acceleration code do 30 i = 1, n ! loop over all points... v(i) = v(i) + 0.5 * dt * a(i) ! and complete vel. step 30 continue end C ACCEL: compute accelerations for harmonic oscillator(s). subroutine accel(a, x, n) real*8 a(*) ! accelerations of points real*8 x(*) ! positions of points integer n ! number of points integer i do 10 i = 1, n ! loop over all points... a(i) = - x(i) ! use linear force law 10 continue end C PRINTSTATE: output system state variables. subroutine printstate(x, v, n, tnow) real*8 x(*) ! positions of all points real*8 v(*) ! velocities of all points integer n ! number of points real*8 tnow ! current value of time integer i do 10 i = 1, n ! loop over all points... write(*, '(f8.4, i4, f12.6, f12.6)') tnow, i, x(i), v(i) 10 continue end