ccccccccccccccccccccccccc     Program 3.2     cccccccccccccccccccccccccc
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                      c
c Please Note:                                                         c
c                                                                      c
c (1) This computer program is part of the book, "An Introduction to   c
c     Computational Physics," written by Tao Pang and published and    c
c     copyrighted by Cambridge University Press in 1997.               c
c                                                                      c
c (2) No warranties, express or implied, are made for this program.    c
c                                                                      c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM PENDULUM
C
C Main Program for a driven pendulum under damping solved with
C the fourth-order Runge-Kutta algorithm.  Parameters: Q, B,
C and W (omega_0).
C
      PARAMETER (N=1000,L=100,M=1)
      DIMENSION Y(2,N)
      COMMON /CONST/ Q,B,W
      PI = 4.0*ATAN(1.0)
      H  = 3.0*PI/L
      Q  = 0.5
      B  = 0.9
      W  = 2.0/3.0
      Y(1,1) = 0.0
      Y(2,1) = 2.0
C
C Using the Runge-Kutta algorithm to integrate the equation
C
      DO      100  I = 1, N-1
        T  = H*I
        Y1 = Y(1,I)
        Y2 = Y(2,I)
        DK11 = H*G1(Y1,Y2,T)
        DK21 = H*G2(Y1,Y2,T)
        DK12 = H*G1((Y1+DK11/2.0),(Y2+DK21/2.0),(T+H/2.0))
        DK22 = H*G2((Y1+DK11/2.0),(Y2+DK21/2.0),(T+H/2.0))
        DK13 = H*G1((Y1+DK12/2.0),(Y2+DK22/2.0),(T+H/2.0))
        DK23 = H*G2((Y1+DK12/2.0),(Y2+DK22/2.0),(T+H/2.0))
        DK14 = H*G1((Y1+DK13),(Y2+DK23),(T+H))
        DK24 = H*G2((Y1+DK13),(Y2+DK23),(T+H))
        Y(1,I+1) = Y(1,I)+(DK11+2.0*(DK12+DK13)+DK14)/6.0
        Y(2,I+1) = Y(2,I)+(DK21+2.0*(DK22+DK23)+DK24)/6.0
C
C Bring theta back to the region [-pi,pi]
C
        Y(1,I+1) = Y(1,I+1)-2.0*PI*NINT(Y(1,I+1)/(2.0*PI))
  100 CONTINUE
      WRITE (6,999) (Y(1,I),Y(2,I),I=1,N,M)
      STOP
  999 FORMAT (2F16.8)
      END
C
      FUNCTION G1(Y1,Y2,T)
      COMMON /CONST/ Q,B,W
        G1 = Y2
      RETURN
      END
C
      FUNCTION G2(Y1,Y2,T)
      COMMON /CONST/ Q,B,W
        G2 = -Q*Y2-SIN(Y1)+B*COS(W*T)
      RETURN
      END