PRO grating_spec, marxdir, Spectrum, GRATING = grating, SCALE = scale, ORDER = order
;;+
;
; PURPOSE: 
;
;  Grating_spec.pro will extract a spectrum from MARX
;  grating simulation output. The spectrum is a histogram,
;  generated by tg_extract.pro, of cts/sec/channel versus 
;  energy (in keV; by default) or wavelength (Angstroms).
;
; INPUT: 
;         
;  marxdir (string) :: Directory containing the marx output
;
;
; OUTPUT:
;
;  spectrum (vector) :: Output vector, of dimension either
;                       2 x N or 4 x N, containing the
;                       binned histogram and the energy
;                       or wavelength of each bin (N = the
;                       number of bins)
;
;   NOTE that the units of the spectrum is cts/sec/bin,
;   *NOT*, e.g., cts/kev/sec or cts/ang/sec
;
;  IF the HETG was used, spectrum is a 4 x N vector - TWO
;  spectra are returned. The first two columns of spectrum
;  contain the energy (or wavelength) and histogram of the
;  MEG data. The last two columns contain the energy (or
;  wavelength) and histogram of the HEG data. It the LETG
;  is used, spectrum is a 2 x N vector containing the energy
;  (or wavelength) and histogram of the LEG data.
;
; KEYWORDS:
;
;   GRATING - Set this equal to the grating; 'HETG'/'hetg'/'H'/'h',
;             or likewise for LEG. If this keyword is not
;             set, the program will determine the grating
;             from the marx.par file.
;
;   SCALE   - This determines the "X-axis" of the spectrum, i.e.
;             whether the scale used is a wavelength scale, or
;             an energy scale. If this keyword is NOT set,
;             the default is to use an ENERGY (keV) scale.
;             Useage: SCALE = 'wave' or SCALE = 'lambda'
;   
;   ORDER   - This determines for which order to extract a spectrum.
;             The program will examine the order.dat file and extract
;             only those events for a given order. If this
;             keyword is NOT set, the default is to use FIRST order
;             ( +/- 1) only. 
;
;  AUTHOR: Chris Baluta, 05/01/97
;
;;-

IF ( N_PARAMS() NE 2) THEN BEGIN
  PRINT,'      '
  PRINT,'   grating_spec  USAGE: '
  PRINT,'      '
  PRINT,'grating_spec, marxdir, Spectrum, GRATING = grating, ' + $
    'SCALE = scale, ORDER = order'
  PRINT,'      '
  PRINT,'   INPUT: '
  PRINT,'         '
  PRINT,'         marxdir: (String) The directory containing the marx output'
  PRINT,'        '
  PRINT,'   OUTPUT: '
  PRINT,'          '
  PRINT,'         spectrum: (vector) A 2xN *or* a 4xN dimensional vector'
  PRINT,'                   containing the output spectrum/spectra'
  PRINT,'                   (N = the number of bins in the spectrum/spectra) '
  PRINT,'                   The first column contains the energy (keV) '
  PRINT,'                   or wavelength (angstroms) of the bin while the '
  PRINT,'                   second column contains the histograms (cts/sec/bin)'
  PRINT,'                   IF the grating is HETG, columns 1,2 are the MEG '
  PRINT,'                   spectrum, columns 3,4 are the HEG spectrum'
  PRINT,'        '
  PRINT,'   KEYWORDS: '
  PRINT,'    '
  PRINT,'      GRATING : The grating to use ("hetg"/"HETG"/"h"/"H", for ex.).'
  PRINT,'                If this keyword is NOT set, the program will determine'
  PRINT,'                the grating from the marx.par file '
  PRINT,'      '
  PRINT,'      SCALE   : The scale of the "X-axis" of the spectrum (energy or'
  PRINT,'                wavelength). Thus, eg, SCALE = "wave" or SCALE =' + $
        ' "lambda"'
  PRINT,'                If this keyword is not set, the default is to use '+ $
         'ENERGY (keV) '
  PRINT,'  '
  PRINT,'      ORDER   : The order for which to extract the spectrum. If this'
  PRINT,'                keyword is not set, the default is to use the FIRST'+$
        ' order'
  PRINT,'      '
  RETURN
ENDIF

TRUE = 1
FALSE = 0

IF ( NOT KEYWORD_SET(ORDER) ) THEN order = 1

