;+
; Project     : Course Computer Project
;
; Name        : PRINTOUT
;
; Purpose     : Printout the polytrope model.
;
; Explanation : This procedure produces the output file with the run data.
;
; Use         : PRINTOUT, OUTFILE, UWR, XI, PHI, THETA, POLY_INDEX, H, $
;                   IFINAL, IFLAG, MODE, TITLE, ALABL
;
; Inputs      : OUTFILE: The name of the output file.  [STRING]
;
;               UWR:     The unit number of the output file.  [INTEGER]
;
;               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]
;
;               H:       The depth step size in units of XI.  [FLOAT]
;
;               IFINAL:  The index of the outermost layer 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.
;
;               MODE:    Let PRINTOUT know what should be done.  [STRING]
;
;               TITLE:   The main title to be placed in the output file.  
;                        [STRING]
;
;               ALABL:   Column labels for output data.  [STRING array]
;
; Opt. Inputs : None.
;
; Outputs     : None.
;
; 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 = ''
;                            PRINTOUT, OUTFILE, UWR, XI, PHI, THETA, $
;                                POLY_INDEX, H, IFINAL, IFLAG, MODE, $
;                                TITLE, ALABL, 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, thermodynamic modeling.
;
; Prev. Hist. : Based on subroutine PRINTOUT 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 PRINTOUT, OUTFILE, UWR, XI, PHI, THETA, POLY_INDEX, H, IFINAL, IFLAG, $
    MODE, TITLE, ALABL, ERRMSG=ERRMSG
;
EMESS = ''
;
; DE:  Ratio of density at XI to that of the average density.
; FR:  Ratio XI(IDEPTH) / XI(SURFACE).
;
;      PARAMETER (DPI4 = 12.56637061436D0)
;
CASE MODE OF
    'head': BEGIN
        OPENW, UWR, OUTFILE, ERROR=IERR, /GET_LUN
        IF IERR NE 0 THEN BEGIN
            EMESS = ['Failed to open output file: '+OUTFILE, '  '+!ERROR_STAT.MSG]
            GOTO, HANDLE_ERROR
        ENDIF
        PRINT, 'Opening file: ', OUTFILE
;
        PRINTF, UWR, TITLE, FORMAT= $
            '(/2X, "Output data from a polytrope.pro run.", /2X, "TITLE:  ", A/)'
        PRINTF, UWR, POLY_INDEX, FORMAT='(5X,"Polytrope model for INDEX: ", F6.2/)'
        PRINTF, UWR, H, FORMAT='(5X,"Initial step size: ", F6.3/)'
        PRINTF, UWR, ALABL[0], ALABL[1], ALABL[2], ALABL[3], ALABL[4], ALABL[5], $
            FORMAT='(/1X, A6, 9X, A2, 12X, A3, 11X, A5, 9X, A7, 6X, A10/)
        END
;
; Print out the error from RK4STEP.
;
    'error': BEGIN
        NERR = N_ELEMENTS(ERRMSG)
        IF NERR GT 0 THEN BEGIN
            FOR I = 0, NERR-1 DO PRINTF, UWR, ERRMSG[I]
            WAIT, 1.0  ; Make sure all contents were written.
            FREE_LUN, UWR
            RETURN
        ENDIF
        END
;
; Close the output file and report success.
;
    'done': BEGIN
;
; Print out the final model with the Lane-Emden constants.
;
        IF PHI[IFINAL] EQ 0.0D0 THEN BEGIN
            PRINT, 'PHI[IFINAL] = 0!  IFINAL = ', IFINAL
            FREE_LUN, UWR
            EMESS = 'PHI(IFINAL) = 0!  Stopping the code.'
            GOTO, HANDLE_ERROR
        ENDIF
        DENSRAT = -XI[IFINAL]/(3.0D0*PHI[IFINAL])
        PMASS = -XI[IFINAL]*XI[IFINAL]*PHI[IFINAL]
