;+
; Project     : Course Computer Project
;
; Name        : POLYTROPE
;
; Purpose     : Driver procedure to generate polytrope models.
;
; Explanation : This program generates a set of polytropes of index N 
;               to represent interior models for stars and white dwarf
;               The models parameters are defined in the main driving 
;               procedure (this code).  A fourth-order Runge-Kutta
;               scheme is used to solve the differential equations used.
;
; Use         : POLYTROPE  (no parameters passed)
;
; Inputs      : None.
;
; Opt. Inputs : None.
;
; Outputs     : None.  Output is stored in various output files.
;
; 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 = ''
;                            POLYTROPE, ERRMSG=ERRMSG
;                            IF ERRMSG(0) NE '' THEN ...
;
; Calls       : PRINTOUT, RDPOLYMODEL, RK4COEFF, RK4STEP.
;
; 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 mypolytrope.pro procedure.
;
; Written     : Donald G. Luttermoser, ETSU/Physics,  5 Oct 2011.
;
; Modified    : Version 1, Donald G. Luttermoser, ETSU/Physics,  5 Oct 2011
;                   Initial program.
;
; Version     : Version 1,  5 October 2011.
;
;-
PRO POLYTROPE, ERRMSG=ERRMSG
;
;ON_ERROR, 2   ; Return to caller if error is encountered.
EMESS = ''     ; Set to non-null string if error is encountered.
AFONT = [-1, 0, 0]  &  PASS = 0
PCHARSIZE = [1.4, 1.3, 2.0]  &  PTHICK = [3, 4, 5] 
;
; Greek letters for plots.
;
SXI = ['!7n', '!7x', '!7x']     &  SRHO = ['!7q', '!7r', '!7r']
STHETA = ['!7h', '!7q', '!7q']  &  SPHI = ['!7u', '!7f', '!7f']
;
; Check the machine type and IDL version release number.
;
TERM = STRUPCASE(!D.NAME)
CASE TERM OF
    'X': BEGIN
         OSTYPE = 'Unix'  &  DELIM = '/'
         END
    'WIN': BEGIN
         OSTYPE = 'Windows'  &  DELIM = '\'
         END
    ELSE: BEGIN
          EMESS = 'Unknown windows and machine type: ' + TERM
          GOTO, HANDLE_ERROR
          END
