PROGRAM MAIN IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PURPOSE: TO COMPUTE A POINCARE SECTION FOR A CHAOTIC PENDULUM. C CONCEPT PRESENTED AS BASIC PROGRAM FROM BAKER&GOLLUB, PG 154 C THIS F77 VERSION WRITTEN 1993 WHILE A STUDENT AT UW-RF. C COPYRIGHT 1993, 2007 JIM HALL, JHALL@FREEDOS.ORG C Permission is hereby granted, free of charge, to any person obtaining a copy C of this software and associated documentation files (the "Software"), to deal C in the Software without restriction, including without limitation the rights C to use, copy, modify, merge, publish, distribute, sublicense, and/or sell C copies of the Software, and to permit persons to whom the Software is C furnished to do so, subject to the following conditions: C C The above copyright notice and this permission notice shall be included in C all copies or substantial portions of the Software. C C THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR C IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, C FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE C AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER C LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, C OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN C THE SOFTWARE. PARAMETER (MSIZE = 60, NFILE = 20, PI = 3.14159265) CHARACTER FNAME*(MSIZE) C PRINT THE NAME OF THE PROGRAM PRINT *, 'THIS PROGRAM WILL COMPUTE A POINCARE SECTION' PRINT *, 'THE FOLLOWING ARE GOOD VALUES TO USE:' PRINT 1010, 'FORCE', 1.15 PRINT 1010, 'INVERSE-DAMPING', 2. PRINT 1010, 'ANGLE', 1. PRINT 1010, 'VELOCITY', 0. PRINT 1010, 'TIME MINIMUM', 2. PRINT 1010, 'TIME MAXIMUM', 10000. PRINT 1010, 'PHASE', 2. PRINT 1010, 'FREQUENCY', .66666 PRINT *, 'ENTER THE DRIVING FORCE:' READ *, G PRINT *, 'ENTER THE INVERSE-DAMPING (FOR NO DAMPING, USE 9999):' READ *, Q PRINT *, 'ENTER THE INITIAL ANGLE AND INITIAL VELOCITY:' READ *, XINT, VINT PRINT *, 'ENTER THE MINIMUM AND MAXIMUM TIMES:' READ *, TMIN, TMAX PRINT *, 'ENTER PHASE IN RADIANS:' READ *, PHI C ALLOW FOR SAVING TO A FILE PRINT *, 'ENTER THE FILE TO SAVE THE OUTPUT:' READ 1000, FNAME 1000 FORMAT (A) OPEN (UNIT=NFILE, FILE=FNAME, FORM='FORMATTED', STATUS='NEW', $ IOSTAT=IE) IF (IE.NE.0) THEN PRINT *, 'OOPS. CANNOT OPEN THAT FILE. QUITTING..' STOP ENDIF C SET SOME INITIAL VALUES W = 2. / 3. EPS = 1.E-5 TSTEP = .5 T = 0. X = XINT V = VINT C BEGIN THE ITERATIONS PRINT *, 'WORKING..' DO 100 I = 1, 10000 CALL RUNGEK (X, V, TSTEP, XNEW, VNEW, T, W, G, Q) TSHALF = TSTEP / 2. C TAKE TWO HALF STEPS CALL RUNGEK (X, V, TSHALF, XNH, VNH, T, W, G, Q) CALL RUNGEK (XNH, VNH, TSHALF, XN, VN, T + TSHALF, W, G, Q) D1 = ABS (XN - XNEW) D2 = ABS (VN - VNEW) DELTA = MAX (D1, D2) C CHECK TO SEE IF DELTA IS LESS THAN EPSILON, AND IF T IS BIGGER C THAN TMIN IF (DELTA.LT.EPS) THEN IF (T.GT.TMIN) THEN TNEW = T + STEP C CHECK IF THIS IS A NEW POINCARE SECTION PHASE1 = INT ((W * T - PHI) / (2. * PI)) PHASE2 = INT ((W * TNEW - PHI) / (2. * PI)) IF (PHASE2.GT.PHASE1) THEN TS = ((2. * PI * PHASE2 + PHI) / 2) - T C COMPUTE THE POINCARE SECTION CALL RUNGEK (X, V, TS, XP, VP, T, W, G, Q) IF (ABS(XP).GT.PI) THEN XP = XP - 2. * PI * ABS(XP) / XP ENDIF PRINT *, 'POINCARE!' WRITE (NFILE,*,IOSTAT=IE) XP, VP IF (IE.NE.0) THEN PRINT *, 'OOPS. CAN''T WRITE TO FILE!' CLOSE (NFILE) STOP ENDIF ENDIF ENDIF C ADVANCE THE VARIABLES.. IF NOT CHANGING MUCH, MAKE TSTEP BIGGER.. C IF CHANGING A LOT, MAKE TSTEP SMALLER X = XNEW V = VNEW T = T + TSTEP TSTEP = TSTEP * .95 * ((EPS / DELTA)**.25) C KEEP THETA IN BOUNDS IF (ABS(X).GT.PI) THEN X = X - 2. * PI * ABS(X) / X ENDIF ELSE TSTEP = TSTEP * .95 * ((EPS / DELTA)**.2) ENDIF IF (T.GT.TMAX) THEN GOTO 110 ENDIF 100 CONTINUE 110 CLOSE (NFILE) PRINT *, 'DONE!' PRINT 1010, 'G', G PRINT 1010, 'Q', Q PRINT 1010, 'PHI', PHI 1010 FORMAT (X, A, '=', T40, G13.5) STOP END CCCCCCCCCCCCCCCCCCCCCCC SUBROUTINE RUNGEK (X, V, TSTEP, XNEW, VNEW, T, W, G, Q) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C PERFORMS A RUNGE KUTTA SOLUTION OVER XNEW AND VNEW C C INPUT: C X -- THE INITIAL ANGLE, IN RADIANS C C V -- THE INITIAL SPEED, IN RADIANS PER SECOND C C TSTEP -- THE SIZE OF A STEP C C T -- THE TOTAL TIME, IN SECONDS C C W -- THE ANGULAR FREQUENCY C C G -- THE DRIVING FORCE C C Q -- THE INVERSE-DAMPING (FOR NO DAMPING, USE 9999) C C OUTPUT: C XNEW -- THE NEW ANGLE, IN RADIANS C C VNEW -- THE NEW SPEED, IN RADIANS PER SECOND * PRINT *, 'RUNGE-KUTTA' C THE FIRST RUNGE KUTTA TERM, ALSO THE EULER TERM XK1 = TSTEP * V VK1 = TSTEP * ACCEL (X, V, T, W, G, Q) C THE SECOND AND THIRD RUNGE KUTTA TERMS, AT THE MIDPOINT XK2 = TSTEP * (V + VK1 / 2.) VK2 = TSTEP * ACCEL (X + XK1 / 2., V + VK1 / 2., T + TSTEP / 2., W $ , G, Q) XK3 = TSTEP * (V + VK2 / 2.) VK3 = TSTEP * ACCEL (X + XK2 / 2., V + VK2 / 2., T + TSTEP / 2., W $ , G, Q) C THE FOURTH RUNGE KUTTA TERM, AT THE END POINT XK4 = TSTEP * (V * VK3) VK4 = TSTEP * ACCEL (X + XK3, V + VK3, T + TSTEP, W, G, Q) C COMBINE THE TERMS FOR THE OUTPUT XNEW = X + (XK1 + 2. * XK2 + 2. * XK3 + XK4) / 6. VNEW = V + (VK1 + 2. * VK2 + 2. * VK3 + VK4) / 6. RETURN END CCCCCCCCCCCCCCCC DOUBLE PRECISION FUNCTION ACCEL (X, V, T, W, G, Q) IMPLICIT DOUBLE PRECISION (A-H,O-Z) C SPECIAL CASE FOR A PHYSICAL PENDULUM ACCELERATION C C INPUT: C X - THE ANGLE OF THE PENDULUM, EVEN UP TO LARGE ANGLES C C V - THE VELOCITY OF THE PENDULUM C C T - THE TIME C C W - THE NATURAL FREQUENCY OF THE PENDULUM C C G - THE DRIVING FORCE C C Q - THE INVERSE-DAMPING (FOR NO DAMPING, USE 999999) * PRINT *, 'ACCEL' ACCEL = -SIN(X) - (V / Q) + (G * COS(W * T)) RETURN END