;
        IF XI[IFINAL] EQ 0.0D0 THEN BEGIN
            PRINT, 'XI[IFINAL] = 0!  IFINAL = ', IFINAL
            FREE_LUN, UWR
            EMESS = 'XI[IFINAL] = 0!  Stopping the code.'
            GOTO, HANDLE_ERROR
        ENDIF
        FR = DBLARR(IFINAL+1)  &  DE = FR
        FOR JJ = 0, IFINAL DO BEGIN
            FR[JJ] = XI[JJ] / XI[IFINAL]
            DE[JJ] = DENSRAT * THETA[JJ]^POLY_INDEX
        ENDFOR
        FOR JJ = 0, IFINAL DO PRINTF, UWR, JJ, XI[JJ], PHI[JJ], THETA[JJ], $
            FR[JJ], DE[JJ], FORMAT='(1X,I5,5F15.7)'
;
        PRINTF, UWR, ' '
        IF IFLAG[2] EQ 50 THEN BEGIN
            PRINTF, UWR, '        The final values calculated are:'
        ENDIF ELSE BEGIN
            PRINTF, UWR, '        The surface values are:'
        ENDELSE
        PRINTF, UWR, IFINAL, XI[IFINAL], PHI[IFINAL], THETA[IFINAL], $  
            FR[IFINAL], DE[IFINAL], FORMAT='(1X,I5,5F15.7)'
;
        PRINTF, UWR, ' '
        PRINTF, UWR, ' '
        PRINTF, UWR, '            The Lane-Emden constants are:'
;
        PRINTF, UWR, ' '
        IF IFLAG[2] EQ 50 THEN BEGIN
            PRINTF, UWR, '  Surface radius approaching infinite size!'
            PRINTF, UWR, ' '
            PRINTF, UWR, '  Mass-Radius constant = 0, all other constants '+ $
                'go to infinity.'
            PRINTF, UWR, ' '
        ENDIF ELSE BEGIN
            PRINTF, UWR, DENSRAT, $
                FORMAT='(8X,"Central Density/Average Density: ", E13.6)'
            PRINTF, UWR, PMASS, $
                FORMAT='(/8X,"Mass Function constant: ", F13.6)'
            IF POLY_INDEX EQ 1.0D0 THEN BEGIN
                PRINTF, UWR, ' '
                PRINTF, UWR, '        Mass-Radius relation constant undefined!'
                CK = ((!DPI*4/(XI[IFINAL]^(POLY_INDEX+1.0D0)* $
                    (-PHI[IFINAL])^(POLY_INDEX-1.0D0))))/(POLY_INDEX+1.0D0)
                WN = 1.0D0/(!DPI*4*(POLY_INDEX+1.0D0)*PHI[IFINAL]*PHI[IFINAL])
                PRINTF, UWR, CK, WN, $
                    FORMAT='(/8X,"K-parameter constant: ", F13.6, //8X,'+ $
                    '"Pressure constant: ", E13.6/)'
            ENDIF ELSE BEGIN
                EXN = (POLY_INDEX+1.0D0)/(POLY_INDEX-1.0D0)
                PMASRAD = -(XI[IFINAL]^EXN) * PHI[IFINAL]
                PRINTF, UWR, PMASRAD, $
                    FORMAT='(/8X,"Mass-Radius relation constant: ", F13.6)'
                IF POLY_INDEX EQ 0.0D0 THEN BEGIN
                    PRINTF, UWR, ' '
                    PRINTF, UWR, '        K-Parameter constant undefined!'
                ENDIF ELSE BEGIN
                    EXNN = 1.0D0/POLY_INDEX
                    CK = ((!DPI*4/(XI[IFINAL]^(POLY_INDEX+1.0D0)* $
                        (-PHI[IFINAL])^(POLY_INDEX-1.0D0)))^EXNN)/ $
                        (POLY_INDEX+1.0D0)
                    PRINTF, UWR, CK, $
                        FORMAT='(/8X,"K-parameter constant: ", F13.6)'
                ENDELSE
                WN = 1.0D0/(!DPI*4*(POLY_INDEX+1.0D0)*PHI[IFINAL]*PHI[IFINAL])
                PRINTF, UWR, WN, FORMAT='(/8X,"Pressure constant: ", E13.6/)'
            ENDELSE
            PRINTF, UWR, 'End of Data'
            WAIT, 0.2  ;  Give the operating system time to complete the write.
;
            FREE_LUN, UWR
            PRINT, 'Output stored in file: ', OUTFILE
        ENDELSE
        END
;
    ELSE:  ; Do nothing.
;
ENDCASE
;
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