      PROGRAM ORBIT
C...
C...  GREENSPAN ORBIT PROBLEM
C...
C...  A SYSTEM OF FOUR NONLINEAR ALGEBRAIC EQUATIONS PROPOSED BY
C...  GREENSPAN, DONALD, A COUNTEREXAMPLE OF THE USE OF ENERGY AS
C...  A MEASURE OF COMPUTATIONAL EFFICIENCY, J. COMPUTATIONAL
C...  PHYSICS, V. 91, 490-494 (1990).  THE NUMERICAL SOLUTION IS
C...  BY NONLINEAR EQUATION SOLVER SNSQE TAKEN FROM KAHANER, DAVID,
C...  CLEVE MOLER AND STEPHEN NASH, NUMERICAL METHODS AND SOFTWARE,
C...  CHAPTER 7, PRENTICE HALL, ENGLEWOOD CLIFFS, NJ, 1989
C...
C...  DECLARATIVE STATEMENTS (NOTE - N SHOULD EQUAL NEQN READ FROM
C...  THE DATA FILE VIA FORMAT 902)
      PARAMETER (N=4,LWA=(3*N**2+13*N)/2)
C...
C...  VARIABLES AND PARAMETERS ARE IN COMMON FOR USE IN SUBORDINATE
C...  SUBROUTINES
      COMMON /T/      T,  NSTOP,  NORUN
     +       /Y/   Y(N)
     +       /F/  YT(N)
     +       /C/    Y10,    Y20,    Y30,    Y40,     TP
     +       /I/     IP,     NP
C...
C...  WORK ARRAY FOR SNSQE
      REAL WA(LWA)
C...
C...  FILE TO HOLD DOCUMENTATION TITLE
      CHARACTER TITLE(20)*4
C...
C...  EXTERNALS CALLED BY SNSQE
      EXTERNAL DERV, JAC
C...
C...  INIITIALIZE THE RUN COUNTER
      DATA NORUN/0/
C...
C...  INPUT/OUTPUT UNIT NUMBERS AND FILES
         NI=5
         NO=6
         OPEN(NI,FILE='DATA'  ,STATUS='OLD')
         OPEN(NO,FILE='OUTPUT',STATUS='NEW')
C...
C...  START THE NEXT RUN
4     NORUN=NORUN+1
C...
C...  READ DATA FILE, STARTING WITH DOCUMENTATION TITLE
      READ(NI,900)(TITLE(I),I=1,20)
C...
C...  CHECK FOR END OF RUNS (IN THE DATA FILE)
      IF((TITLE(1).EQ.'END ').AND.
     +   (TITLE(2).EQ.'OF R').AND.
     +   (TITLE(3).EQ.'UNS '))THEN
C...
C...     END OF RUNS WAS READ, SO TERMINATE EXECUTION
         STOP
C...
C...     CONTINUE RUN WITH CURRENT SET OF DATA
      END IF
C...
C...  READ INITIAL TIME, FINAL TIME AND PRINT INTERVAL OF TIME
      READ(NI,901)T0,TF,TP
C...
C...  READ NUMBER OF EQUATIONS, ERROR TOLERANCE
      READ(NI,902)NEQN,TOL
C...
C...  PRINT DATA SUMMARY
      WRITE(NO,903)NORUN,(TITLE(I),I=1,20),T0,TF,TP,NEQN,TOL
C...
C...  INITIAL CONDITIONS
         T=T0
         CALL INITAL
C...
C...  PRINT THE INITIAL CONDITIONS
      CALL PRINT(NI,NO)
C...
C...  PARAMETERS FOR SNSQE
         IOPT=2
         NPRINT=0
C...
C...  NUMBER OF STEPS AND POINTS ALONG THE SOLUTION
         NSTEPS=INT((TF-T0)/TP)
         NP=NSTEPS+1
C...
C...  STEP THROUGH THE NSTEPS STEPS
      DO 3 NSTEP=1,NSTEPS
C...
C...  TAKE A STEP ALONG THE SOLUTION
      CALL SNSQE(DERV,JAC,IOPT,N,Y,YT,TOL,NPRINT,INFO,WA,LWA)
C...
C...  UPDATE TIME
      T=T+TP
C...
C...  PRINT THE SOLUTION
      CALL PRINT(NI,NO)
C...
C...  UPDATE THE EQUATIONS
         Y10=Y(1)
         Y20=Y(2)
         Y30=Y(3)
         Y40=Y(4)
C...
C...  TAKE THE NEXT TIME STEP
3     CONTINUE
C...
C...  ALL TIME STEPS HAVE BEEN TAKEN (RUN IS COMPLETE), SO GO ON TO
C...  THE NEXT RUN BY READING THE NEXT SET OF DATA
      GO TO 4
