;+
; Project     :	Computational Physics
;
; Name        :	PROJ_RK4_EQN
;
; Purpose     :	Supplies PROJ_TUTORIAL and RK4 with derivative equations.
;
; Explanation :	This function calculates the derivatives for distance and
;               velocity of a projectile near the surface of a primary body.  
;               This function is designed to be called from PROJ_TUTORIAL and 
;               RK4 (which is called by PROJ_TUTORIAL).
;
; Use         :	DYDT = PROJ_RK4_EQN(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 = PROJ_RK4_EQN(T, Y, GRAV, ERRMSG=ERRMSG)
;                           IF ERRMSG[0] NE '' THEN ...
;
; Calls       :	None.
;
; Common      :	GRAVSURF.
;
; Restrictions:	None.
;
; Side effects:	None.
;
; Category    :	Computation physics.
;
; Prev. Hist. :	Based on ORBEQN.
;
; Written     :	Donald G. Luttermoser, ETSU/Physics, 9 November 2015
;
; Modified    :	Version 1, Donald G. Luttermoser, ETSU/Physics, 9 November 2015
;                   Initial program.
;
; Version     :	Version 1,  9 November 2015.
;
;-
;
FUNCTION PROJ_RK4_EQN, T, Y, ERRMSG=ERRMSG
;
COMMON GRAVSURF, GRAV
;
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
;
;       [x(m), y(m)]         [v_x(m/s), v_y(m/s)]
POSRK = [Y[0], Y[1]]  &  VELRK = [Y[2], Y[3]]
ACCEL = [0.D0, -GRAV]
;
; (0) dx/dt = v(x); (1) dy/dt = v(y); (2) dv(x)/dt = 0; (3) dv(y)/dt = -g
;
DYDT = [VELRK, 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