      PROGRAM ORBIT2
C...
C...  FEYNMAN ORBIT PROBLEM - SOLUTION BY SUBROUTINE RKF45
C...
C...  SEE SUBROUTINE INITAL FOR A DISCUSSION OF THE PROBLEM.
C...
C...  THIS MAIN PROGRAM ESSENTIALLY CALLS SUBROUTINE RKF45 FOR THE
C...  SOLUTION OF A SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS (ODES).
C...
C...  THE MODEL INITIAL CONDITIONS ARE SET IN SUBROUTINE INITAL, AND
C...  THE MODEL DERIVATIVES ARE PROGRAMMED IN SUBROUTINE DERV.  THE
C...  NUMERICAL SOLUTION IS PRINTED AND PLOTTED IN SUBROUTINE PRINT.
C...
C...  SUBROUTINES INITAL, DERV AND PRINT.  THE FOLLOWING CODING IS FOR
C...  250 ORDINARY DIFFERENTIAL EQUATIONS (ODES).  IF MORE ODES ARE TO
C...  BE INTEGRATED, ALL OF THE 250*S SHOULD BE CHANGED TO THE REQUIRED
C...  NUMBER
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    Y(250)
     2      /F/    F(250)
C...
C...  THE NUMBER OF DIFFERENTIAL EQUATIONS IS IN COMMON/N/ FOR USE IN
C...  SUBROUTINE FCN
      COMMON/N/      NEQN
C...
C...  COMMON AREA TO PROVIDE THE INPUT/OUTPUT UNIT NUMBERS TO OTHER
C...  SUBROUTINES
      COMMON/IO/       NI,        NO
C...
C...  ABSOLUTE DIMENSIONING OF THE ARRAYS REQUIRED BY RKF45
      DIMENSION YV(250), WORK(1503), IWORK(5)
C...
C...  EXTERNAL THE DERIVATIVE ROUTINE CALLED BY RKF45
      EXTERNAL FCN
C...
C...  ARRAY FOR THE TITLE (FIRST LINE OF DATA), CHARACTERS  END OF RUNS
      CHARACTER TITLE(20)*4, ENDRUN(3)*4
C...
C...  DEFINE THE CHARACTERS  END OF RUNS
      DATA ENDRUN/'END ','OF R','UNS '/
C...
C...  DEFINE THE INPUT/OUTPUT UNIT NUMBERS
         NI=5
         NO=6
C...
C...  DEFINE INPUT/OUTPUT FILES
         OPEN(NI,FILE='DATA'  ,STATUS='OLD')
         OPEN(NO,FILE='OUTPUT',STATUS='NEW')
C...
C...  INITIALIZE THE RUN COUNTER
      NORUN=0
C...
C...  BEGIN A RUN
1     NORUN=NORUN+1
C...
C...  INITIALIZE THE RUN TERMINATION VARIABLE
      NSTOP=0
C...
C...  READ THE FIRST LINE OF DATA
      READ(NI,1000,END=999)(TITLE(I),I=1,20)
C...
C...  TEST FOR  END OF RUNS  IN THE DATA
      DO 2 I=1,3
         IF(TITLE(I).NE.ENDRUN(I))GO TO 3
2     CONTINUE
C...
C...  AN END OF RUNS HAS BEEN READ, SO TERMINATE EXECUTION
999   STOP
C...
C...  READ THE SECOND LINE OF DATA
3     READ(NI,1001,END=999)T0,TF,TP
C...
C...  READ THE THIRD LINE OF DATA
      READ(NI,1002,END=999)NEQN,TOL
C...
C...  PRINT A DATA SUMMARY
      WRITE(NO,1003)NORUN,(TITLE(I),I=1,20),T0,TF,TP,NEQN,TOL
C...
C...  INITIALIZE TIME
      T=T0
C...
C...  SET THE INITIAL CONDITIONS
      CALL INITAL
C...
C...  SET THE INITIAL CONDITIONS FOR SUBROUTINE RKF45
      DO 5 I=1,NEQN
         YV(I)=Y(I)
5     CONTINUE
C...
C...  SET THE PARAMETERS FOR SUBROUTINE RKF45
C...
C...     INITIAL TIME
         TV=T0
C...
C...     ERROR TOLERANCE
         ABSERR=0.
         RELERR=TOL
C...
C...     INITIALIZATION FOR MULTI-STEP MODE (RKF45 NOT CONSTRAINED
C...     BY OUTPUT INTERVAL TP)
         IFLAG=1
