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