ENDCASE
DPOS = STRPOS(!VERSION.RELEASE, '.')
IDLVERS = FLOAT(STRMID(!VERSION.RELEASE, 0, DPOS+2))
;
; Input parameters need for the model.
;
POLY_INDEX_ARRAY = [1.0, 1.5, 3.0]  ; Regular Undergraduate Students.
POLY_INDEX_HONORS_ARRAY = [0., 0.5, 1.0, 1.5, 2.5, 3.0, 3.5, 4.0]  ; Honors Students
TITLEDEF = 'Student Polytrope Model'  ;  Default title.
HI = 0.01   ;  Initial step size of radius parameter.
OUTROOTNAME = 'polyout'  & OUTSUFFIXNAME = '.txt'
PSROOTNAME = 'ptrope'
NMODELDEPTHMAX = 0  ; This parameter will keep track of the maximum depth points computed.
;
; Set a maximum number of depths to insure that we don't integrate forever.
;
NDMAX = 30000L
;
; ********* NOTE!  Read the Instructions for details of the required 
; ********* additional coding you will be doing in this procedure.
; ********* 
; ********* First define the integer scalar NPOLY here (perhaps making use
; ********* of the N_ELEMENTS function of IDL.  Then define a string array
; ********* to store the names of your various output files (say MYOUTFILES).
; ********* 
; ********* At this position you need to set up a FOR loop to 
; ********* cycle through the POLY_INDEX arrays (POLY_INDEX_ARRAY for 
; ********* undergrads or POLY_INDEX_HONORS for honors undergrads).
; ********* Note that the first line after your FOR loop should associate
; ********* the parameter POLY_INDEX to the approriate POLY_INDEX array 
; ********* for the current index counter, e.g.,
; *********   POLY_INDEX = POLY_INDEX_ARRAY[J]
; ********* or
; *********   POLY_INDEX = POLY_INDEX_HONORS_ARRAY[J]
; ********* depending upon whether you are a regular undergraduate or an 
; ********* honors undergraduate/graduate student, where J is the loop counter.
; *********
; ********* Just after that line, you need to define a specific output 
; ********* file name for the given POLY_INDEX run, e.g.,
; *********   OUTFILE = OUTROOTNAME + STRING(J, FORMAT='(I1)') + OUTSUFFIXNAME
; ********* then store this name for future use in MYOUTFILES, e.g.,
; *********   MYOUTFILES[J] = OUTFILE
; *********
; ********* Finally at this point, you will need to create a specific title for
; ********* the current output file making use of the default title similar to
; *********   TITLE = TITLEDEF + ' # ' + STRING(J, FORMAT='(I1)')
;
; Define internal parameters.
;
IDEPTH = 0L  ;  Depth indicie, IDEPTH = 0 is the center of the polytrope.
IFINAL = 0L  ; If non-zero, the integration is complete.
H = HI       ; Current depth step size -- set to initial value here.
KCOEFF = DBLARR(4)   ;  RK4 K coefficients (see write up).
LCOEFF = DBLARR(4)   ;  RK4 L coefficients (see write up).
ALABL = ['IDEPTH', 'XI', 'PHI', 'THETA', 'XI/XI_s', 'RHO/RHOBAR']
;
; 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.
;
IFLAG = INTARR(3)
;
; Define initial values of output arrays.
;
XI = [0.D0]  &  PHI = [0.D0]  &  THETA = [1.D0]
;
; Open the output file and write file header.
;
PRINTOUT, OUTFILE, UWR, XI, PHI, THETA, POLY_INDEX, HI, IFINAL, IFLAG, $
    'head', TITLE, ALABL, ERRMSG=EMESS
IF EMESS[0] NE '' THEN GOTO, HANDLE_ERROR
;
; Set the step size H to the initial step size:
;
H = HI
;
; Carry out the integration using 4th-order Runge-Kutta.
;
WHILE IFINAL EQ 0 DO BEGIN
;
;  Calculate the Runge-Kutta coefficients.
;
    RK4COEFF, IDEPTH, IFLAG, H, XI, PHI, THETA, POLY_INDEX, KCOEFF, $
        LCOEFF, ERRMSG=EMESS
    IF EMESS NE '' THEN BEGIN
        FREE_LUN, UWR
        GOTO, HANDLE_ERROR
    ENDIF
;
;  Calculate the next step in the Runge-Kutta integration.
;
    RK4STEP, IDEPTH, IFLAG, ALABL, NDMAX, IFINAL, H, XI, PHI, $
        THETA, POLY_INDEX, KCOEFF, LCOEFF, ERRMSG=EMESS
    IF EMESS NE '' THEN BEGIN
        PRINTOUT, OUTFILE, UWR, XI, PHI, THETA, POLY_INDEX, H, $
            IFINAL, IFLAG, 'error', TITLE, ALABL, ERRMSG=EMESS
        FREE_LUN, UWR
        GOTO, HANDLE_ERROR
    ENDIF
ENDWHILE
;
; Print the model in the output file and close it.
;
PRINTOUT, OUTFILE, UWR, XI, PHI, THETA, POLY_INDEX, H, IFINAL, $
    IFLAG, 'done', TITLE, ALABL, ERRMSG=EMESS
;
; Keep track of the maximum number of depth points in these models.
;
NMODELDEPTHMAX = MAX([IFINAL+1, NMODELDEPTHMAX])
;
; ********* It might be a good idea to put a little print statement here letting the
; ********* user know that the current model has been computed.
;
; ********* You need to close out the FOR\DO block at this point:
; *********   ENDFOR
;
;
; ********* Now we create the needed plots. (Did you define NPOLY?).
;
XIARRAY = DBLARR(NPOLY, NMODELDEPTHMAX)  &  PHIARRAY = XIARRAY
THETAARRAY = XIARRAY  &  RHO_RATIOARRAY = XIARRAY
XI2XI_SARRAY = XIARRAY  &  CEN2AVEDENS = DBLARR(NPOLY)
MFCONST = CEN2AVEDENS  &  MRCONST = MFCONST  &  KPARAM = MFCONST
PCONST = MFCONST  &  NDEPTH = LONARR(NPOLY)
;
; Read in and save the data from the output files for plotting.
;
FOR J = 0, NPOLY-1 DO BEGIN
    RDPOLYMODEL, MYOUTFILES[J], FTITLE, FPOLY_INDEX, FHI, FXI, FPHI, FTHETA, $
        FXI2XI_S, FRHO_RATIO, FCEN2AVEDENS, FMFCONST, FMRCONST, FKPARAM, $
        FPCONST, ERRMSG=EMESS
    IF EMESS[0] NE '' THEN GOTO, HANDLE_ERROR
;
    NDEPTH[J] = N_ELEMENTS(FXI)
    IF N_ELEMENTS(FCEN2AVEDENS) GT 0 THEN CEN2AVEDENS[J] = FCEN2AVEDENS
    IF N_ELEMENTS(FMFCONST) GT 0 THEN MFCONST[J] = FMFCONST
    IF N_ELEMENTS(FMRCONST) GT 0 THEN MRCONST[J] = FMRCONST
    IF N_ELEMENTS(FKPARAM) GT 0 THEN KPARAM[J] = FKPARAM
    IF N_ELEMENTS(FPCONST) GT 0 THEN PCONST[J] = FPCONST
    XIARRAY[J, 0:NDEPTH[J]-1] = FXI
    PHIARRAY[J, 0:NDEPTH[J]-1] = FPHI
    THETAARRAY[J, 0:NDEPTH[J]-1] = FTHETA
    XI2XI_SARRAY[J, 0:NDEPTH[J]-1] = FXI2XI_S
    RHO_RATIOARRAY[J, 0:NDEPTH[J]-1] = FRHO_RATIO
ENDFOR
;
; Determine minima and maxima for these arrays.
;
XIMAX = MAX(XIARRAY)  &  PHIMIN = MIN(PHIARRAY)
RHORATMAX = MAX(RHO_RATIOARRAY)
;
; Make the first plot.
;
PLOTIT1:
;
PLOT, XIARRAY[0, 0:NDEPTH[0]-1], THETAARRAY[0, 0:NDEPTH[0]-1], $
    XTITLE=SXI[PASS], YTITLE=STHETA[PASS], FONT=AFONT[PASS], XSTYLE=1, $
    YSTYLE=1, CHARSIZE=PCHARSIZE[PASS], XRANGE=[0,XIMAX+1.], $
    YRANGE=[0,1.1], THICK=PTHICK[PASS]
FOR JJ = 1, NPOLY-1 DO BEGIN
    OPLOT, XIARRAY[JJ, 0:NDEPTH[JJ]-1], THETAARRAY[JJ, 0:NDEPTH[JJ]-1], $
        THICK=PTHICK[PASS]
ENDFOR
WAIT, 0.2  ; Let the operating system finish drawing the graph.
;
; ********* Can you think of a way to label each curve on the plot?
; ********* Note that if you use the CURSOR command that you only do this 
; ********* if PASS = 0.
;
;
; Save to hardcopy file.  First make an encapsulated PS plot then a standard PS plot.
;
CASE PASS OF
    0:  BEGIN             ; Encapsulated postscript file.
        SET_PLOT, 'PS'
        DEVICE, /PORTRAIT, /ENCAPSULATE, FILENAME=PSROOTNAME+'1.eps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 1
        GOTO, PLOTIT1
        END
    1:  BEGIN              ; Normal postscript file.
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Encapsulated postscript file stored in: "'+ $
            PSROOTNAME+'1.eps"'
        PRINT, '!!! Use this file in your LaTeX manuscript.'
        PRINT, ' '
        SET_PLOT, 'PS'
        DEVICE, /LANDSCAPE, FILENAME=PSROOTNAME+'1.ps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 2
        GOTO, PLOTIT1
        END
    2:  BEGIN
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Normal postscript file stored in: "'+ $
            PSROOTNAME+'1.ps"'
        PRINT, '!!! You can send this plot directly to a printer.'
        PRINT, ' '
        SET_PLOT, TERM
        PASS = 0
        END
ENDCASE
;
; Make the second plot.
;
PLOTIT2:
;
PLOT, XI2XI_SARRAY[0, 0:NDEPTH[0]-1], THETAARRAY[0, 0:NDEPTH[0]-1], $
    XTITLE=SXI[PASS]+'!6/'+SXI[PASS]+'!6!Ds!N', YTITLE=STHETA[PASS], $
    FONT=AFONT[PASS], CHARSIZE=PCHARSIZE[PASS], THICK=PTHICK[PASS], $
    XSTYLE=1, YSTYLE=1, XRANGE=[0,1.1], YRANGE=[0,1.1]
FOR JJ = 1, NPOLY-1 DO BEGIN
    OPLOT, XI2XI_SARRAY[JJ, 0:NDEPTH[JJ]-1], $
        THETAARRAY[JJ, 0:NDEPTH[JJ]-1], THICK=PTHICK[PASS]
ENDFOR
WAIT, 0.2  ; Let the operating system finish drawing the graph.
;
; ********* Can you think of a way to label each curve on the plot?
; ********* Note that if you use the CURSOR command that you only do this 
; ********* if PASS = 0.
;
;
; Save to hardcopy file.  First make an encapsulated PS plot then a standard PS plot.
;
CASE PASS OF
    0:  BEGIN             ; Encapsulated postscript file.
        SET_PLOT, 'PS'
        DEVICE, /PORTRAIT, /ENCAPSULATE, FILENAME=PSROOTNAME+'2.eps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 1
        GOTO, PLOTIT2
        END
    1:  BEGIN              ; Normal postscript file.
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Encapsulated postscript file stored in: "'+ $
            PSROOTNAME+'2.eps"'
        PRINT, '!!! Use this file in your LaTeX manuscript.'
        PRINT, ' '
        SET_PLOT, 'PS'
        DEVICE, /LANDSCAPE, FILENAME=PSROOTNAME+'2.ps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 2
        GOTO, PLOTIT2
        END
    2:  BEGIN
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Normal postscript file stored in: "'+ $
            PSROOTNAME+'2.ps"'
        PRINT, '!!! You can send this plot directly to a printer.'
        PRINT, ' '
        SET_PLOT, TERM
        PASS = 0
        END
ENDCASE
;
PLOTIT3:
;
PLOT, XIARRAY[0, 0:NDEPTH[0]-1], PHIARRAY[0, 0:NDEPTH[0]-1], $
    XTITLE=SXI[PASS], YTITLE=SPHI[PASS], FONT=AFONT[PASS], $
    CHARSIZE=PCHARSIZE[PASS], THICK=PTHICK[PASS], XSTYLE=1, YSTYLE=1, $
    XRANGE=[0,XIMAX+1.], YRANGE=[PHIMIN-0.1,0.1]
FOR JJ = 1, NPOLY-1 DO BEGIN
    OPLOT, XIARRAY[JJ, 0:NDEPTH[JJ]-1], PHIARRAY[JJ, 0:NDEPTH[JJ]-1], $
        THICK=PTHICK[PASS]
ENDFOR
WAIT, 0.2  ; Let the operating system finish drawing the graph.
;
; ********* Can you think of a way to label each curve on the plot?
; ********* Note that if you use the CURSOR command that you only do this 
; ********* if PASS = 0.
;
;
; Save to hardcopy file.  First make an encapsulated PS plot then a standard PS plot.
;
CASE PASS OF
    0:  BEGIN             ; Encapsulated postscript file.
        SET_PLOT, 'PS'
        DEVICE, /PORTRAIT, /ENCAPSULATE, FILENAME=PSROOTNAME+'3.eps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 1
        GOTO, PLOTIT3
        END
    1:  BEGIN              ; Normal postscript file.
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Encapsulated postscript file stored in: "'+ $
            PSROOTNAME+'3.eps"'
        PRINT, '!!! Use this file in your LaTeX manuscript.'
        PRINT, ' '
        SET_PLOT, 'PS'
        DEVICE, /LANDSCAPE, FILENAME=PSROOTNAME+'3.ps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 2
        GOTO, PLOTIT3
        END
    2:  BEGIN
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Normal postscript file stored in: "'+ $
            PSROOTNAME+'3.ps"'
        PRINT, '!!! You can send this plot directly to a printer.'
        PRINT, ' '
        SET_PLOT, TERM
        PASS = 0
        END
ENDCASE
;
PLOTIT4:
;
PLOT, XIARRAY[0, 0:NDEPTH[0]-1], RHO_RATIOARRAY[0, 0:NDEPTH[0]-1], $
    XTITLE=SXI[PASS], YTITLE=SRHO[PASS]+'!6/'+SRHO[PASS]+'!6!Dave!N', $
    FONT=AFONT[PASS], CHARSIZE=PCHARSIZE[PASS], THICK=PTHICK[PASS], $
    XSTYLE=1, YSTYLE=1, XRANGE=[0,XIMAX+1.], YRANGE=[0,RHORATMAX+2]
FOR JJ = 1, NPOLY-1 DO BEGIN
    OPLOT, XIARRAY[JJ, 0:NDEPTH[JJ]-1], $
        RHO_RATIOARRAY[JJ, 0:NDEPTH[JJ]-1], THICK=PTHICK[PASS]
ENDFOR
WAIT, 0.2  ; Let the operating system finish drawing the graph.
;
; ********* Can you think of a way to label each curve on the plot?
; ********* Note that if you use the CURSOR command that you only do this 
; ********* if PASS = 0.
;
;
; Save to hardcopy file.  First make an encapsulated PS plot then a standard PS plot.
;
CASE PASS OF
    0:  BEGIN             ; Encapsulated postscript file.
        SET_PLOT, 'PS'
        DEVICE, /PORTRAIT, /ENCAPSULATE, FILENAME=PSROOTNAME+'4.eps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 1
        GOTO, PLOTIT4
        END
    1:  BEGIN              ; Normal postscript file.
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Encapsulated postscript file stored in: "'+ $
            PSROOTNAME+'4.eps"'
        PRINT, '!!! Use this file in your LaTeX manuscript.'
        PRINT, ' '
        SET_PLOT, 'PS'
        DEVICE, /LANDSCAPE, FILENAME=PSROOTNAME+'4.ps'
        DEVICE, /SCHOOLBOOK, /BOLD, FONT_INDEX=17
        DEVICE, /SYMBOL, FONT_INDEX=7
        DEVICE, /BKMAN, /DEMI, FONT_INDEX=6
        DEVICE, /HELVETICA, /BOLD, FONT_INDEX=5
        PASS = 2
        GOTO, PLOTIT4
        END
    2:  BEGIN
        DEVICE, /CLOSE_FILE
        PRINT, ' '
        PRINT, '!!! Normal postscript file stored in: "'+ $
            PSROOTNAME+'4.ps"'
        PRINT, '!!! You can send this plot directly to a printer.'
        PRINT, ' '
        SET_PLOT, TERM
        PASS = 0
        END
ENDCASE
;
; Label the curves.  (To be inserted later.)
;
;
; Let the user know that the program is complete.
;
PRINT, ' '
PRINT, '**********************************************'
PRINT, 'POLYTROPE procedure has successfully finished!'
PRINT, '**********************************************'
;
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