;+
; Project : Computational Physics
;
; Name : ORBEQNHW
;
; Purpose : Supplies ORBITSHW and RK4 with derivative equations.
;
; Explanation : This function calculates the derivatives for distance and
; velocity of a planet in orbit about the Sun. This function is
; designed to be called from ORBITSHW and RK4 (which is
; called by ORBITSHW).
;
; Use : DYDT = ORBEQNHW(T, Y)
;
; Inputs : T: The time (independent variable) where solution is
; sought.
;
; Y: The dependent variables which are used to calculate the
; derivatives.
;
; Opt. Inputs : None.
;
; Outputs : DYDT: The calculated derivatives.
;
; Opt. Outputs: None.
;
; Keywords : ERRMSG: If defined and passed, then any error messages will
; be returned to the user in this parameter rather than
; being printed internbally by this procedure. 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 = ''
; DYDT = ORBEQNHW(T, Y, ERRMSG=ERRMSG)
; IF ERRMSG(0) NE '' THEN ...
;
; Calls : None.
;
; 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 December 2002
; Initial program.
; Version 2, Donald G. Luttermoser, ETSU/Physics, 3 May 2007
; Replaced use of the IDL MESSAGE utility with a simple print
; statement.
; Version 3, Donald G. Luttermoser, ETSU/Physics, 23 Nov 2009
; Changed function name from ORBEQN to ORBEQNFINAL 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 4, Donald G. Luttermoser, ETSU/Physics, 8 Nov 2011
; Change function name to ORBEQNHW due to modifications
; made for using this as a homework problem.
;
; Version : Version 4, 8 November 2011.
;
;-
;
FUNCTION ORBEQNHW, T, Y, ERRMSG=ERRMSG
;
EMESS = '' ; Set to non-null string if error is encountered.
;
NEQN = N_ELEMENTS(Y)
DYDT = DBLARR(NEQN)
;
IF NEQN NE 4 THEN BEGIN
EMESS = 'Number of equations (i.e., Y variable) must be 4!'
GOTO, HANDLE_ERROR
ENDIF
;
RCART = [Y(0), Y(1)] ; [x(AU), y(AU)]
VCART = [Y(2), Y(3)] ; [v_x(AU/yr), v_y(AU/yr)]
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.
;
; Students: Insert an equation for r cubed (RFORM3) as a function of
; x and y in AU:
;
RNORM3 =
;
; Students: Insert an equation for the acceleration (ACCEL) stored in
; the following 2-element array (remember, RCART is a 2-element array
; containing x and y) in units of AY/yr.
;
ACCEL = ; ACCEL in AY/yr^2.
;
; DYDT[0] = dx/dt = v(x); DYDT[1] = dy/dt = v(y);
; DYDT[2] = dv(x)/dt = a(x); DYDT[3] = dv(y)/dt = a(y).
;
DYDT = [VCART, ACCEL]
;
IF EMESS NE '' THEN GOTO, HANDLE_ERROR
IF N_ELEMENTS(ERRMSG) GT 0 THEN ERRMSG = EMESS
RETURN, DYDT
;
; Error handling portion of the procedure.
;
HANDLE_ERROR:
IF N_ELEMENTS(ERRMSG) NE 0 THEN ERRMSG = EMESS ELSE PRINT, EMESS
;
RETURN, DYDT
END