PROGRAM PENDULUM ! -------------------------------------------------------------------------------------- ! ! Copyright 2006 Matthew Norton ! ! -------------------------------------------------------------------------------------- ! This program is free software; you can redistribute it and/or modify ! it under the terms of the GNU General Public License as published by ! the Free Software Foundation; either version 2 of the License, or ! (at your option) any later version. ! ! This program is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU General Public License for more details. ! ! You should have received a copy of the GNU General Public License ! along with this program; if not, write to the Free Software ! Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA ! ! NOTE: This is in FORTRAN 90 and you will need a FORTRAN 90 complier to run this code ! -------------------------------------------------------------------------------------- ! A Simulation of a Simple Pendulum ! ! This program simulates the motion of a simple pendulum. Using defined initial conditions. ! This program uses the arccos(-1) to obtain a value for Pi. ! ! -------------------------------------------------------------------------------------- ! ! PROGRAM HISTORY ! ! AUTHOR: Matthew S. Norton ! CREATION DATE: January 26, 2006 ! UPDATED: November 22, 2006 ! CREATED FOR: PHY 112 Lab ! UPDATED for use in Differential Equations project ! -------------------------------------------------------------------------------------- 1 IMPLICIT NONE ! Declare variables 2 INTEGER :: I, Icnt 3 Double Precision Pi, g, l, X, V, T, A, dt, Xold, Vold, Told, Xolder, Volder, dtold, CL2, TV, TVold 4 OPEN (30, FILE = "Euler.txt") ! Open the output file for Euler method 5 OPEN (40, FILE = "Euler-Cromer.txt") ! Open the output file for Euler-Cromer method ! Declare initial constants 6 Pi = ACOS(-1.) ! Claculate Pi from the arccos function, arccos(-1) = Pi 7 g = 9.8 ! Gravitational constant in m / s ** 2 8 l = 1. ! Length of pendulm in meters 9 T = 0 ! Start at time = zero 10 Told = 0 ! No negative time 11 CL2 = ((2.*Pi)/SQRT(g/l))*(1./1000.) 12 X = 0 ! Start at zero displacement 13 PRINT *, "Input initial angular displacement in degrees (it will be converted to radians)" 14 READ *, Xolder 15 Xold = Xolder * Pi / 180 !Initial displacement 16 V = 0 ! No inital velocity 17 PRINT *, "Input initial angular velocity (in radians per second)" 18 READ *, Volder 19 Vold = Volder 20 A = 0 21 dt = 0 22 WRITE (*,*) "Please input dt, the time increment" 23 WRITE (*, '(1X, A, 1X)', ADVANCE = "NO") "NOTE!! dt for CL-2 was 0.00200709! Enter 0 to use CL-2 value" 24 Read *, dtold 25 IF (dtold < 0.000001) THEN 26 dt = CL2 27 ELSE 28 dt = dtold 29 END IF 30 Icnt = 0 ! Euler Algorithim 31 DO I=1, 10000 32 Icnt = Icnt + 1 33 A = -g*sin(Xold)/L 34 V = Vold + A*dt 35 X = Xold +Vold*dt 36 T = Told +dt 37 IF (Icnt == 10) THEN 38 WRITE (30,1000) T, X, V, A 39 Icnt = 0 40 END IF 41 Vold = V 42 Xold = X 43 Told = T 44 END DO ! Reset variables 45 T = 0 ! Start at time = zero 46 Told = 0 ! No negative time 47 X = 0 ! Start at zero displacement 48 Xold = Xolder * Pi / 180 49 V = 0 ! No inital velocity 50 Vold = Volder 51 A = 0 52 Icnt = 0 ! Euler-Cromer Algorrithim 53 DO I=1, 10000 54 Icnt = Icnt + 1 55 A = -g*sin(Xold)/L 56 V = Vold + A*dt 57 X = Xold +V*dt 58 T = Told +dt 59 IF (Icnt == 10) THEN 60 WRITE (40,1000) T, X, V, A 61 Icnt = 0 62 END IF 63 Vold = V 64 Xold = X 65 Told = T 66 END DO 67 OPEN (50, FILE = "F-N.txt") ! Open the output file for Feynman-Newton method ! Reset variables 68 A = -( Xold + Vold) 69 Xold = Xolder * Pi / 180 70 Vold = Volder 71 Icnt = 0 72 X = 0 73 V = 0 74 A = 0 75 T = 0 76 TV = 0 77 Told = 0 78 Icnt = 0 79 TVold = 0 80 DO I = 1, 10000 81 Icnt = Icnt + 1 82 X = Xold + dt * V 83 A = -g*sin(X)/L 84 V = Vold + Dt * A 85 T = Told + Dt 86 TV = TVold + Dt / 2. 87 IF (Icnt == 10) THEN 88 WRITE (50, 1001) T, X, TV, V, A 89 Icnt = 0 90 END IF 91 Xold = X 92 Vold = V 93 Told = T 94 TVold = TV 95 END DO 1000 FORMAT (4F20.10) 1001 FORMAT (5F20.10) 96 CLOSE (30) ! Close the output file 97 CLOSE (40) ! Close the output file 98 CLOSE (50) ! Close the output file 99 PRINT *, "Output files, Euler.txt and Euler-Cromer.txt have been closed, program is complete!" 100 PRINT *, "Have a great day!" END PROGRAM