;+ ; Project : Course Computer Project ; ; Name : RK4STEP ; ; Purpose : Set up polytrope 4th order Runge-Kutta coefficients. ; ; Explanation : This procedure calculates values for the next step in ; the 4th-order Runge-Kutta integration. Results for every ; 10th step is printed to the screen. ; ; Use : RK4STEP, IDEPTH, IFLAG, ALABL, NDMAX, IFINAL, 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. ; ; ALABL: Column labels for output data. [STRING array] ; ; NDMAX: Maximum number of depths to prevent infinite loops. ; [LONG] ; IFINAL: The depth of the surface. Note that IFINAL = 0 ; ; until the surface is reached. [LONG] ; ; 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] ; ; 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. Inputs : None. ; ; Outputs : 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] ; ; 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, ALABL, NDMAX, IFINAL, $ ; 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 RK4STEP 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 RK4STEP, IDEPTH, IFLAG, ALABL, NDMAX, IFINAL, H, XI, PHI, $ THETA, POLY_INDEX, KCOEFF, LCOEFF, ERRMSG=ERRMSG ; EMESS = '' ; XI = [XI, XI[IDEPTH] + H] PHI = [PHI, PHI[IDEPTH] + (KCOEFF[0] + $ 2.0D0*(KCOEFF[1]+KCOEFF[2]) + KCOEFF[3]) / 6.0D0] THETA = [THETA, THETA[IDEPTH] + (LCOEFF[0] + $ 2.0D0*(LCOEFF[1]+LCOEFF[2]) + LCOEFF[3]) / 6.0D0] ; ; Print to screen. Uncomment these for debugging. ; ;IF IDEPTH EQ 1 THEN BEGIN ; PRINT, ALABL[0], ALABL[1], ALABL[2], ALABL[3], $ ; FORMAT='(/A5, 9X, A7, 8X, A7, 7X, A7/)' ;ENDIF ; ;IF (IDEPTH MOD 10) EQ 0 THEN BEGIN ; PRINT, IDEPTH, XI[IDEPTH], PHI[IDEPTH], THETA[IDEPTH], $ ; FORMAT='(I5,3F15.7)' ;ENDIF ; NEXT: IDEPTH = IDEPTH + 1 ; STATEMENT 100 ; ; Compare Runge-Kutta value for THETA with series expansion when XI ; is small. If discrepancy is greater than a factor of 2, ; integration is stopped and discrepancy and step position are ; printed. ; IF XI[IDEPTH] LT 0.1D0 THEN BEGIN T = 1.0D0 - (XI[IDEPTH] * XI[IDEPTH]) / 6.0D0 + $ (XI[IDEPTH]^4) / 120.0D0 RATIO_T = T / THETA[IDEPTH] IF (RATIO_T LT 0.5D0) OR (RATIO_T GT 2.0D0) THEN BEGIN PRINT, 'Series approximation does not agree when XI << 1.' IFLAG[0] = IFLAG[0] + 1 EMESS = ['Series approximation does not agree when XI << 1.', $ STRING(IDEPTH, FORMAT='("IDEPTH = ", I3)')+STRING(XI[IDEPTH], $ FORMAT='(", XI = ", F12.6)')+STRING(RATIO_T, $ FORMAT='(", RATIO = ", F12.6)'), $ 'Program is aborting in RK4STEP.'] GOTO, HANDLE_ERROR ENDIF ENDIF ; ; Surface check: This condition is met if |THETA| < 1.0E10. When ; the surface is found, the polytrope model is written to the output ; file upon exit from this procedure and the program successfully ; ends. ; IF ABS(THETA[IDEPTH]) LT 1.0D-10 THEN BEGIN IFINAL = IDEPTH - 1 PRINT, XI[IFINAL], FORMAT='(/" XI(surface) = ", F15.7)' RETURN ENDIF ; ; Stop the integration if IDEPTH >= NDMAX. ; IF IDEPTH GE NDMAX THEN BEGIN EMESS = 'IDEPTH >= NDMAX! Increase step size and rerun.' GOTO, HANDLE_ERROR RETURN ENDIF ; ; Check to see if surface radius is becoming infinitely large. If so, ; stop the integration, write the model to the output file and let the ; user know the model is infinitely large. ; SLOPE = ABS(PHI[IDEPTH]) IF (SLOPE LE 1.0D-5) AND (THETA[IDEPTH] LE 5.0D-3) THEN BEGIN IFLAG[2] = 50 IFINAL = IDEPTH - 1 PRINT, 'Polytrope surface is becoming infinitely large.' PRINT, SLOPE, XI[IFINAL], FORMAT='(" SLOPE = ", F9.5, " at XI = ", F10.6)' ; 150 FORMAT(' SLOPE = ', F9.5, ' at XI = ', F10.6) RETURN ENDIF ; ; Check to see if the surface has been passed (i.e., THETA < 0). ; If so, halve the step size and go back to the previous step and ; integrate again. Make IDEPTH stay a positive integer. ; IF THETA[IDEPTH] LT 0.0D0 THEN BEGIN H = H / 2.0D0 IDEPTH = (IDEPTH - 1L) > 0L XI = XI[0:IDEPTH] & THETA = THETA[0:IDEPTH] & PHI = PHI[0:IDEPTH] XI[IDEPTH] = XI(IDEPTH) - H IFLAG2 = IFLAG2 + 1 ENDIF ; 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) EQ 0 THEN PRINT, EMESS ERRMSG = EMESS RETURN END