900   FORMAT(20A4)
901   FORMAT(3E10.0)
902   FORMAT(I5,20X,E10.0)
903   FORMAT(1H1,/,' RUN NO. ',I3,' - ',20A4,//,
     + '  INITIAL VALUE OF TIME = ',E10.3,//,
     + '    FINAL VALUE OF TIME = ',E10.3,//,
     + ' PRINT INTERVAL OF TIME = ',E10.3,//,
     + ' SOLUTION BY SUBROUTINE SNSQE',//,
     + '    NUMBER OF EQUATIONS = ',I3,//,
     + '    ERROR TOLERANCE = ',E10.3,//)
      END
      SUBROUTINE INITAL
C...
C...  DECLARATIVE STATEMENTS
      PARAMETER (N=4)
C...
C...  COMMON AREAS DEFINED IN MAIN PROGRAM ORBIT
      COMMON /T/      T,  NSTOP,  NORUN
     +       /Y/   Y(N)
     +       /F/  YT(N)
     +       /C/    Y10,    Y20,    Y30,    Y40,     TP
     +       /I/     IP,     NP
C...
C...  INITIAL CONDITIONS
         Y10=0.5
         Y20=0.
         Y30=0.
         Y40=1.63
C...
C...  ESTIMATES OF THE SOLUTION FOR SUBROUTINE SNSQE AT T = 0
         Y(1)=Y10
         Y(2)=Y20
         Y(3)=Y30
         Y(4)=Y40
C...
C...  INITIALIZE A COUNTER FOR THE PLOTTED SOLUTION
      IP=0
      RETURN
      END
      SUBROUTINE DERV(N,Y,YT,IFLAG)
      REAL Y(N),  YT(N)
C...
C...  COMMON AREA FOR THE SOLUTION AT THE PRECEDING TIME STEP
      COMMON /C/    Y10,    Y20,    Y30,    Y40,     DT

C...  EQUATION (2.26)
         YT(1)=(Y(1)-Y10)/DT-(Y(2)+Y20)/2.0
C...
C...  EQUATION (2.27)
         YT(2)=(Y(3)-Y30)/DT-(Y(4)+Y40)/2.0
C...
C...  EQUATION (2.28)
         TERM0=SQRT( Y10**2+ Y30**2)
         TERM1=SQRT(Y(1)**2+Y(3)**2)
         TERM2=TERM0*TERM1
         TERM3=0.5*(TERM0+TERM1)
         YT(3)=(Y(2)-Y20)/DT+(1.0/TERM2)*(0.5*(Y(1)+Y10))/TERM3
C...
C...  EQUATION (2.29)
         YT(4)=(Y(4)-Y40)/DT+(1.0/TERM2)*(0.5*(Y(3)+Y30))/TERM3
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
C...
C...  DECLARATIVE STATEMENTS
      PARAMETER (N=4)
C...
C...  COMMON AREAS DEFINED IN MAIN PROGRAM ORBIT
      COMMON /T/      T,  NSTOP,  NORUN
     +       /Y/   Y(N)
     +       /F/  YT(N)
     +       /C/    Y10,    Y20,    Y30,    Y40,     TP
     +       /I/     IP,     NP
C...
C...  PRINT THE NUMERICAL SOLUTION
C...
C...     ENERGY, EQUATION (2.38)
         EN=ENERGY(N,Y)
C...
C...     L2 NORM USING SUBROUTINE DERV
         FNORM=0.
         CALL DERV(N,Y,YT,IFLAG)
         DO 2 I=1,N
            FNORM=FNORM+YT(I)**2
2        CONTINUE
         FNORM=SQRT(FNORM)
C...
C...     PRINT NUMERICAL SOLUTION, RESIDUAL, ENERGY AND NORM
         WRITE(NO,1)T,(I,Y(I),YT(I),I=1,N),EN,FNORM
1        FORMAT('     T = ',F7.3,//,
     +     4(2X,  '   I = ',I2,    '   Y = ',F9.6,2X,
     +            '  RE = ',E11.3,/),/,
     +       2X,  '  EN = ',F9.6,7X,'   FNORM = ',E11.3,///)
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 (2.38)
         R=SQRT(Y(1)**2+Y(3)**2)
         ENERGY=0.5*(Y(2)**2+Y(4)**2)-1.0/R
      RETURN
      END
      SUBROUTINE JAC
C...
C...  JAC IS A DUMMY ROUTINE TO SATISFY THE LOADER.  IT IS NOT CALLED
C...  BY SNSQE SINCE IOPT = 2
      RETURN
      END
      SUBROUTINE PLOT
C...
C...  DECLARATIVE STATEMENTS
      PARAMETER (N=4)
C...
C...  COMMON AREAS DEFINED IN MAIN PROGRAM ORBIT
      COMMON /T/      T,  NSTOP,  NORUN
     +       /Y/   Y(N)
     +       /F/  YT(N)
     +       /C/    Y10,    Y20,    Y30,    Y40,     TP
     +       /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 SNSQE - 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