C...
C...     INITIALIZATION FOR ONE-STEP MODE (RKF45 CONSTRAINED BY
C...     OUTPUT INTERVAL TP)
C...     IFLAG=-1
C...
C...     FINAL TIME
         TEND=T0
C...
C...  CALL SUBROUTINE RKF45 TO START THE SOLUTION FROM THE INITIAL
C...  CONDITION (IFLAG = 1 OR -1) OR COMPUTE THE SOLUTION TO THE NEXT
C...  PRINT POINT (IFLAG = 2 OR -2)
4     CALL RKF45(FCN,NEQN,YV,TV,TEND,RELERR,ABSERR,IFLAG,WORK,IWORK)
C...
C...  PRINT THE INITIAL CONDITION (IFLAG = 1 OR -1) OR THE SOLUTION AT
C...  THE NEXT PRINT POINT (IFLAG = 2 OR -2)
      T=TV
      DO 6 I=1,NEQN
         Y(I)=YV(I)
6     CONTINUE
      CALL PRINT(NI,NO)
C...
C...  TEST FOR AN ERROR CONDITION
      IF((IFLAG.NE.2).AND.(IFLAG.NE.-2))THEN
C...
C...     PRINT A MESSAGE INDICATING AN ERROR CONDITION
         WRITE(NO,1004)IFLAG
C...
C...     GO ON TO THE NEXT RUN
         GO TO 1
      END IF
C...
C...  CHECK FOR A RUN TERMINATION
      IF(NSTOP.NE.0)GO TO 1
C...
C...  CHECK FOR THE END OF THE RUN
         TEND=TV+TP
         IF(TV.LT.(TF-0.5*TP))GO TO 4
C...
C...  THE CURRENT RUN IS COMPLETE, SO GO ON TO THE NEXT RUN
      GO TO 1
C...
C...  *****************************************************************
C...
C...  FORMATS
C...
1000  FORMAT(20A4)
1001  FORMAT(3E10.0)
1002  FORMAT(I5,20X,E10.0)
1003  FORMAT(1H1,
     1 ' RUN NO. - ',I3,2X,20A4,//,
     2 ' INITIAL T - ',E10.3,//,
     3 '   FINAL T - ',E10.3,//,
     4 '   PRINT T - ',E10.3,//,
     5 ' NUMBER OF DIFFERENTIAL EQUATIONS - ',I3,//,
     6 ' INTEGRATION ALGORITHM - RKF45 ',//,
     7 ' MAXIMUM INTEGRATION ERROR - ',E10.3,//,
     8 1H1)
1004  FORMAT(1H ,//,' IFLAG = ',I3,//,
     1 ' INDICATING AN INTEGRATION ERROR, SO THE CURRENT RUN'     ,/,
     2 ' IS TERMINATED.  PLEASE REFER TO THE DOCUMENTATION FOR'   ,/,
     3 ' SUBROUTINE',//,25X,'RKF45',//,
     4 ' FOR AN EXPLANATION OF THESE ERROR INDICATORS'            )
      END
      SUBROUTINE FCN(TV,YV,YDOT)
C...
C...  SUBROUTINE FCN IS AN INTERFACE ROUTINE BETWEEN SUBROUTINES RKF45
C...  AND DERV
C...
C...  ODE COMMON AREA
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(1)
     2      /F/      F(1)
C...
C...  THE NUMBER OF DIFFERENTIAL EQUATIONS IS AVAILABLE THROUGH COMMON
C...  /N/
      COMMON/N/      NEQN
C...
C...  ABSOLUTE DIMENSION THE DEPENDENT, DERIVATIVE VECTORS
      DIMENSION YV(250), YDOT(250)
C...
C...  TRANSFER THE INDEPENDENT VARIABLE, DEPENDENT VARIABLE VECTOR
C...  FOR USE IN SUBROUTINE DERV
      T=TV
      DO 1 I=1,NEQN
         Y(I)=YV(I)
1     CONTINUE
C...
C...  EVALUATE THE DERIVATIVE VECTOR
      CALL DERV
C...
C...  TRANSFER THE DERIVATIVE VECTOR FOR USE BY SUBROUTINE RKF45
      DO 2 I=1,NEQN
         YDOT(I)=F(I)
2     CONTINUE
      RETURN
      END
      SUBROUTINE INITAL
