;+
; Project     :	Computational Physics Final Exam
;
; Name        :	ORBITSFINAL
;
; Purpose     :	Calculate the orbit of a planet/comet of the Solar System.
;
; Explanation :	This procedure is designed to solve the differential equation
;               that describes the orbit of a planet, comet, asteroid about the
;               Sun.  It solves the equation of motion for Newton's law of
;               gravity, using the Sun as the center mass.  In order to avoid
;               numeric overflows and underflows, calculations are made
;               with the radius vector in units of AU and the orbital velocity
;               in units of AU/yr.  Then, the constant G*M-sun = v^2 r = 4\pi^2
;               AU^3/yr^2, where we have used the fact that the Earth's
;               average velocity is C/T = (2\pi a)/(1 yr) = 2\pi AU/yr, where
;               C is the circumference of the Earth's orbit (assuming that it
;               is circular), T is the orbital period, and a = 1 AU is the
;               Earth's semimajor axis (or radius if a circular orbit is
;               assumed).  Note that even though we use these units here, we
;               are not assuming that the Earth follows an exactly circular
;               orbit.
;
; Use         :	ORBITS [, M, A, ECC, T]
;
; Inputs      :	None.
;
; Opt. Inputs :	M:    Mass of planet in Earth masses.  (Default = 1).
;
;               A:    Semimajor axis of planet in AU.  (Default = 1).
;
;               ECC:  Eccentricty of orbit.  (Default = 0.0167 [Earth]).
;
; Outputs     :	None.
;
; Opt. Outputs:	T:    Period of the orbit (in days).
;
; Keywords    :	NSTEPS:	The number of steps per orbit in the integration.
;                       (Default = 5000).
;
;               NORBIT: Number of orbits the calculation is to make.  (Default
;                       = 5).  Note: NSTEPS*NORBIT cannot exceed 32,000!
;
;               PLANET: Possible values:
;
;                       'MERCURY': Compute for planet Mercury.
;
;                       'VENUS':   Compute for planet Venus.
;
;                       'EARTH':   Compute for planet Earth.  (Default)
;
;                       'MARS':	   Compute for planet Mars.
;
;                       'JUPITER': Compute for planet Jupiter.
;
;                       'SATURN':  Compute for planet Saturn.
;
;                       'URANUS':  Compute for planet Uranus.
;
;                       'NEPTUNE': Compute for planet Neptune.
;
;                       'PLUTO':   Compute for planet Pluto.
;
;                       'HALLEY':  Compute for comet Halley.
;
;                   INTEG:  Integration scheme:
;
;                       'RK':   Runge-Kutta 4th Order Integration (note that
;                               you must supply additional function ORBEQN
;                               that calculates the derivatives of your
;                               equations).  [DEFAULT]
;
;                       'EULER': Euler method (round-off errors cause nonclosed
;                               orbits).
;
;                       'EC':   Euler-Cromer method (still causes open orbits
;                               unless NSTEPS is large (which means TAU is
;                               small).
;
;                   ERRMSG:	If defined and passed, then any error messages will
;                           be returned to the user in this parameter rather
;                           than being printed to the screen.  If no errors are
;                           encountered, then a null string is returned.  In
;                           order to use this feature, the string ERRMSG must
;                           be defined first, e.g.,
;
;                               ERRMSG = ''
;                               ORBITSFINAL, M, A, ECC, ERRMSG=ERRMSG
;                               IF ERRMSG(0) NE '' THEN ...
;
; Calls       :	ORBEQNFINAL (if Runge-Kutta method selected), ORBPLOT.  Note
;               ORBPLOT also calls procedure ORBTICKS.
;
; Common      :	None.
;
; Restrictions:	None.
;
; Side effects:	None.
;
; Category    :	Computation physics.
;
; Prev. Hist. :	None.
;
; Written     :	Donald G. Luttermoser, ETSU/Physics, 9 December 2002
;
; Modified    :	Version 1, Donald G. Luttermoser, ETSU/Physics, 9 Dec 2002
;                   Initial program.
;               Version 2, Donald G. Luttermoser, ETSU/Physics, 3 May 2007
;                   Replaced the IDL MESSAGE utility with a PRINT statement.
;                   Added a print-out of the total energy of the orbit.
;               Version 3, Donald G. Luttermoser, ETSU/Physics, 9 Dec 2007
;                   Use A.U. for distances and years for time instead of
;                   meters and seconds to improve round-off and overflow/
;                   underflow errors.
;               Version 4, Donald G. Luttermoser, ETSU/Physics, 23 Nov 2009
;                   Changed procedure name from ORBITS to ORBITSFINAL for
;                   distribution to the Computational Physics students.
;                   Changed the default integration scheme from Euler-Cromer
;                   to 4th Order Runge-Kutta.  Changed variable array indicies
;                   symbols from "()" to "[]".
;
; Version     :	Version 4,  23 November 2008.
;
;-
;
PRO ORBITSFINAL, M, A, ECC, T, NSTEPS=NSTEPS, NORBIT=NORBIT, PLANET=PLANET, $
    INTEG=INTEG, ERRMSG=ERRMSG
;
EMESS = ''      ; Set to non-null string if error is encountered.
TERM = !D.NAME  ; Either 'X' for X-Windows or 'WIN' for Microsoft Windows.
;
; Internal variables.
;
MASS = 1.0D0  ; Default planet mass in Earth masses.
ASEMI = 1.0D0 ; Default semimajor axis in A.U. (Astronomical Units).
ORBECC = 0.0167D0   ; Default orbit eccentricity (= to Earth's)
MEARTH = 5.9763D24  ; Mass of the Earth in kilograms.
REARTH = 6.367E6    ; Radius of the Earth in meters.
MSUN = 1.9891D30    ; Mass of the Sun in kilograms.
G = 6.6720D-11      ; Newton's Universal Gravity Constant in SI units.
;
TYPE = 3   ; Assumes that the Earth is the default planet.
QINT = 0   ; Uses 4th-order Runge-Kutta for the default integration scheme.
;
; Check input.
;
CASE N_PARAMS() OF
    4:  BEGIN
        MASS = DOUBLE(M)
        ASEMI = DOUBLE(A)
        ORBECC = DOUBLE(ECC)
        TYPE = -3
        END
    3:  BEGIN
        MASS = DOUBLE(M)
        ASEMI = DOUBLE(A)
        ORBECC = DOUBLE(ECC)
        TYPE = -3
        END
;
    2:  BEGIN
        MASS = DOUBLE(M)
        ASEMI = DOUBLE(A)
        TYPE = -2
        END
;
    1:  BEGIN
        MASS = DOUBLE(M)
        TYPE = -1
        END
;
    0:  TYPE = 3
;
    ELSE:  BEGIN
        EMESS = 'Syntax:  ORBITS [, M, A, ECC, T]'
        GOTO, HANDLE_ERROR
        END
ENDCASE
;
; Planet information.
;
THEPLANETS = ['Hypothetical', 'Mercury', 'Venus', 'Earth', 'Mars', $
    'Jupiter', 'Saturn', 'Uranus', 'Neptune', 'Pluto', 'Halley']
PLECC = [0.0167, 0.2056, 0.068, 0.0167, 0.0934, 0.0485, 0.0556, 0.0472, $
    0.0086, 0.250, 0.967]
PLMASS = [1., 0.0558, 0.815, 1.000, 0.107, 317.89, 95.184, 14.536, 17.148, $
    0.0022, 1.70E-11]   ; In units of Earth masses.
PLSEMIMAJOR = [1., 0.387, 0.723, 1.000, 1.524, 5.203, 9.555, 19.22, 30.11, $
    39.44, 17.95]       ; In A.U.
;
; Check keywords.
;
IF KEYWORD_SET(NSTEPS) THEN NPTS = NSTEPS ELSE NPTS = 5000
IF KEYWORD_SET(NORBIT) THEN NORB = NORBIT ELSE NORB = 5
IF KEYWORD_SET(ECC) THEN ORBECC = DOUBLE(ECC) ELSE ORBECC = 0.0167D0
IF KEYWORD_SET(PLANET) THEN NPLANET = STRUPCASE(PLANET) ELSE NPLANET = 'EARTH'
TYPE = WHERE(STRUPCASE(THEPLANETS) EQ NPLANET, NTYPE)
IF NTYPE EQ 0 THEN BEGIN
    EMESS = 'Planet ' + NPLANET + ' does not exist in the data base!'
    GOTO, HANDLE_ERROR
ENDIF ELSE TYPE = TYPE[0]
IF KEYWORD_SET(INTEG) THEN BEGIN
    CASE STRUPCASE(INTEG) OF
        'EULER':  QINT = 1
        'EC':     QINT = 2
        ELSE:     QINT = 0  ; 'RK' Runge-Kutta integration.
    ENDCASE
ENDIF
;
; Internal parameters.
;
ASK = 'Y'
PASS = 0   ; PASS = 0:  graphics on terminal,
;            PASS = 1:  graphics to postscript file.
;
INTYPE = ['4th-Order Runge-Kutta', 'Euler', 'Euler-Cromer']
FNAME = ['orbrk', 'orbeuler', 'orbec']
FSUFFIX = ['.ps', '.eps']
INTSZ = [0.8, 1.1, 1.1]
INTUSE = INTYPE[QINT]
NTOT = NPTS * NORB                ; Total number of orbit calculations.
GMCONST = 1.327393D20             ; Grav. constant * mass of Sun in m^3/s^2.
GM = 4.0D0 * !DPI^2               ; Grav. constant * mass of Sun in AU^3/yr^2.
AU = 1.4959787061D11              ; Astronomical Unit in meters.
SECPERYR = 3.15569259747D7        ; Number of seconds per year.
SECPERDAY = 8.6400D4              ; Number of seconds per day.
VCONVER = AU * 1.0D-3 / SECPERYR  ; Conversion from AU/yr to km/s.
;
MLAB = 'Solar Orbit'              ; Make sure to set font 17 in main title.
TLAB = '!7h !6(degrees)'          ; Use font 6 in x-axis label.
RLAB = '!6r (AU)'                 ; Use font 6 in y-axis label.
VLAB = '!6Velocity (km/s)'        ; Use font 6 in y-axis label.
;
; Set orbital parameters.  For negative types, assume unpassed parameters are
; the same as Earth's parameters.
;
CASE TYPE OF
    -3: TYPE = 0	; Do nothing.
;
    -2: BEGIN
        ORBECC = DOUBLE(PLECC[3])
        TYPE = 0
        END
;
    -1: BEGIN
        ASEMI = DOUBLE(PLSEMIMAJOR[3])
        ORBECC = DOUBLE(PLECC[3])
        TYPE = 0
        END
;
    ELSE: BEGIN
          MASS = DOUBLE(PLMASS[TYPE])
          ASEMI = DOUBLE(PLSEMIMAJOR[TYPE])
          ORBECC = DOUBLE(PLECC[TYPE])
          END
ENDCASE
;
IF TYPE EQ 0 THEN PLNAME = 'Hypothetical' ELSE $
    PLNAME = THEPLANETS[TYPE] + "'s"
MLAB = '!17' + PLNAME + ' ' + MLAB
;
; Perihelion and aphelion distances.  Make a temporary circular orbit.
;
RPERI = ASEMI * (1.0D0 - ORBECC)
RAP = ASEMI * (1.0D0 + ORBECC)
BSEMI = ASEMI * SQRT(1.0D0 - ORBECC^2)
IF ORBECC GE 0.7 THEN SCPERI = 10. ELSE SCPERI = 1.5
;
; Total energy of the orbit in Joules.
;
TOTENERGY = GMCONST * MASS * MEARTH / (2 * ASEMI * AU)
PRINT, 'Total orbital energy of ' + PLNAME + ' is ' + $
    STRING(TOTENERGY, FORMAT='(E10.3)') + ' Joules.'
;
; Initial conditions (at the aphelion point), all calcs done in AU and yrs.
;
PERCALC = SQRT(ASEMI^3)	; Period (in years) based upon Kepler's 3rd Law.
TAU = (PERCALC / NPTS)	; Time step of integration in years.
TAUPR = TAU * 365.24	; Time step to be printed in days.
;
; Students:  To address this maximum value problem, think about why I have the
; second condition of this product being less than zero.  This should give you
; a hint on how to overcome this maximum number of time steps problem.
;
IF (NPTS*NORB GT 32000) OR (NPTS*NORB LT 0) THEN BEGIN
    EMESS = 'NSTEPS*NORBIT must be less than 32000!'
    GOTO, HANDLE_ERROR
ENDIF
;
; Variable arrays.  Note that here, THETA = 0 occurs when r = R-aphelion.
;
RCART = DBLARR(2)  	; r stored in Cartesian x & y coordinates.
VCART = RCART	    ; v stored in Cartesian x & y coordinates.
;
; Define internal arrays.
;
VEL = DBLARR(NTOT)
DIST = VEL
TIME = VEL
THETA = VEL    ; Angle in radians wrt aphelion distance.
DIST[0] = RAP  ; Initial distance from Sun in A.U.s.
VEL[0] = SQRT(GM * (2.D0/DIST[0] - 1.D0/ASEMI))  ; Initial velocity in AU/yr.
VMAX = SQRT(GM * (2.D0/RPERI - 1.D0/ASEMI)) * VCONVER	; Store as km/s.
VMIN = SQRT(GM * (2.D0/RAP - 1.D0/ASEMI)) * VCONVER	; Store as km/s.
;
; RCART[0] contains the x distance from the Sun in AUs.
; RCART[1] contains the y distance from the Sun in AUs.
;
RCART(0) = DIST[0]
RCART[1] = 0.D0
;
; VCART[0] contains the v(x) orbital velocity in AU/yr.
; VCART[1] contains the v(y) orbital velocity in AU/yr.
;
VCART[0] = 0.D0
VCART[1] = VEL[0]
;
; Set up plot to display orbit in real time.
;
IF !D.NAME EQ 'WIN' THEN BEGIN
   WINDOW, 0, XSIZE=400, YSIZE=320
ENDIF ELSE BEGIN
   WINDOW, 0, XSIZE=500, YSIZE=400
ENDELSE
;
ORBPLOT, ASEMI, ORBECC, MLAB, 'AU', HARDCOPY=PASS, ERRMSG=EMESS
XS0 = !X.S
YS0 = !Y.S
IF EMESS NE '' THEN GOTO, HANDLE_ERROR
;
IF !D.NAME EQ 'WIN' THEN BEGIN
   WINDOW, 1, XSIZE=400, YSIZE=320
ENDIF ELSE BEGIN
   WINDOW, 1, XSIZE=500, YSIZE=400
ENDELSE
;
PLOT, [0.95*RPERI, 1.05*RAP], [0.95*VMIN, 1.05*VMAX], XSTYLE=1, YSTYLE=1+16, $
    TITLE='!17Orbital Velocity', XTITLE=RLAB, YTITLE=VLAB, /NODATA
XS1 = !X.S
YS1 = !Y.S
;
IF !D.NAME EQ 'WIN' THEN BEGIN
   WINDOW, 2, XSIZE=400, YSIZE=320
ENDIF ELSE BEGIN
   WINDOW, 2, XSIZE=500, YSIZE=400
ENDELSE
;
PLOT, [-10., 370.], [0.95*VMIN, 1.05*VMAX], TITLE='!17Orbital Velocity', $
    XSTYLE=1, YSTYLE=1+16, XTITLE=TLAB, YTITLE=VLAB, /NODATA
XS2 = !X.S
YS2 = !Y.S
XYOUTS, 0.2, 0.85, '!6N(steps/orbit) = '+STRTRIM(STRING(NPTS),2), /NORM, $
    SIZE=1.6
XYOUTS, 0.2, 0.80, '!6N(orbit) = '+STRTRIM(STRING(NORB),2), /NORM, SIZE=1.6
;
IF !D.NAME EQ 'WIN' THEN BEGIN
   WINDOW, 3, XSIZE=400, YSIZE=320
ENDIF ELSE BEGIN
   WINDOW, 3, XSIZE=500, YSIZE=400
ENDELSE
;
PLOT, [0., NORB*PERCALC], [0.95*RPERI, 1.05*RAP], XSTYLE=1, YSTYLE=1+16, $
    TITLE='!17Planetary Distance', XTITLE='!6Time (yr)', YTITLE=RLAB, $
    /NODATA
XS3 = !X.S
YS3 = !Y.S
;
; Calculate angles, velocities and distances of NORB orbits (NORB*NTOT)
; using any of the three listed differential equation solvers.
;
TSTART = SYSTIME(1)
IQPRINT = 0
LASTRAT = 10L
FOR ISTEP = 1, NTOT-1 DO BEGIN
    ROLD = RCART
    VOLD = VCART
    TIME[ISTEP] = TIME[ISTEP-1] + TAU
    RNORM = SQRT(RCART[0]^2 + RCART[1]^2)
    ACCEL = -GM * (RCART / RNORM^3)
;
    CASE QINT OF
;
; Euler's method.
;
        1:  BEGIN
            RCART = RCART + TAU * VCART
            VCART = VCART + (TAU * ACCEL)
            END
;
; Euler-Cromer's method.
;
        2:  BEGIN
            VCART = VCART + (TAU * ACCEL)
            RCART = RCART + TAU * (VCART + VOLD) / 2.
            END
;
; 4-th Order Runge-Kutta (must supply separate function 'ORBEQN').
;
        ELSE: BEGIN
              STATE = [RCART, VCART]
              DYDT = ORBEQNFINAL(TIME[ISTEP], STATE, ERRMSG=EMESS)
              IF EMESS NE '' THEN GOTO, HANDLE_ERROR
              RESULT = RK4(STATE, DYDT, TIME[ISTEP], TAU, 'orbeqnfinal', $
                  /DOUBLE)
              RCART[0] = RESULT[0]
              RCART[1] = RESULT[1]
              VCART[0] = RESULT[2]
              VCART[1] = RESULT[3]
              END
    ENDCASE
;
; Store the current angle (in radians), distance, and orbital velocity.
;
    THETA[ISTEP] = ATAN(RCART[1], RCART[0])
    IF THETA[ISTEP] LT 0. THEN THETA[ISTEP] = 2.0D0*!DPI + THETA[ISTEP]
    DIST[ISTEP] = SQRT(RCART[0]^2 + RCART[1]^2)
    VEL[ISTEP] = SQRT(VCART[0]^2 + VCART[1]^2)
;
    WSET, 0
    !X.S = XS0
    !Y.S = YS0
    OPLOT, [ROLD[0], RCART[0]], [ROLD[1], RCART[1]], PSYM=2, SYMSIZE=0.2
    WSET, 1
    !X.S = XS1
    !Y.S = YS1
    OPLOT, [DIST[ISTEP-1], DIST[ISTEP]], $
        [VEL[ISTEP], VEL[ISTEP-1]]*VCONVER, PSYM=2, SYMSIZE=0.2
    WSET, 2
    !X.S = XS2
    !Y.S = YS2
    OPLOT, [THETA[ISTEP-1]*!RADEG, THETA[ISTEP]*!RADEG], $
        [VEL[ISTEP], VEL[ISTEP-1]]*VCONVER, PSYM=2, SYMSIZE=0.2
    WSET, 3
    !X.S = XS3
    !Y.S = YS3
    OPLOT, [TIME[ISTEP-1], TIME[ISTEP]], [DIST[ISTEP-1], $
        DIST[ISTEP]], PSYM=2, SYMSIZE=0.2
;
; Let the user know where the integration is current at.
;
    RAT = 100L * ISTEP / NTOT
    IF RAT GE LASTRAT THEN IQPRINT = 1
    IF IQPRINT EQ 1 THEN BEGIN
        PRINT, 'The integration of the solution is ' + $
            STRING(RAT, FORMAT='(I3)') + '% complete.'
        IQPRINT = 0
        LASTRAT = LASTRAT + 10L
    ENDIF
ENDFOR
TEND = SYSTIME(1)
CPUTIME = TEND - TSTART
TCPU = STRTRIM(STRING(CPUTIME, FORMAT='(F10.2)'), 2) + ' seconds.'
PRINT, INTUSE + ' integration scheme took ' + TCPU
;
IF STRUPCASE(TERM) EQ 'TEK' THEN BEGIN
    PRINT, '*************************************************************'
    PRINT, '               Hit <return> to continue!'
    PRINT, '*************************************************************'
    READ, ASK
ENDIF ELSE BEGIN
    QRES = DIALOG_MESSAGE('Click OK to continue!', /INFORMATION, $
        TITLE='Calculation Done')
ENDELSE
;
; Calculate orbital period, perihelion position, aphelion position, v (1AU)
; from the calculated values.
;
RPERIOUT = MIN(DIST, IMIN)
RAPOUT = MAX(DIST, IMAX)
IBIG = WHERE(DIST GE 1.D0, NBIG)
IF NBIG EQ 0 THEN BEGIN
    VELPOS = VEL[IMAX] * VCONVER
    PLPOS = RAPOUT
ENDIF ELSE BEGIN
    VELTMP = VEL[IBIG]
    DISTTMP = DIST[IBIG]
    TMP = MIN(DISTTMP, ITMP)
    VELPOS = VELTMP[ITMP] * VCONVER
    PLPOS = 1.D0
ENDELSE
NP = N_ELEMENTS(TIME)
PERNUM = -1.D0
FOR I = 2, NP-3 DO BEGIN
    IF (DIST[I]-DIST[I-1] GE 0.D0) AND (DIST[I+1]-DIST[I] LT 0.D0) THEN $
        PERNUM = [PERNUM, TIME[I]]
ENDFOR
CALCPER = 0.D0
NCP = N_ELEMENTS(PERNUM)
IF NCP GT 1 THEN BEGIN
    FOR I = 1, NCP-2 DO CALCPER = CALCPER + (PERNUM[I+1] - PERNUM[I])
    CALCPER = CALCPER / (NCP-2)
ENDIF
;
WDELETE, 0, 1, 2, 3
WAIT, 1
;
; Set plotting characteristics.
;
PLOTIT:
!P.CHARSIZE = 0.8
!X.CHARSIZE = 1.0
!Y.CHARSIZE = 1.0
!P.MULTI = [0, 2, 2]
;
; Plot 1 (upper left):
; Plot the integration.  First the orbit in Cartesian coordinates.
;
ORBPLOT, ASEMI, ORBECC, MLAB, 'AU', HARDCOPY=PASS, ERRMSG=EMESS
IF PASS EQ 0 THEN BEGIN
    XPL = DIST * COS(THETA)
    YPL = DIST * SIN(THETA)
ENDIF
OPLOT, XPL, YPL, THICK=2*PASS+2
;
; Plot 2 (upper right):
; Now v versus r to make the shape is linear (i.e., a converged solution).
;
VPL = VEL * VCONVER	; Convert to km/s from AU/yr.
PLOT, DIST, VPL, TITLE='!17'+PLNAME+' Orbital Velocity', $
    XTITLE=RLAB, YTITLE=VLAB, XSTYLE=2, YSTYLE=2+16
X2L = !X.CRANGE[0]
X2R = !X.CRANGE[1]
X2DIFF = X2R - X2L
Y2B = !Y.CRANGE[0]
Y2T = !Y.CRANGE[1]
Y2DIFF = Y2T - Y2B
XYOUTS, X2L+0.5*X2DIFF, Y2T-0.10*Y2DIFF, '!6N(steps/orbit) = ' + $
    STRTRIM(STRING(NPTS), 2), SIZE=1.0
XYOUTS, X2L+0.5*X2DIFF, Y2T-0.17*Y2DIFF, '!6N(orbit) = ' + $
    STRTRIM(STRING(NORB), 2), SIZE=1.0
XYOUTS, X2L+0.5*X2DIFF, Y2T-0.24*Y2DIFF, '!7s !6= ' + STRTRIM(STRING(TAUPR, $
    FORMAT='(F10.4)'), 2) + ' days', SIZE=1.0
XYOUTS, X2L+0.5*X2DIFF, Y2T-0.31*Y2DIFF, '!6E = ' + STRING(TOTENERGY, $
    FORMAT='(E10.3)') + ' Joules', SIZE=1.0
XYOUTS, X2L+0.05*X2DIFF, Y2B+0.41*Y2DIFF, '!6Output from ODE', SIZE=1.0
XYOUTS, X2L+0.05*X2DIFF, Y2B+0.34*Y2DIFF, '!6Numerical Solution', SIZE=1.0
XYOUTS, X2L+0.05*X2DIFF, Y2B+0.25*Y2DIFF, '!6r!Dp!N = ' + $
    STRTRIM(STRING(RPERIOUT,FORMAT='(F9.2)'), 2) + ' AU', SIZE=1.0
XYOUTS, X2L+0.05*X2DIFF, Y2B+0.18*Y2DIFF, '!6r!Da!N = '+STRTRIM(STRING(RAPOUT, $
    FORMAT='(F9.2)'), 2) + ' AU', SIZE=1.0
XYOUTS, X2L+0.05*X2DIFF, Y2B+0.11*Y2DIFF, '!6P = '+STRTRIM(STRING(CALCPER, $
    FORMAT='(F9.2)'),2) + ' yr', SIZE=1.0
XYOUTS, X2L+0.05*X2DIFF, Y2B+0.04*Y2DIFF, '!6v(at '+STRTRIM(STRING(PLPOS, $
    FORMAT='(F6.1)'), 2) + ' AU) = ' + STRTRIM(STRING(VELPOS, $
    FORMAT='(F9.2)'), 2) + ' km/s', SIZE=1.0
;
; Plot 3 (lower left):
; Now the v(theta).
;
PLOT, THETA*!RADEG, VPL, TITLE='!17Orbital Velocity', $
    XSTYLE=2, YSTYLE=2+16, XTITLE=TLAB, YTITLE=VLAB, PSYM=2, SYMSIZE=0.5
X3L = !X.CRANGE[0]
X3R = !X.CRANGE[1]
X3DIFF = X3R - X3L
Y3B = !Y.CRANGE[0]
Y3T = !Y.CRANGE[1]
Y3DIFF = Y3T - Y3B
XYOUTS, X3L+0.05*X3DIFF, Y3T-0.10*Y3DIFF, '!6Integration Method: '+INTUSE, $
    SIZE=INTSZ(QINT)
XYOUTS, X3L+0.05*X3DIFF, Y3B+0.06*Y3DIFF, '!6Integration Time = '+TCPU, $
    SIZE=INTSZ[QINT]
;
; Plot 4 (lower right):
; Finally, plot r(t) to see the accuracy of this technique.
;
PLOT, TIME, DIST, XSTYLE=2, YSTYLE=2+16, TITLE='!17'+PLNAME+$
    ' Distance from Sun', XTITLE='!6Time (yr)', YTITLE=RLAB
X4L = !X.CRANGE[0]
X4R = !X.CRANGE[1]
X4DIFF = X4R - X4L
Y4B = !Y.CRANGE[0]
Y4T = !Y.CRANGE[1]
Y4DIFF = Y4T - Y4B
XYOUTS, X4L+0.05*X4DIFF, Y4T-0.10*Y4DIFF, '!6Output from', SIZE=1.0
XYOUTS, X4L+0.05*X4DIFF, Y4T-0.17*Y4DIFF, '!6Analytic Solution', SIZE=1.0
XYOUTS, X4L+0.05*X4DIFF, Y4T-0.26*Y4DIFF, '!6r!Dp!N = ' + $
    STRTRIM(STRING(RPERI,FORMAT='(F9.2)'), 2) + ' AU', SIZE=1.0
XYOUTS, X4L+0.05*X4DIFF, Y4T-0.33*Y4DIFF, '!6r!Da!N = '+STRTRIM(STRING(RAP, $
    FORMAT='(F9.2)'), 2) + ' AU', SIZE=1.0
XYOUTS, X4L+0.05*X4DIFF, Y4T-0.41*Y4DIFF, '!6P = '+STRTRIM(STRING(PERCALC, $
    FORMAT='(F9.2)'),2) + ' yr', SIZE=1.0
;
; Save as hardcopy?
;
IF PASS EQ 0 THEN BEGIN
    IF STRUPCASE(TERM) EQ 'TEK' THEN BEGIN
        READ, 'Do you wish to save this plot in a postscript file ' + $
            '(Y/N)? [Y] ', ASK
        IF STRUPCASE(ASK) EQ 'Y' THEN READ, 'Do you wish this postscript '+$
            'plot file to be encapsulated (Y/N)? [N] ', APS
        IF STRUPCASE(APS) EQ 'Y' THEN QPS = 1 ELSE QPS = 0
    ENDIF ELSE BEGIN
        ASK = DIALOG_MESSAGE(['Do you wish to save this', $
            'plot in a postscript file?'], /QUESTION, $
            TITLE='Save Resulting Plots?')
        IF STRUPCASE(ASK) EQ 'YES' THEN BEGIN
            APS = DIALOG_MESSAGE(['Do you wish to encapsulate your', $
                'postscript plot file?'], /QUESTION, $
                TITLE='Encapsulate plot?')
            IF STRUPCASE(APS) EQ 'YES' THEN QPS = 1 ELSE QPS = 0
        ENDIF ELSE QPS = 0
    ENDELSE
    IF STRUPCASE(STRMID(STRTRIM(ASK, 2), 0, 1)) NE 'N' THEN BEGIN
        SET_PLOT, 'PS'
        FNOUT = FNAME(QINT)
        DEVICE, /LANDSCAPE, ENCAPSULATE=QPS, FILENAME=FNOUT+FSUFFIX[QPS]
        PRINT, 'Plots are being saved to:  ', FNOUT+FSUFFIX[QPS]
        PASS = 1
        GOTO, PLOTIT
    ENDIF
ENDIF ELSE BEGIN
    DEVICE, /CLOSE
    PASS = 0
    SET_PLOT, TERM    ; Defined above.
ENDELSE
!P.MULTI = 0
PRINT, 'Program ORBITSFINAL is done.'
;
IF EMESS NE '' THEN GOTO, HANDLE_ERROR
IF N_ELEMENTS(ERRMSG) GT 0 THEN ERRMSG=EMESS
RETURN
;
; Error handling portion of the procedure.
;
HANDLE_ERROR:
IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = EMESS ELSE PRINT, EMESS
;
RETURN
END