PROGRAM pendula
*     animation of linear and nonlinear pendula
      IMPLICIT NONE
      REAL*8 dt, dt_2, omega2, r, t, theta1, theta2, v1, v2
      INTEGER ball(2), IC1_, IC2_, IR_
      LOGICAL TBkeyinput
      DATA IC1_/1/
      DATA IC2_/2/
      CALL GWopen(IR_, 0)
      CALL initial(IC1_,IC2_,
     +             theta1,v1,theta2,v2,omega2,t,dt,dt_2,ball,r)
      DO WHILE(.NOT.TBkeyinput())
         CALL linear(theta1,v1,omega2,dt,dt_2)
         CALL nonlinear(theta2,v2,omega2,dt,dt_2)
         t = t + dt
         CALL animation(IC1_,theta1,ball(1))
         CALL animation(IC2_,theta2,ball(2))
      END DO
      CALL GWquit(IR_)
      END
C
      SUBROUTINE initial(IC1_,IC2_,
     +                   theta1,v1,theta2,v2,omega2,t,dt,dt_2,ball,r)
      IMPLICIT NONE
      REAL*8 theta1, v1, theta2, v2, omega2, t, dt, dt_2, r
      REAL xmax, xwin, ymin, ywin, GWaspect
      INTEGER ball(2), IR_, IC1_, IC2_
      t = 0
      theta1 = 0.25D0               ! position (radians) of linear oscillator
      theta2 = theta1               ! position of nonlinear oscillator
      v1 = 0
      v2 = v1
      omega2 = 9                            ! g/L ratio
      dt = 0.01D0
      dt_2 = 0.5D0*dt
      xmax = theta1                         ! dimension of window
      CALL GWvport(IR_, 0.0, 0.0, GWaspect(1), 0.5) ! bottom half of screen
      CALL compute_aspect_ratio(xmax,xwin,ywin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      CALL GWputtxt(IR_, -xwin*0.9, -ywin*0.9, 'linear oscillator')
*     draw circle
      CALL GWsetpen(IR_, 13, -1, -1, -1)
      r = 0.1D0                             ! radius of circle
      CALL GWcmbmrk(ball(1), 0, 6, REAL(2*r), 13, -1, 4)
      CALL GWsavevp(IR_, IC1_)
      CALL GWvport(IR_, 0.0, 0.5, GWaspect(1), 1.0) ! top half of screen
      CALL compute_aspect_ratio(xmax,xwin,ymin)
      CALL GWindow(IR_, -xwin, -ywin, xwin, ywin)
      CALL GWputtxt(IR_, -xwin*0.9, -ywin*0.9, 'nonlinear oscillator')
      CALL GWsetpen(IR_, 16, -1, -1, -1)
      CALL GWcmbmrk(ball(2), 0, 6, REAL(2*r), 16, -1, 4)
      CALL GWsavevp(IR_, IC2_)
      CALL GWanchor(IR_, 1)
      END
C
      SUBROUTINE linear(theta,v,omega2,dt,dt_2)
      IMPLICIT NONE
      REAL*8 a, amid, theta_mid, vmid, theta, v, omega2, dt, dt_2
      a = -omega2*theta
      vmid = v + a*dt_2
      theta_mid = theta + v*dt_2
      amid = -omega2*theta_mid
      v = v + amid*dt
      theta = theta + vmid*dt
      END
C
      SUBROUTINE nonlinear(theta,v,omega2,dt,dt_2)
      IMPLICIT NONE
      REAL*8 a, amid, theta_mid, vmid, theta, v, omega2, dt, dt_2
      a = -omega2*sin(theta)
      vmid = v + a*dt_2
      theta_mid = theta + v*dt_2
      amid = -omega2*sin(theta_mid)
      v = v + amid*dt
      theta = theta + vmid*dt
      END
C
      SUBROUTINE animation(IC9_,theta,ball)
      IMPLICIT NONE
      REAL*8 theta
      INTEGER ball, IR_, IC9_
      CALL GWselvp(IR_, IC9_)
      CALL GWsetogn(IR_, 0, -IC9_*10)
      CALL GWputcmb(IR_, ball, REAL(theta), 0.0, 0.0, 0.0, 0)
      CALL GWflush(IR_, -IC9_*10)
      END