!!!!!!!!!!!!!!!!!!!!!!!!!!!   Program 5.A   !!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!                                                                       !
! Please Note:                                                          !
!                                                                       !
! (1) This computer program is written by Tao Pang in conjunction with  !
!     his book, "An Introduction to Computational Physics," published   !
!     by Cambridge University Press in 1997.                            !
!                                                                       !
! (2) No warranties, express or implied, are made for this program.     !
!                                                                       !
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!
MODULE CB
  REAL :: Q,B,W
END MODULE CB
!
PROGRAM PENDULUM
!
! Program for the power spectra of a driven pendulum under damping with
! the fourth order Runge-Kutta algorithm. Given parameters: Q, B, and W
! (omega_0).  Copyright (c) Tao Pang 1997.
!
  USE CB
  IMPLICIT NONE
  INTEGER, PARAMETER :: N=65536,L=128,M=16,MD=16
  INTEGER :: I,J
  REAL :: PI,F1,H,OD,T,Y1,Y2,G1,GX1,G2,GX2
  REAL :: DK11,DK21,DK12,DK22,DK13,DK23,DK14,DK24
  REAL, DIMENSION (N) :: AR,AI,WR,WI,O
  REAL, DIMENSION (2,N) :: Y

!
  PI = 4.0*ATAN(1.0)
  F1 = 1.0/SQRT(FLOAT(N))
  W  = 2.0/3.0
  H  = 2.0*PI/(L*W)
  OD = 2.0*PI/(N*H*W*W)
  Q  = 0.5
  B  = 1.15
  Y(1,1) = 0.0
  Y(2,1) = 2.0
!
! Runge-Kutta algorithm to integrate the equation
!
  DO I = 1, N-1
    T  = H*I
    Y1 = Y(1,I)
    Y2 = Y(2,I)
    DK11 = H*GX1(Y1,Y2,T)
    DK21 = H*GX2(Y1,Y2,T)
    DK12 = H*GX1((Y1+DK11/2.0),(Y2+DK21/2.0),(T+H/2.0))
    DK22 = H*GX2((Y1+DK11/2.0),(Y2+DK21/2.0),(T+H/2.0))
    DK13 = H*GX1((Y1+DK12/2.0),(Y2+DK22/2.0),(T+H/2.0))
    DK23 = H*GX2((Y1+DK12/2.0),(Y2+DK22/2.0),(T+H/2.0))
    DK14 = H*GX1((Y1+DK13),(Y2+DK23),(T+H))
    DK24 = H*GX2((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
!
! Bring theta back to region [-pi,pi]
!
     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)
     END IF
  END DO
!
  DO I = 1, N
    AR(I) = Y(1,I)
    WR(I) = Y(2,I)
    AI(I) = 0.0
    WI(I) = 0.0
  END DO
  CALL FFT (AR,AI,N,M)
  CALL FFT (WR,WI,N,M)
!
  DO 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))
  END DO
  WRITE(6,"(3F16.10)") (O(I),AR(I),WR(I),I=1,(L*MD),4)
END PROGRAM PENDULUM
!
  SUBROUTINE FFT (AR,AI,N,M)
!
! An example of the fast Fourier transform subroutine with N = 2**M.
! AR and AI are the real and imaginary part of data in the input and
! corresponding Fourier coefficients in the output.
! Copyright (c) Tao Pang 1997.
!
  IMPLICIT NONE
  INTEGER, INTENT (IN) :: N,M
  INTEGER :: N1,N2,I,J,K,L,L1,L2
  REAL :: PI,A1,A2,Q,U,V
  REAL, INTENT (INOUT), DIMENSION (N) :: AR,AI
!
  PI = 4.0*ATAN(1.0)
  N2 = N/2
!
  N1 = 2**M
  IF(N1.NE.N) STOP 'Indices do not match'
!
! Rearrange the data to the bit reversed order
!
  L = 1
  DO 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
    END IF
    J   = N2
    DO WHILE (J.LT.L)
      L = L-J
      J = J/2
    END DO
    L = L+J
  END DO
!
! Perform additions at all levels with reordered data
!
  L2 = 1
  DO L = 1, M
    Q  =  0.0
    L1 =  L2
    L2 =  2*L1
    DO K = 1, L1
      U   =  COS(Q)
      V   = -SIN(Q)
      Q   =  Q + PI/L1
      DO 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
      END DO
    END DO
  END DO
END SUBROUTINE FFT
!
FUNCTION GX1 (Y1,Y2,T) RESULT (G1)
!
  G1 = Y2
END FUNCTION GX1
!
FUNCTION GX2 (Y1,Y2,T) RESULT (G2)
  USE CB
!
  G2 = -Q*Y2-SIN(Y1)+B*COS(W*T)
END FUNCTION GX2