C...
C...  THE FOLLOWING SYSTEM OF ORDINARY DIFFERENTIAL EQUATIONS
C...  (ODES) FOR PLANETARY MOTION IS DISCUSSED BY FEYNMAN ET AL
C...  (1963)
C...
C...      2    2
C...     D X/DT = - X/R**3                                    (1)
C...
C...      2    2
C...     D Y/DT = - Y/R**3                                    (2)
C...
C...     R = (X**2 + Y**2)**0.5                               (3)
C...
C...     X(0) = 0.5, DX(0)/DT =  0                         (4)(5)
C...
C...     Y(0) = 0,   DY(0)/DT = 1.63                       (6)(7)
C...
C...  THE TOTAL ENERGY FOR THE SYSTEM, WHICH IS CONSTANT, IS GIVEN
C...  BY
C...
C...     E = (1/2)*(V **2 + V **2) - 1/R                      (8)
C...                 X       Y
C...  WHERE
C...
C...     V  = DX/DT, V  = DY/DT
C...      X           Y
C...
C...  E FROM EQUATION (8) IS COMPUTED AND PRINTED IN SUBROUTINE
C...  PRINT TO TEST THE OPERATION OF THE INTEGRATOR.
C...
C...  EQUATIONS (1) AND (2) CAN BE STATED IN TERMS OF NEW VARIABLES
C...  DEFINED AS
C...
C...     Y1 = X, Y2 = DX/DT
C...
C...     Y3 = Y, Y4 = DY/DT
C...
C...  THE FOLLOWING PROGRAMMING IS IN TERMS OF X1, X2, Y1 AND Y2.
C...
C...  EQUATIONS (1) TO (7) HAVE AN EXACT SOLUTION (GREENSPAN, 1990),
C...  WHICH, IF WRITTEN IN PARAMETRIC FORM (STENGLE, 1991), IS
C...
C...     X(T) = -C + A*COS(THETA(T))                          (9)
C...
C...     Y(T) = B*SIN(THETA(T))                              (10)
C...
C...  WHERE THETA(T) IS GIVEN BY THE ODE (2)
C...
C...     DTHETA/DT = SQRT((1/TERM1)*TERM2)                   (11)
C...
C...     TERM1 = (1/2)*((A*SIN(THETA))**2 + (B*COS(THETA))**2)   (12)
C...
C...     TERM2 = -1/(2*A) + 1/SQRT(TERM3)                    (13)
C...
C...     TERM3 = SQRT((A*COS(THETA) - C)**2 + (B*SIN(THETA))**2) (14)
C...
C...  THE INITIAL CONDITION FOR EQUATION (11) IS
C...
C...     THETA(0) = 0                                        (15)
C...
C...  THE SOLUTION TO EQUATIONS (1), (2) AND (11) IS BY ODE SOLVER
C...  RKF45 DISCUSSED IN FORSYTHE (1976).
C...
C...  REFERENCES:
C...
C...  FEYNMAN, RICHARD P., ROBERT B. LEIGHTON AND MATTHEW SANDS, LEC-
C...  TURES ON PHYSICS, ADDISON-WESLEY PUBLISHING COMPANY, READING,
C...  MA, PP 9-6 TO 9-9 (1963)
C...
C...  FORSYTHE, GEORGE E., MICHAEL A. MALCOLM, AND CLEVE B. MOLER,
C...  COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS, PRENTICE-HALL,
C...  TATIONS, PRENTICE-HALL, ENGLEWOOD CLIFF, NJ, PP 129-147 (1976)
C...
C...  GREENSPAN, DONALD, A COUNTEREXAMPLE OF THE USE OF ENERGY AS A
C...  MEASURE OF COMPUTATIONAL EFFICIENCY, J. COMPUTATIONAL PHYSICS,
C...  V. 91, 490-494 (1990)
C...
C...  STENGLE, GILBERT, PRIVATE COMMUNICATION, OCTOBER 8, 1991
C...
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)
     2      /F/   DYDT(N)
     3      /C/         R,         A,         B,         C
     4      /I/        IP,        NP
C...
C...  MODEL PARAMETERS
         A=10000.0/13431.0
         B=SQRT(664225.0/1343100.0)
         C=6569.0/26862.0
C...
C...  INITIAL CONDITIONS, EQUATIONS (4) TO (7), (11)
         Y(1)=0.5
         Y(2)=0.
         Y(3)=0.
         Y(4)=SQRT(2.0*(-1.0/(2.0*A)+1.0/(A-C)))
         Y(5)=0.
C...
C...  COMPUTE THE INITIAL DERIVATIVES
      CALL DERV
C...
C...  INITIALIZE INTEGER VARIABLES TO CONTROL THE PRINTING
         IP=0
         NP=201
      RETURN
      END
      SUBROUTINE DERV
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)
     2      /F/   DYDT(N)
     3      /C/         R,         A,         B,         C
     4      /I/        IP,        NP