IF ( KEYWORD_SET (SCALE)) THEN BEGIN
 wave_scale = TRUE 
ENDIF ELSE BEGIN
  wave_scale = FALSE 
ENDELSE

IF ( KEYWORD_SET(GRATING)) THEN BEGIN
  CASE STRUPCASE(STRMID(grating, 0, 1)) OF
       'H': grat = 'HETG'
       'L': grat = 'LETG'
       'N': grat = 'NONE'
     ELSE : BEGIN
             PRINT,'  '
             PRINT,'   Illegal value for KEYWORD GRATING'
             PRINT,'   GRATING='+grating
             PRINT,'  '
             PRINT,' GRATING must be "h","H","HETG","hetg" or'
             PRINT,' "l","L","letg", or "LETG" '
             PRINT,'   '
             RETURN
            END
  ENDCASE
ENDIF ELSE BEGIN
  grat = read_parfile ('GratingType', MARXDIR = marxdir)
ENDELSE

IF (STRMID(grat, 0, 1)  EQ '"') THEN grat = STRMID(grat, 1, 4)

;;; Now get the grating parameters

IF ( KEYWORD_SET(GRATING) ) THEN BEGIN     ;;; Assume that there is no .par file
 hegPeriod = 2000.                    ;;; NEED DEFAULT STRUCTURE VALUES HERE
 megPeriod = 4000. 
 legPeriod = 9912.
 nogPeriod = 0.0
 rowland = 8634.    ;; mm
ENDIF ELSE BEGIN
 hegPeriod = ( read_parfile ('hegPeriod', MARXDIR = marxdir) ) * 1.0e4
 megPeriod = ( read_parfile ('megPeriod', MARXDIR = marxdir) ) * 1.0e4
 legPeriod = ( read_parfile ('legPeriod', MARXDIR = marxdir) ) * 1.0e4
 rowland = read_parfile ('RowlandDiameter', MARXDIR = marxdir)
ENDELSE


;;; Some other constants
detector = read_parfile ('DetectorType', MARXDIR = marxdir)
IF (STRMID(detector, 0, 1) EQ 'A') THEN BEGIN
  pixel = read_parfile('ACISPixelSize', MARXDIR = marxdir)  ;;;.024 
ENDIF ELSE BEGIN
  pixel = read_parfile('HRCPixelSize', MARXDIR = marxdir) 
ENDELSE

kev_angstroms = 12.39854            ;;; NEED DEFAULT STRUC. VALUES HERE

;;;

;;; Now read in the data
PRINT,'  '
PRINT,'    Reading in data from ' + marxdir + ' now.'
PRINT,'  '

;;;;;
ypos = read_marx_file (marxdir + 'ypos.dat')
IF (N_ELEMENTS(ypos) EQ -1) THEN RETURN

zpos = read_marx_file (marxdir + 'zpos.dat') 
IF (N_ELEMENTS(zpos) EQ -1) THEN RETURN

IF (grat NE 'NONE')  THEN BEGIN
  orders = read_marx_file (marxdir + 'order.dat')
  IF (N_ELEMENTS(orders) EQ -1) THEN RETURN
ENDIF 

time = read_marx_file (marxdir + 'time.dat')
IF (N_ELEMENTS(time) EQ -1) THEN RETURN
;;;;;

PRINT,'    Successfully read in ypos, zpos, order, and time data.'
PRINT,'   '

;; Now filter on the chosen order
IF (grat EQ 'NONE') THEN BEGIN
  ypos_ord = ypos
  zpos_ord = zpos
ENDIF ELSE BEGIN
  ord = WHERE (ABS(orders) EQ order)
  ypos_ord = ypos (ord)
  zpos_ord = zpos (ord)
  time_ord = time (ord)
ENDELSE
;;

y0 = 0.0
z0 = 0.0        ;;; Assume 0-order is at (0.0, 0.0)
maxdist = 0.2    ;;; mm 

delta_t = MAX (time) - MIN (time)

