;+ ; Project : Computational Physics Homework Problem ; ; Name : ORBITSHW ; ; 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 = '' ; ORBITSHW, M, A, ECC, ERRMSG=ERRMSG ; IF ERRMSG(0) NE '' THEN ... ; ; Calls : ORBEQNHW (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 5, Donald G. Luttermoser, ETSU/Physics, 8 Nov 2011 ; Change procedure name to ORBITSHW due to modifications ; made for using this as a homework problem. ; ; Version : Version 5, 8 November 2011. ; ;- ; PRO ORBITSHW, 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: ORBITSHW [, 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 'ORBEQNHW'). ; ELSE: BEGIN STATE = [RCART, VCART] DYDT = ORBEQNHW(TIME[ISTEP], STATE, ERRMSG=EMESS) IF EMESS NE '' THEN GOTO, HANDLE_ERROR RESULT = RK4(STATE, DYDT, TIME[ISTEP], TAU, 'orbeqnhw', $/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 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 ORBITSHW 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