;+
; Project     :	Computational Physics
;
; Name        :	ORBEQNFINAL
;
; Purpose     :	Supplies ORBITSFINAL 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 ORBITSFINAL and RK4 (which is
;               called by ORBITSFINAL).
;
; Use         :	DYDT = ORBEQNFINAL(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 = ORBEQNFINAL(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     :	Version 3,  23 November 2009.
;
;-
;
FUNCTION ORBEQNFINAL, 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