;;; Now extract the spectrum
CASE grat OF 
 'HETG': BEGIN
           tg_extract, 'h', maxdist, ypos_ord, zpos_ord, y0, z0, dd_heg, dxd_heg 
           tg_extract, 'm', maxdist, ypos_ord, zpos_ord, y0, z0, dd_meg, dxd_meg

           heg_min = 0.0
           meg_min = 0.0

           heg_spec = ( HISTOGRAM (dd_heg, BINS=pixel, MIN=heg_min) ) ;;/ delta_t
           meg_spec = ( HISTOGRAM (dd_meg, BINS=pixel, MIN=meg_min) ) ;;/ delta_t

           num_heg_bins = N_ELEMENTS (heg_spec) 
           num_meg_bins = N_ELEMENTS (meg_spec) 

           IF (num_heg_bins GT num_meg_bins) THEN BEGIN
             max_bins = num_heg_bins - 1 
             fill = num_meg_bins - 1   
           ENDIF ELSE BEGIN
             max_bins = num_meg_bins - 1
             fill = num_heg_bins - 1 
           ENDELSE
 
           spectrum = FLTARR( 4, max_bins)

           const_heg = hegPeriod / rowland
           const_meg = megPeriod / rowland
          
           x_heg = FINDGEN (num_heg_bins) * pixel + heg_min
           x_meg = FINDGEN (num_meg_bins) * pixel + meg_min 

           wave_heg = const_heg * x_heg(1:*) / order
           wave_meg = const_meg * x_meg(1:*) / order
 
           IF (NOT wave_scale) THEN BEGIN
             ene_heg = REVERSE (kev_angstroms / wave_heg)
             ene_meg = REVERSE (kev_angstroms / wave_meg)
             heg_spec = REVERSE (heg_spec(1:*))
             meg_spec = REVERSE (meg_spec(1:*))
            
             IF (max_bins EQ (num_heg_bins-1) ) THEN BEGIN             
               spectrum (0, 0:fill-1) = ene_meg
               spectrum (0, fill:*) = 1e-6 ;;; Arbitrary small number 
               spectrum (2, *) = ene_heg  
             ENDIF ELSE BEGIN
               spectrum (0, *) = ene_meg
               spectrum (2, 0:fill-1) = ene_heg
               spectrum (2, fill:*) = 1e-6 ;;; Arbitrary small number 
             ENDELSE
           ENDIF

           IF (max_bins EQ (num_heg_bins-1) ) THEN BEGIN
             spectrum (1, 0:fill-1) = meg_spec
             spectrum (1, fill:*) = 1e-6 
             spectrum (3, *) = heg_spec
           ENDIF ELSE BEGIN
             spectrum (1, *) = meg_spec
             spectrum (3, 0:fill-1) = heg_spec
             spectrum (3, fill:*) =  1e-6
           ENDELSE  

         END                ;;; CASE grat = 'HETG'
 'LETG': BEGIN
           tg_extract, 'l', maxdist, ypos_ord, zpos_ord, y0, z0, dd_leg, dxd_leg


           leg_min = 0.0
           leg_spec = ( HISTOGRAM (dd_leg, BINS=pixel, MIN=leg_min) ) ;;/ delta_t

           num_leg_bins = N_ELEMENTS (leg_spec)
           spectrum = FLTARR (2, num_leg_bins-1)

           const_leg = legPeriod / rowland
           x_leg = FINDGEN (num_leg_bins) * pixel + leg_min
           wave_leg = const_leg * x_leg(1:*) / order

           IF (NOT wave_scale) THEN BEGIN
             ene_leg = REVERSE (kev_angstroms / wave_leg)
             leg_spec = REVERSE (leg_spec(1:*))
             spectrum (0, *) = ene_leg
           ENDIF ELSE BEGIN
             spectrum = FLTARR (2, num_leg_bins)
             spectrum (0, *) = wave_leg
           ENDELSE 
             spectrum (1, *) = leg_spec
         END
 'NONE': BEGIN
             pha = (read_marx_file (marxdir + 'pha.dat')) * 12.0 / 4096.
             spec = HISTOGRAM(pha, BINS=.01)
             numbins = N_ELEMENTS(spec)

             ene_pha = FINDGEN(numbins) * .01
             wave_pha = REVERSE(kev_angstroms / ene_pha)

             spectrum = FLTARR(2, numbins)
             IF (NOT wave_scale) THEN BEGIN
               spectrum(0, *) = ene_pha
               spectrum(1, *) = spec
             ENDIF ELSE BEGIN
               spectrum(0, *) = wave_pha
               spectrum(1, *) = REVERSE(spec)
             ENDELSE 
         END
ENDCASE


END