ccccccccccccccccccccccccc     Program 5.A     ccccccccccccccccccccccccccc
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c                                                                       c
c Please Note:                                                          c
c                                                                       c
c (1) This computer program is written by Tao Pang in conjunction with  c
c     his book, "An Introduction to Computational Physics," published   c
c     by Cambridge University Press in 1997.                            c
c                                                                       c
c (2) No warranties, express or implied, are made for this program.     c
c                                                                       c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
      PROGRAM PENDULUM
C
C Program for the power spectra of a driven pendulum under
C damping with the fourth order Runge-Kutta algorithm. Given
C parameters: Q, B, and W (omega_0).
C Copyright (c) Tao Pang 1997.
C
      PARAMETER (N=65536, L=128, M=16, MD=16)
      DIMENSION Y(2,N),AR(N),AI(N),WR(N),WI(N),O(N)
      COMMON /CONST/ Q,B,W
      PI = 4.*ATAN(1.)
      F1 = 1.0/SQRT(FLOAT(N))
      W  = 2./3.
      H  = 2.*PI/(L*W)
      OD = 2*PI/(N*H*W*W)
      Q  = 0.5
      B  = 1.15
      Y(1,1) = 0.0
      Y(2,1) = 2.0
C
C 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.),(Y2+DK21/2.),(T+H/2.))
        DK22 = H*G2((Y1+DK11/2.),(Y2+DK21/2.),(T+H/2.))
        DK13 = H*G1((Y1+DK12/2.),(Y2+DK22/2.),(T+H/2.))
        DK23 = H*G2((Y1+DK12/2.),(Y2+DK22/2.),(T+H/2.))
        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.*(DK12+DK13)+DK14)/6.
        Y(2,I+1) = Y(2,I)+(DK21+2.*(DK22+DK23)+DK24)/6.
C
C Bring theta back to region [-pi,pi]
C
        IF(ABS(Y(1,I+1)).GT.PI) THEN
          Y(1,I+1) = Y(1,I+1) - 2.*PI*ABS(Y(1,I+1))/Y(1,I+1)
        ENDIF
  100 CONTINUE
C
      DO    200  I = 1, N
        AR(I) = Y(1,I)
        WR(I) = Y(2,I)
        AI(I) = 0.0
        WI(I) = 0.0
  200 CONTINUE
      CALL FFT(AR,AI,N,M)
      CALL FFT(WR,WI,N,M)

      DO    300 I = 1, N
        O(I)  =  (I-1)*OD
        AR(I) =  (F1*AR(I))**2 + (F1*AI(I))**2
        WR(I) =  (F1*WR(I))**2 + (F1*WI(I))**2
        AR(I) = ALOG10(AR(I))
        WR(I) = ALOG10(WR(I))
  300 CONTINUE
      WRITE(6,999) (O(I),AR(I),WR(I),I=1,(L*MD),4)
      STOP
  999 FORMAT (3F16.10)
      END
C
      SUBROUTINE FFT (AR,AI,N,M)
C
C An example of the fast Fourier transform subroutine with
C N = 2**M. AR and AI are the real and imaginary part of
C data in the input and corresponding Fourier coefficients
C in the output. 
C
      DIMENSION AR(N),AI(N)
C
      PI = 4.0*ATAN(1.0)
      N2 = N/2
C
      N1 = 2**M
      IF(N1.NE.N) STOP 'Indices do not match'
C
C Rearrange the data to the bit reversed order
C
      L = 1
      DO      150  K = 1, N-1
        IF(K.LT.L) THEN
          A1    = AR(L)
          A2    = AI(L)
          AR(L) = AR(K)
          AR(K) = A1
          AI(L) = AI(K)
          AI(K) = A2
        ENDIF
        J   = N2
        DO      100  WHILE (J.LT.L)
          L = L - J
          J = J/2
  100   END DO
        L = L + J
  150 CONTINUE
C
C Perform additions at all levels with reordered data
C
      L2 = 1
      DO      200  L = 1, M
        Q  =  0.0
        L1 =  L2
        L2 =  2*L1
        DO    190  K = 1, L1
          U   =  COS(Q)
          V   = -SIN(Q)
          Q   =  Q + PI/L1
          DO  180  J = K, N, L2
            I     =  J + L1
            A1    =  AR(I)*U - AI(I)*V
            A2    =  AR(I)*V + AI(I)*U
            AR(I) =  AR(J) - A1
            AR(J) =  AR(J) + A1
            AI(I) =  AI(J) - A2
            AI(J) =  AI(J) + A2
  180     CONTINUE
  190   CONTINUE
  200   CONTINUE
        RETURN
        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