      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