;+
; Project     : Course Computer Project
;
; Name        : RK4COEFF
;
; Purpose     : Set up polytrope 4th order Runge-Kutta coefficients.
;
; Explanation : This procedure calculates the 4th-order Runge-Kutta 
;               integration coefficients.  If XI=0, PHI/XI=-1/3, and 
;               THETA=1, each argument that is raised to the POLY_INDEX
;               power is checked to determine if the argument is positive 
;               or negative.  If the argument is negative, the step size
;               H is halved and the coefficients are recalculated.
;
; Use         : RK4COEFF, IDEPTH, IFLAG, H, XI, PHI, THETA, POLY_INDEX, $
;                   KCOEFF, LCOEFF
;
; Inputs      : IDEPTH:  Depth indicie, IDEPTH = 0 is the center of the 
;                        polytrope.  [LONG]
;
;               IFLAG:   Integration flags:  [3-element INTEGER array]
;
;                        IFLAG[0]:  If > 1 indicates discrepancy in THETA 
;                                   when XI << 1.
;                        IFLAG[1]:  Counts the number of times H has been 
;                                   halved.
;                        IFLAG[2]:  A value of 50 indicates the model radius 
;                                   is approaching infinity.
;
;               H:       The depth step size in units of XI.  [FLOAT]
;
;               XI:      The dimensionless depth scale factor of the Fortran
;                        run.  [DOUBLE array]
;
;               PHI:     The derivative dTHETA/dXI.  [DOUBLE array]
;
;               THETA:   Dimensionless scale factor of density (RHO) defined 
;                        by RHO = RHO_CENTRAL * THETA^POLY_INDEX.  
;                        [DOUBLE array]
;
;               POLY_INDEX:  The index of the polytrope in the Fortran run, 
;                        where 0 <= POLY_INDEX < 5.  Note that POLY_INDEX >= 
;                        5 produces a polytrope with infinite radius.  [FLOAT]
;
; Opt. Inputs : None.
;
; Outputs     : KCOEFF:  RK4 integration K coefficients.  See write-up for 
;                        for details.  [4-element DOUBLE array]
;
;               LCOEFF:  RK4 integration L coefficients.  See write-up for 
;                        for details.  [4-element DOUBLE array]
;
; 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 internally 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 = ''
;                            RK4COEFF, IDEPTH, IFLAG, H, XI, PHI, THETA, $
;                                POLY_INDEX, KCOEFF, LCOEFF, ERRMSG=ERRMSG
;                            IF ERRMSG(0) NE '' THEN ...
;
; Calls       : None.
;
; Common      : None.
;
; Restrictions: IDL Version 5.0 or later must be used for the polytrope procedures.
;
; Side effects: None.
;
; Category    : Class project, stellar modeling.
;
; Prev. Hist. : Based on subroutine RK4COEFF in ptrope.f Fortran code.
;
; Written     : Donald G. Luttermoser, ETSU/Physics,  4 Aug 2011.
;
; Modified    : Version 1, Donald G. Luttermoser, ETSU/Physics,  4 Aug 2011
;                   Initial program.
;
; Version     : Version 1,  4 August 2011.
;
;-
PRO RK4COEFF, IDEPTH, IFLAG, H, XI, PHI, THETA, POLY_INDEX, KCOEFF, $
    LCOEFF, ERRMSG=ERRMSG
;
STEP1:
IF XI[IDEPTH] EQ 0.D0 THEN KCOEFF[0] = -H/3.0D0 ELSE $
    KCOEFF[0] = -H*((2*PHI[IDEPTH]/XI[IDEPTH]) + $
        (THETA[IDEPTH]^POLY_INDEX))
LCOEFF[0] = H*PHI[IDEPTH]
;
THETAN = THETA[IDEPTH] + LCOEFF[0]/2.0D0
IF THETAN LT 0.0D0 THEN BEGIN
    H = H/2.0D0
    GOTO, STEP1
ENDIF ELSE KCOEFF[1] = -H*(2.0D0*(PHI[IDEPTH] + KCOEFF[0]/2.0D0) $
    / (XI[IDEPTH] + H/2.0D0) + (THETAN^POLY_INDEX))
LCOEFF[1] = H*(PHI[IDEPTH] + KCOEFF[0]/2.0D0)
;
THETAN = THETA[IDEPTH] + LCOEFF[1]/2.0D0
IF THETAN LT 0.0D0 THEN BEGIN
    H = H/2.0D0
    GOTO, STEP1
ENDIF ELSE KCOEFF[2] = -H*(2.0D0*(PHI[IDEPTH] + KCOEFF[1]/2.0D0) / $ 
    (XI[IDEPTH] + H/2.0D0) + (THETAN^POLY_INDEX))
LCOEFF[2] = H*(PHI[IDEPTH] + KCOEFF[1]/2.0D0)
;
THETAN = THETA[IDEPTH] + LCOEFF[2]
IF THETAN LT 0.0D0 THEN BEGIN
    H = H/2.0D0
    GOTO, STEP1
ENDIF ELSE KCOEFF[3] = -H*(2.0D0*(PHI[IDEPTH] + KCOEFF[2]) / $ 
    (XI[IDEPTH] + H) + (THETAN^POLY_INDEX))
LCOEFF[3] = H*(PHI[IDEPTH] + KCOEFF[2])
;
RETURN
END