PROGRAM PENDULUM C C This program calculates the oscillation of a simple pendulum (small C angle oscillation) using RK4. C C23456789012345678901234567890123456789012345678901234567890123456789012 C EXTERNAL PENDEQN C INTEGER NEQN, IOSC, NOSC, NSTEPS, IFINAL, LFN, LOSCHOLD(1200) REAL*8 TIME(2000), OMEGA(1200), THETA(1200), TSTEP, RLENGTH REAL*8 T, THETA0, GRAV, PI REAL*8 PERIOD, X, Y(2), YOUT(2), DYDT(2) CHARACTER ERRMSG*80, OUTFILE*20 C C Pass GRAV and RLENGTH in a common block. C COMMON /PEND/ GRAV, RLENGTH C C Define constant values. C DATA GRAV, PI / 9.8036D0, 3.141592653589793D0 / DATA NEQN, NOSC, NSTEPS / 2, 5, 1200 / C C Set the time step. C TSTEP = 0.010D0 C C Read in the input values C PRINT *, 'Enter the length of the pendulum in meters:' READ *, RLENGTH PRINT *, 'Enter the initial offset angle in degrees:' READ *, THETA0 C C Define initial values and convert the angle into radians. C OUTFILE = ' ' THETA(1) = THETA0 * PI / 180.0D0 OMEGA(1) = 0.0D0 LOSCHOLD(1) = IOSC DO I = 1, 80 ERRMSG(I:I) = ' ' ENDDO C C Read in outout data file name. C LFN = 0 PRINT *, 'Enter the rootname of the output file: ' READ *, OUTFILE DO I = 1, 20 IF (LFN .EQ. 0) THEN IF (OUTFILE(I:I) .EQ. ' ') LFN = I - 1 IF (OUTFILE(I:I) .EQ. '.') LFN = I - 1 ENDIF ENDDO PRINT *, 'Output will be stored in file: ', OUTFILE(1:LFN), '.txt' OPEN(12, FILE=OUTFILE(1:LFN) // '.txt', STATUS='NEW') C C Create the time array. C DO I = 1, NSTEPS TIME(I) = DFLOAT(I-1) * TSTEP ENDDO C C Calculate the period of oscillation. C PERIOD = 2.0D0 * PI * DSQRT(RLENGTH / GRAV) C C Note that our 2 differential equations are C d-theta/dt = omega = DYDT(1) C d-omega/dt = -(grav/rlength) * theta = DYDT(2) C Set up the initial values for DYDT and for Y: C Y(1) = theta(i) C Y(2) = omega(i) C C Set the initial values of X and Y. C T = TIME(1) Y(1) = THETA(1) Y(2) = OMEGA(1) C C Calculate the initial values of DYDT. C CALL PENDEQN(T, Y, DYDT, NEQN, ERRMSG) IF (ERRMSG(1:1) .NE. ' ') GOTO 999 C C Initialize counters and switches. C IOSC = 1 IFINAL = 1 LOSCHOLD(1) = IOSC C C Integrate through the NOSC oscillations. C DO I = 2, NSTEPS T = TIME(I) C C Check to see if NOSC oscillations have been completed. C IF (T .GE. PERIOD*DFLOAT(IOSC)) IOSC = IOSC + 1 LOSCHOLD(I) = IOSC IF (IOSC .GT. NOSC) THEN IFINAL = I - 1 GOTO 900 ENDIF C CALL RK4_PENDULUM(Y, DYDT, NEQN, T, TSTEP, YOUT, PENDEQN, + ERRMSG) IF (ERRMSG(1:1) .NE. ' ') GOTO 999 THETA(I) = YOUT(1) OMEGA(I) = YOUT(2) Y(1) = YOUT(1) Y(2) = YOUT(2) ENDDO C C Output the data to a file. C 900 PRINT *, 'Storing data to output file.' WRITE (12, 905) RLENGTH, THETA0, PERIOD, TSTEP, GRAV 905 FORMAT('Pendulum motion for RLENGTH = ', F6.4, ' m, THETA0 = ', + F6.3, ' deg', /, 'PERIOD = ', F6.3, ' s, TSTEP = ', F6.3, + ' s, GRAVITY = ', F6.4, /, 'STEP', 3X, 'TIME (s)', 2X, + 'Osc. Num.', 2X, 'THETA (deg)', 2x, 'OMEGA (rad/s)') DO I = 1, IFINAL WRITE(12, 910) I, TIME(I), LOSCHOLD(I), THETA(I)*180.0D0/PI, + OMEGA(I) 910 FORMAT(I4, 2X, F8.4, 7X, I1, 6X, F8.3, 5X, 1PE12.4) ENDDO C C Close the output file. C CLOSE (12, STATUS='KEEP') PRINT *, 'Program PENDULUM done.' STOP C C Error handling portion of the code. C 999 PRINT *, 'Fatal error has occured in PENDULUM:' PRINT *, ' ', ERRMSG C STOP END