C...
C...  RADIUS (WITH CUBING TO MINIMIZE CALCULATIONS)
         R=SQRT(Y(1)**2+Y(3)**2)**3
C...
C...  EQUATION (1)
         DYDT(1)=Y(2)
         DYDT(2)=-Y(1)/R
C...
C...  EQUATION (2)
         DYDT(3)=Y(4)
         DYDT(4)=-Y(3)/R
C...
C...  EQUATION (11)
         SINE=SIN(Y(5))
         COSE=COS(Y(5))
         TERM1=0.5*((A*SINE)**2+(B*COSE)**2)
         TERM2=-1.0/(2.0*A)+1.0/SQRT((A*COSE-C)**2+(B*SINE)**2)
         DYDT(5)=SQRT((1.0/TERM1)*TERM2)
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)
     2      /F/   DYDT(N)
     3      /C/         R,         A,         B,         C
     4      /I/        IP,        NP
      REAL YE(N),   ER(N)
C...
C...  PRINT THE NUMERICAL SOLUTION
C...
C...     ENERGY, EQUATION (8)
         EN=ENERGY(N,Y)
C...
C...     EXACT SOLUTION AND ERROR IN NUMERICAL SOLUTION
C...     (THETA = Y(5), DTHETA/DT = DYDT(5) IN EQUATION (11))
         THETA=   Y(5)
         DTHDT=DYDT(5)
         YE(1)=-C+A*COS(THETA)
         YE(2)=-A*SIN(THETA)*DTHDT
         YE(3)=   B*SIN(THETA)
         YE(4)= B*COS(THETA)*DTHDT
         DO 2 I=1,N-1
            ER(I)=Y(I)-YE(I)
2        CONTINUE
C...
C...     PRINT NUMERICAL AND EXACT SOLUTIONS, NUMERICAL ERROR AND
C...     ENERGY
         WRITE(NO,1)T,THETA,(I,Y(I),YE(I),ER(I),I=1,N-1),EN
1        FORMAT('     T = ',F6.3,'   THETA = ',F6.3,//,
     +     4(2X,  '   I = ',I2,      '   Y = ',F9.6,
     +            '  YE = ',F9.6,    '  ER = ',E9.3,/),/,2X,
     +            '  EN = ',F9.6,///)
C...
C...  STORE THE NUMERICAL SOLUTION FOR PLOTTING
         IP=IP+1
         CALL PLOT
      RETURN
      END
      FUNCTION ENERGY(N,Y)
C...
C...  FUNCTION ENERGY COMPUTES THE TOTAL ENERGY OF THE SYSTEM
C...
      REAL Y(N)
C...
C...  ENERGY, EQUATION (8)
         R=SQRT(Y(1)**2+Y(3)**2)
         ENERGY=0.5*(Y(2)**2+Y(4)**2)-1.0/R
      RETURN
      END
      SUBROUTINE PLOT
      PARAMETER (N=5)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/      Y(N)
     2      /F/   DYDT(N)
     3      /C/         R,         A,         B,         C
     4      /I/        IP,        NP
C...
C...  INITIATE PLOTTING
      IF(IP.EQ.1)THEN
C...
C...     OPEN FILE FOR TOP DRAWER PLOTTING      
         IT=4
         OPEN(IT,FILE= 'ORBIT.TOP',STATUS='UNKNOWN')
C...
C...     SCALE THE AXES
         WRITE(IT,7)
7        FORMAT(' SET LIMITS X FROM -1 TO 1 Y FROM -1 TO 1',/,
     1          ' SET FONT DUPLEX')
         WRITE(IT,4)
4        FORMAT(' SET WINDOW X 2 TO 8 Y 2 TO 8')
      END IF
C...
C...  WRITE Y(3) VS Y(1) FOR PLOTTING
         WRITE(IT,5)Y(1),Y(3)
5        FORMAT(2F10.4)
C...
C...  COMPLETE PLOTTING
      IF(IP.EQ.NP)THEN
C...
C...     CONNECT POINTS
         WRITE(IT,6)
6        FORMAT(' JOIN 1')
C...
C...     LABEL THE AXES
         WRITE(IT,8)NP
8        FORMAT(
     1   ' TITLE 2.75 8.5 "',I3,'-point orbit from RKF45 - SP"'
     2   ,/,' TITLE LEFT "       y3(t)"'
     3   ,/,' TITLE BOTTOM "y1(t)"'
     4   ,/,' TITLE 2.5 0.7 "Greenspan/Feynman Orbit Problem"')
C...
C...     END PLOT
         WRITE(IT,10)
10       FORMAT(' NEW FRAME')
      END IF
      RETURN
      END
