PRO marx_plots, marxdir, PARFILE = parfile, PSFILE = psfile
;;+
;
;
; PURPOSE: marx_plots will take as input a directory
;          containing marx output and will create
;          a series of postscript plots of the output.
;          Included among the plots are: 
;
;            - a plot of the simulation input spectrum
;            - a plot of pulse-height versus dispersion 
;              distance (1 plot per grating)
;            - a plot of the "full image", i.e. ypos versus zpos
;            - a plot of the first order (+/-1) spectrum
;              (1 plot per grating)
;
;
; INPUT:   marxdir  (string) : The directory containing the marx output
;
; OUTPUT: A series of postscript plots (names to be determined)        
;
; KEYWORD: parfile (string) : Set PARFILE = "parfile-name.par", if the
;          the parameter file used was NOT marx.par
;
; AUTHOR/DATE: C. Baluta, 5/5/97
;;-

TRUE = 1
FALSE = 0

ckev = 12.39854               ;;; NEED DEFAULT STRUCTURE VARIABLES HERE
angstrom = STRING("305B)

PRINT,'   '

IF ( N_PARAMS() NE 1) THEN BEGIN
 PRINT,'   '
 PRINT,'   marx_plots USAGE: '
 PRINT,'   '
 PRINT,'   marx_plots, marxdir, PARFILE = parfile'
 PRINT,'   '
 PRINT,'   marxdir (string) = directory containing marx output'
 PRINT,'   '
 PRINT,'   PARFILE (Keyword; string) = name of the parameter file'
 PRINT,'                               IF the file was NOT marx.par'
 PRINT,'   '
 RETURN
ENDIF

;; Get the parameter file name, concatenate the directory onto the name
IF ( NOT KEYWORD_SET(PARFILE) ) THEN parfile_name = 'marx.par'

put_slash, marxdir
parfile_name = marxdir + parfile_name
;;

;; Read in GratingType, DetectorType and other parameters
grating = read_parfile( 'GratingType', PARFILE = parfile_name)

IF (grating EQ '-1') THEN BEGIN
  PRINT,'  ' 
  PRINT,'  marx_plots.pro: Aborting.'
  PRINT,'  '
  RETURN
ENDIF

IF (grating EQ 'NONE') THEN BEGIN
  nograt = TRUE
ENDIF ELSE BEGIN
  nograt = FALSE
ENDELSE
 
;;; Can't find .par file. (See read_parfile.pro)
IF (STRCOMPRESS(STRING(grating),/REMOVE_ALL) EQ '-1') THEN RETURN 

detector = read_parfile( 'DetectorType', PARFILE = parfile_name)
mirror = read_parfile( 'MirrorType', PARFILE = parfile_name)
spec_type = read_parfile( 'SpectrumType', PARFILE = parfile_name)
numrays = read_parfile( 'NumRays', PARFILE = parfile_name)

IF (spec_type EQ 'FILE') THEN BEGIN
 spec_file = read_parfile( 'SpectrumFile', PARFILE = parfile_name)
 min_energy = 0.0
 max_energy = 0.0
 ray_flux = 0.0
ENDIF ELSE BEGIN
 spec_file = ''
 min_energy = read_parfile( 'MinEnergy', PARFILE = parfile_name)
 max_energy = read_parfile( 'MaxEnergy', PARFILE = parfile_name)
 ray_flux = read_parfile( 'RayFlux', PARFILE = parfile_name)
ENDELSE 

;;

PRINT,'    '
PRINT,'       Simulation Information:'
PRINT,'    '
PRINT,'                      Grating: ' + grating
PRINT,'                     Detector: ' + detector
PRINT,'                       Mirror: ' + mirror
PRINT,'    '
PRINT,'                     Spectrum: ' + spec_type
 IF (spec_type EQ 'FILE') THEN BEGIN
PRINT,'                Spectrum File: ' + spec_file
 ENDIF ELSE BEGIN
PRINT,'                  Min. Energy: ' + STRING(min_energy)
PRINT,'                  Max. Energy: ' + STRING(max_energy)
PRINT,'                         Flux: ' + STRING(ray_flux)
 ENDELSE
PRINT,'    '

;; Now read in or generate it (if SpectrumType='FLAT')
PRINT,'   '
PRINT,'   Reading in Input Spectrum now'
PRINT,'   '

get_input_spec, spec_type, spec_file, min_energy, max_energy, ray_flux, $
 Inputspec, Ffile

IF (Ffile EQ FALSE) THEN RETURN     ;; Ffile=Logical variable
                                    ;; =TRUE if file was found, FALSE otherwise

PRINT,'   Input spectrum read in/created successfully'
PRINT,'  '

;; Now read in the marx output vectors (ypos, zpos, time, order and pha)
read_marx_data, marxdir, ypos, zpos, time, order, pha, PARFILE = parfile_name


IF ( N_ELEMENTS(ypos) EQ 1 AND ypos(0) EQ -1) THEN BEGIN
  PRINT,'    '
  PRINT,'  Error in marx_plots.pro: No ypos data! '
  PRINT,'  Aborting.'
  PRINT,'    '
  RETURN
ENDIF
 
counts = N_ELEMENTS (ypos)

IF ( N_ELEMENTS(time) EQ 1 AND time(0) EQ -1) THEN BEGIN
  ttime = 1.0
ENDIF ELSE BEGIN
  ttime = MAX(time) - MIN(time)
ENDELSE


;;
IF (NOT KEYWORD_SET (PSFILE)) THEN BEGIN
 psfile = 'marx_plots.ps' 
ENDIF
  
SET_PLOT, 'ps'
DEVICE, FILENAME = psfile 
;;

;; First PLOTS: ypos versus zpos; PHA versus dispersion distance
!P.MULTI = [0, 1, 2]
xmax = 100.
xmin = -100.

tit = STRCOMPRESS(STRING(counts),/REMOVE_ALL) + ' cts; ' + $
 STRCOMPRESS(STRING(ttime),/REMOVE_ALL) + ' seconds' 
xtit = 'ypos (mm)'
ytit = 'zpos (mm)'

PLOT, ypos, zpos, PSYM = 3, TITLE = tit, XTITLE = xtit, YTITLE = ytit, $
 YRA = [xmin, xmax]  


;; Now create the first order spectrum
grating_spec, marxdir, Spectrum_first, GRATING = grating

!P.MULTI = [2, 2, 2]

IF (grating EQ 'LETG') THEN BEGIN
 
 cts = STRCOMPRESS(STRING(FIX(TOTAL(spectrum_first(1,*))) ), /REMOVE_ALL)
 tit = ' LEG +/- 1st Order Spectrum: ' + cts + ' cts'
 xtit = 'E (keV)'
 ytit = 'Counts'
 ymax = MAX(spectrum_first(1,*))

 PLOT_OO, spectrum_first(0,*), spectrum_first(1,*), XRA = [0.1,10], $
  XTIT = xtit, YTIT = ytit, YRA = [1, ymax], XSTY = 8, $
  YSTY = 8, XMARGIN = [8, 8], YMARGIN = [4, 4], SUBTITLE = tit

 AXIS, XAXIS = 1, XRANGE = ckev / (10. ^ !X.CRANGE), XSTYLE = 1, $
   XTITLE = 'Wavelength (' + angstrom +')' 

 AXIS, YAXIS = 1, YRANGE = (10. ^!Y.CRANGE) / ttime, YSTYLE =1, $
   YTITLE = 'Counts/s' 

;; Now create the PHA vs. dd plot
 y0 = 0.0
 z0 = 0.0
 maxdist = 0.2
 tg_extract, 'l', maxdist, ypos, zpos, y0, z0, dd_leg, dxd_leg, list_leg

 xx = WHERE(dd_leg GE 0.0)
 cts = STRCOMPRESS(STRING(N_ELEMENTS(dd_leg(xx))),/REMOVE_ALL)
 tit = 'LEG Pha vs. dd: ' + cts + ' Counts'
 
 PLOT, dd_leg, pha(list_leg), TITLE = tit, XTIT = 'dd_leg (mm)', $
   YTIT = 'Pha', XRA = [0.0, MAX(dd_leg)], PSYM = 3

ENDIF ELSE BEGIN            ;;;; HETG
   IF (grating EQ 'HETG') THEN BEGIN

      cts = FIX(TOTAL(spectrum_first(1,*)))
      cts_meg = STRCOMPRESS(STRING(cts),/REMOVE_ALL)

      cts = FIX(TOTAL(spectrum_first(3,*)))
      cts_heg = STRCOMPRESS(STRING(cts), /REMOVE_ALL)

      zero = WHERE(spectrum_first(1,*) EQ 0.0)
      spectrum_first(1,zero) = 1e-6

      zero = WHERE(spectrum_first(3,*) EQ 0.0)
      spectrum_first(3,zero) = 1e-6

      tit = 'MEG +/- 1st Order Spectrum: ' + cts_meg + ' cts'
      xtit = 'E (keV)'
      ytit = 'Counts'

      ymax = MAX( spectrum_first(1,*) )
      PLOT_OO, spectrum_first(0,*), spectrum_first(1,*), XRA=[0.1,10], $
        XTIT = xtit, YTIT = ytit, YRA = [1, ymax], XSTY = 8, YMARGIN = [4, 4], $
        XMARGIN = [8, 8], SUBTITLE = tit

      AXIS, XAXIS = 1, XRANGE = ckev / (10. ^ !X.CRANGE), XSTYLE = 1, $
        XTIT =  ' Wavelength (' + angstrom + ')' 
      AXIS, YAXIS = 1, YRANGE = (10. ^ !Y.CRANGE) / ttime, YSTYLE = 1, $
        YTIT = 'Counts/s'

      ymax = MAX( spectrum_first(3,*) )
      tit = 'HEG +/- 1st Order Spectrum: ' + cts_heg + ' cts'
      PLOT_OO, spectrum_first(2,*), spectrum_first(3,*), XRA=[0.1,10], $
        XTIT = xtit, YTIT = ytit, YRA = [1, ymax], $ 
        XSTY = 8, YSTY = 8, YMARGIN = [4, 4], XMARGIN = [8, 8], SUBTITLE = tit

      AXIS, XAXIS = 1, XRANGE = ckev / (10. ^ !X.CRANGE), XSTYLE = 1, $
        XTIT =  ' Wavelength (' + angstrom + ')'
      AXIS, YAXIS = 1, YRANGE = (10. ^ !Y.CRANGE) / ttime, YSTYLE = 1, $
        YTIT = 'Counts/s'

     ;; Next PLOT: PHA vs. dd
      !P.MULTI = [0, 2, 2]

      y0 = 0.0
      z0 = 0.0
      maxdist = 0.2
      tg_extract, 'h', maxdist, ypos, zpos, y0, z0, dd_heg, dxd_heg, list_heg
      tg_extract, 'm', maxdist, ypos, zpos, y0, z0, dd_meg, dxd_meg, list_meg
     
      xx = WHERE(dd_meg GE 0.0)
      cts = STRCOMPRESS(STRING(N_ELEMENTS(dd_meg(xx))),/REMOVE_ALL)

      tit = 'MEG Pha vs. dd: ' + cts + ' Counts'
      PLOT, dd_meg, pha(list_meg), TITLE = tit, XTIT = 'dd_meg (mm)', $
        YTIT = 'Pha', XRA = [0.0, MAX( dd_meg)], PSYM = 3

      xx = WHERE(dd_heg GE 0.0)
      cts = STRCOMPRESS(STRING(N_ELEMENTS(dd_heg(xx))),/REMOVE_ALL)

      tit = 'HEG Pha vs. dd: ' + cts + ' Counts'
      PLOT, dd_heg, pha(list_heg), TITLE = tit, XTIT = 'dd_heg (mm)', $
        YTIT = 'Pha', XRA = [0.0, MAX (dd_heg)], PSYM = 3

   ENDIF ELSE BEGIN       ;;; grating = "NONE" plots
      ymax = MAX(spectrum_first(1,*))
      xtit = 'E (keV)'
      ytit = 'Counts'
      tit = ' PHA spectrum '

      PLOT_OO, spectrum_first(0,*), spectrum_first(1,*), XRA=[0.1,10], $
        XTIT = xtit, YTIT = ytit, YRA = [1, ymax], XSTY = 8, YMARGIN = [4, 4], $
        XMARGIN = [8, 8], SUBTITLE = tit

      AXIS, XAXIS = 1, XRANGE = ckev / (10. ^ !X.CRANGE), XSTYLE = 1, $
        XTIT =  ' Wavelength (' + angstrom + ')'
      AXIS, YAXIS = 1, YRANGE = (10. ^ !Y.CRANGE) / ttime, YSTYLE = 1, $
        YTIT = 'Counts/s'

   ENDELSE
ENDELSE

tit = 'Input Spectrum'
maxflux = MAX(inputspec(1,*))

;;
PLOT_OO, inputspec(0,*), inputspec(1,*), XTIT = 'E (keV)', $
 YTIT = 'Flux (photons cm!E-2!n s!E-1!n keV!E-1!n!x)', $
 XRA = [0.03, 10], $
 YRA = [1.0e-8*maxflux, maxflux], XSTY = 8, YMARGIN = [4, 4], $
 SUBTITLE = tit

AXIS, XAXIS = 1, XRANGE = ckev / (10. ^ !X.CRANGE), XSTYLE = 1, $
  XTIT = 'Wavelength (' + angstrom + ') ' 
;;


;; Now print the simulation information on the plot
total_eff = FLOAT(counts) / FLOAT(numrays)

string1 = ' !11YE SIMULATION INFORMATION!x '
string2 = '  Grating: ' + grating
string3 = '  Detector: ' + detector
string4 = '  NumRays: ' + $
  STRCOMPRESS(STRING(numrays), /REMOVE_ALL) + ' phot'
string5 = '  Exposure Time: ' + $
  STRCOMPRESS(STRING(ttime), /REMOVE_ALL) + ' sec'
string6 = '  Detected Counts: ' + $
  STRCOMPRESS(STRING(counts), /REMOVE_ALL)
string7 = '  Total Efficiency: ' + $
  STRCOMPRESS(STRING(total_eff), /REMOVE_ALL)

PLOT, [0, 1], [0, 1], /NODATA, YSTYLE = 4, XSTYLE = 4, TITLE = string1
;XYOUTS, 0.0, 0.9, string1
XYOUTS, 0.0, 0.8, string2
XYOUTS, 0.0, 0.7, string3
XYOUTS, 0.0, 0.6, string4
XYOUTS, 0.0, 0.5, string5
XYOUTS, 0.0, 0.4, string6
XYOUTS, 0.0, 0.3, string7

DEVICE, /CLOSE
SET_PLOT,'x'

PRINT,' '
PRINT,'        The output plots have been written to ' + psfile 
PRINT,'  '

!P.MULTI = [0, 1, 1]
END