PRO main_hetg_m, mord, READ_DATA = read_data, OLD_FORMAT=old_format, $
	USERSUPPORT=usersupport, ACISI = acisi
;% Time-stamp: <97/05/20 10:46:02 dph> 
;% MIT Directory: ~dph/h1/HETG/Aeff_calc
;% CfA Directory: /proj/asc/dph/Aeff_calc
;% File: main_hetg_m.pro
;% Author: D. Huenemoerder
;% Original version: 970416
;%=====================================================================
; 10/2/97 Arg!  Dan Dewey here to modify this for:
;  - Use rdb tables that are in /spectra/d6/CIP/ for input
; 10/7/97 Add USERSUPPORT option to output PIMMS tables
; 10/7/97 Add ACISI option
;-----------------------------------------------------------------------
; main for doing eff area tables.

; usage:
; set data source: use_cips = 1  ; uses /nfs/spectra/d6/CIP/ data
; set order: mord=1 (valid values for HETG are 0 to +11)
; to read data tables, set do_read=1
; to restore previously saved data tables (by this program), set
; do_restore=1 and do_read=0
; Then: .run main_hetg
;       See  the data file definitions for the input tables. (just
;       after "IF do_read THEN BEGIN"

; This program writes several ascii tables, of kev, value:
; EA_HETG-AS.tbl - total HETG-ACIS-S effective area
; EA_HEG-AS.tbl  - HEG-ACIS-S effective area
; EA_MEG-AS.tbl  - MEG-ACIS-S effective area
; ACIS_OBF.tbl   - acis optical blocking filter transmission (in case
; it changes, you can fix the EA tables easily).

; after running this file, you can see the eff area via:
; plot,ehetg,ea_hetgs
; oplot,ehetg,ea_heg_mrg
; oplot,ehetg,ea_meg_mrg
; (and variations on that theme)

@AXAF_params

IF n_elements(mord) EQ 0 THEN mord = 1

IF KEYWORD_SET(READ_DATA) THEN begin 
  do_read = 1
  do_restore = 0
end else begin
  do_read = 0
  do_restore = 1
end

; compute limits of spectra in mm along dispersion - use for
; fundamental axis.
; The following assumes that heg and meg angles are symmetric about
; detector x-axis (axaf y) and uses the average angle as the angle.

IF KEYWORD_SET(ACISI) then begin
  acisi_geometry,  rgaps, lgaps
end else begin
  aciss_geometry,  rgaps, lgaps
end

; generate an array of positions on a single-pixel grid, from
; min(lgaps) - the left edge at negative r, to max(rgaps), the right
; edge at positive r.

pix = 0.024                     ; mm/pix for acis
rmin = min(lgaps)
rmax = max(rgaps)
npix = fix( (rmax-rmin) / pix)
rr = dindgen(npix)/(npix-1) * (rmax-rmin)+rmin

; Use the dispersion relations to assign energies to
; these locations
; For zero-order go to low energies by changing the rmin, rmax, rr
if mord EQ 0 then begin
  ; special for zero-order...  array is +/- 0.05 keV for MEG
  pix = 0.2                     ; mm/pix for acis
  rmin = -1.*!x_tg*(!hc/0.07)/!meg_period
  rmax = !x_tg*(!hc/0.07)/!meg_period
  npix = fix( (rmax-rmin) / pix)
  rr = dindgen(npix)/(npix-1) * (rmax-rmin)+rmin
  emax = 10.0                     ; truncate at 10 keV
  emeg = 1 * !hc * !x_tg / !meg_period / abs(rr)
  emeg = emeg < emax
  ; ditto for HEG:
  eheg = 1 * !hc * !x_tg / !heg_period / abs(rr)
  eheg = eheg < emax
end else begin
  emax = 10.0                     ; truncate at 10 keV
  emeg = mord * !hc * !x_tg / !meg_period / abs(rr)
  emeg = emeg < emax
  ; ditto for HEG:
  eheg = mord * !hc * !x_tg / !heg_period / abs(rr)
  eheg = eheg < emax
end

;; compute "visibility" array: =1 where there is a detector, =0 where
;;                             not, using gaps and rr

vis = fltarr(npix)+1.0

; skip if zero-order
if mord NE 0 then begin
  off_array = where( (rr LT rmin) OR (rr GT rmax) )
  gap_S0_S1 = where( (rr GT lgaps[5]) AND (rr LT lgaps[4]) )
  gap_S1_S2 = where( (rr GT lgaps[3]) AND (rr LT lgaps[2]) )
  gap_S2_S3 = where( (rr GT lgaps[1]) AND (rr LT lgaps[0]) )
  gap_S3_S4 = where( (rr GT rgaps[0]) AND (rr LT rgaps[1]) )
  gap_S4_S5 = where( (rr GT rgaps[2]) AND (rr LT rgaps[3]) )

  IF (size(off_array))[0] NE 0 THEN vis[off_array] = 0
  vis[gap_S0_S1] = 0
  vis[gap_S1_S2] = 0
  vis[gap_S2_S3] = 0
  vis[gap_S3_S4] = 0
  vis[gap_S4_S5] = 0
end

; read data files, and interpolate to meg, heg acis energy grids:
!Sav_DIR = './Sav/'
!Tbl_DIR = './Tbl/'
IF !version.os EQ 'MacOS' THEN !Sav_DIR = ':Sav:'
IF !version.os EQ 'MacOS' THEN !Tbl_DIR = ':Tbl:'

IF do_read THEN BEGIN
   if KEYWORD_SET(OLD_FORMAT) then BEGIN
   ;-----------------------------   

   ;; could reinsert Dave's read stuff here if desired...

   
   ;-----------------------------   
   ; end of old Data read in
   end else begin
     ; Here get data from dd's /nfs/spectra/d6/CIPs
     cip_dir = !CIP_Dir  ; defined in AXAF_params.pro
     ; Use rdb_read

   IF !version.os EQ 'MacOS' THEN BEGIN
      Data_dir = ':Data:'
   ENDIF ELSE Data_dir = cip_dir
   
   DQE_file = Data_dir + 'acis_qe.rdb' ; acis QE(E)
   DFT_file = Data_dir + 'acis_fil.rdb' ; acis filter transmission(E)
   HRMA_EA_file = Data_dir + 'hrma.rdb' ; HRMA effective area
   
   MEG_1_file = Data_dir + 'HETG_shell1_effic.rdb' ; dan's and john's membrane effic.
   MEG_3_file = Data_dir + 'HETG_shell3_effic.rdb' ;   by shell, includes vign.
   HEG_4_file = Data_dir + 'HETG_shell4_effic.rdb'
   HEG_6_file = Data_dir + 'HETG_shell6_effic.rdb'
   
   acis_qe = rdb_read(dqe_file)
   det_e = acis_qe.energy
   dqe_bi = acis_qe.bi
   dqe_fi = acis_qe.fi

   acis_fil = rdb_read(dft_file)
   filt_e = acis_fil.energy
   t_filt = acis_fil.s   
   ti_filt = acis_fil.i

   hrma = rdb_read(HRMA_EA_file)
   hrma_e = hrma.energy
   hrma_wav = !hc/hrma_e
   atot = hrma.all
   a1 = hrma.shell1
   a3 = hrma.shell3
   a4 = hrma.shell4
   a6 = hrma.shell6

; HETG tables have energy, then effic from -11 to +11 orders; 
; same for rdb tables
   maxord = 11
   
   message, 'Begin reading HETG effic files...', /info
foo = systime(1)
   meg1 = rdb_read(meg_1_file)
   m1_e = meg1.energy
   meg_eff1 = meg1(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', meg_1_file

   message, '... continue reading HETG effic files...', /info
foo = systime(1)
   meg3 = rdb_read(meg_3_file)
   m3_e = meg3.energy
   meg_eff3 = meg3(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', meg_3_file

   message, '... continue reading HETG effic files...', /info
foo = systime(1)
   heg4 = rdb_read(heg_4_file)
   h4_e = heg4.energy
   heg_eff4 = heg4(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', heg_4_file

   message, '... continue reading HETG effic files...', /info
foo = systime(1)
   heg6 = rdb_read(heg_6_file)
   h6_e = heg6.energy
   heg_eff6 = heg6(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', heg_6_file

   message, 'Done reading all tables!', /info

; save for quick restore:
   save, det_e, dqe_bi, dqe_fi, file = !Sav_DIR+'DQE.sav'
   save, filt_e, t_filt, ti_filt, file = !Sav_DIR+'DFT.sav'
   save, hrma_wav, hrma_e, atot, a1, a3, a4, a6, file = !Sav_DIR+'HRMA_EA.sav'
   save, m1_e, meg_eff1, m3_e, meg_eff3, h4_e, heg_eff4, h6_e, heg_eff6, $
    file=!Sav_DIR+'HETGS_eff'+'_m'+strcompress(string(mord),/rem)+'.sav'
   end   
ENDIF

IF do_restore THEN BEGIN
   restore, !Sav_DIR+'DQE.sav'
   restore, !Sav_DIR+'DFT.sav'
   restore, !Sav_DIR+'HRMA_EA.sav'
   restore, !Sav_DIR+'HETGS_eff'+'_m'+strcompress(string(mord),/rem)+'.sav'
ENDIF


; interpolate to common grid.

; "cm_" => Common Meg
;
message, 'Interpolating to MEG grid...', /info

cm_dqe_fi = interpol_sort(dqe_fi, det_e, emeg, /sil)
cm_dqe_bi = interpol_sort(dqe_bi, det_e, emeg, /sil)
if KEYWORD_SET(ACISI) then begin
  cm_t_filt = interpol_sort(ti_filt, filt_e, emeg, /sil)
end else begin
  cm_t_filt = interpol_sort(t_filt, filt_e, emeg, /sil)
end
cm_a1 = interpol_sort(a1, !hc/hrma_wav, emeg, /sil)
cm_a3 = interpol_sort(a3, !hc/hrma_wav, emeg, /sil)
cm_meg_eff1 = interpol_sort(meg_eff1, m1_e, emeg, /sil)
cm_meg_eff3 = interpol_sort(meg_eff3, m3_e, emeg, /sil)

; "ch_" => Common Heg
;
message, 'Interpolating to HEG grid...', /info

ch_dqe_fi = interpol_sort(dqe_fi, det_e, eheg, /sil)
ch_dqe_bi = interpol_sort(dqe_bi, det_e, eheg, /sil)
ch_t_filt = interpol_sort(t_filt, filt_e, eheg, /sil)
ch_a4 = interpol_sort(a4, !hc/hrma_wav, eheg, /sil)
ch_a6 = interpol_sort(a6, !hc/hrma_wav, eheg, /sil)
ch_heg_eff4 = interpol_sort(heg_eff4, h4_e, eheg, /sil)
ch_heg_eff6 = interpol_sort(heg_eff6, h6_e, eheg, /sil)


; merge

; determine where backside (BI) chips are in rr, emeg, and eheg space:
; BI's are in S1 and S3 on ACIS-S

; ACIS-S
if mord EQ 0 then begin
  ; all backside for zero-order
  ll_BI = indgen(n_elements(rr))
end else begin
  ; S1 is from lgap[3:4]
  ; S3 is from lgap[0] to rgap[0]
  ; calculate a list of pixel indices where array is BI:
  ;
  ll_BI = where(  ( (rr LE lgaps[3]) AND (rr GE lgaps[4]) ) $
                OR $
                ( (rr LE rgaps[0]) AND (rr GT lgaps[0]) ) $
             )
end

IF KEYWORD_SET(ACISI) THEN begin
  cm_dqe_hybrid = cm_dqe_fi
  ch_dqe_hybrid = ch_dqe_fi
end else begin
  ; ACIS-S: merge the detector QE arrays:
  cm_dqe_hybrid = cm_dqe_fi
  cm_dqe_hybrid[ll_BI] = cm_dqe_bi[ll_BI]
  ch_dqe_hybrid = ch_dqe_fi
  ch_dqe_hybrid[ll_BI] = ch_dqe_bi[ll_BI]
end

; apply vignetting terms, sum 1&3 into MEG Aeff, 4&6 into HEG Aeff:

EA_MEG = (cm_a1*cm_meg_eff1 + cm_a3*cm_meg_eff3)*!hrma_vign * !HETG_vign $
 * cm_dqe_hybrid * cm_t_filt * vis

EA_HEG = (ch_a4*ch_heg_eff4 + ch_a6*ch_heg_eff6)*!hrma_vign * !HETG_vign $
 * ch_dqe_hybrid * ch_t_filt * vis


; need to "fold" in half about rr=0.  Define a final grid, min([emeg,eheg]) to
; max([emeg,eheg]). 
; Interpolate left and right emeg, eheg onto that grid.
; Sum left and right Ae's.

ll_lt0 = where(rr LT 0)
ll_ge0 = reverse(where(rr GE 0)) ; reverse puts in ascending egrid order


; make an r-grid from 0 to max(abs(rr)) at 1/2 pixels, then convert to
; e-grid at MEG scale.
; This is the final output grid - do log spacing for 0 order
if mord EQ 0 then begin
  egrid = 10.0^( FLOAT(indgen(4001))/2000.0 -1.0 )
end else begin
  rgrid = findgen(npix)/(npix-1) * max(abs(rr))
  egrid = reverse ( (mord * !hc * !x_tg / !meg_period / rgrid) <  emax )  
end

EA_meg_plus = EA_MEG(ll_ge0)
EA_meg_minus = EA_MEG(ll_lt0)
EA_heg_plus = EA_HEG(ll_ge0)
EA_heg_minus = EA_HEG(ll_lt0)

em_plus = emeg(ll_ge0)
em_minus = emeg(ll_lt0)
eh_plus = eheg(ll_ge0)
eh_minus = eheg(ll_lt0)

; interpolate and add
message,'Interpolate and add',/info
m_ea_p = interpol_sort(ea_meg_plus,em_plus,egrid, /sil)
m_ea_m = interpol_sort(ea_meg_minus,em_minus,egrid, /sil)
h_ea_p = interpol_sort(ea_heg_plus,eh_plus,egrid, /sil)
h_ea_m = interpol_sort(ea_heg_minus,eh_minus,egrid, /sil)

; "fix" heg areas - interpol extrapolates, need to zero off-chip
; areas.
ll = where(egrid LT min(eh_plus),nll)
if nll GE 1 then begin
  h_ea_p(ll) = 0
end
ll = where(egrid LT min(eh_minus),nll)
if nll GE 1 then begin
  h_ea_m(ll) = 0
end
; for MEG too
ll = where(egrid LT min(em_plus),nll)
if nll GE 1 then begin
  m_ea_p(ll) = 0
end
ll = where(egrid LT min(em_minus),nll)
if nll GE 1 then begin
  m_ea_m(ll) = 0
end

; Merge and set negative values to zero
if mord EQ 0 then begin
  EA_MEG_mrg = ((m_ea_p) > 0.0)
  EA_HEG_mrg = ((h_ea_p) > 0.0)
  EA_HETGS = ea_meg_mrg + ea_heg_mrg
end else begin
  EA_MEG_mrg = ((m_ea_p + m_ea_m) > 0.0)
  EA_HEG_mrg = ((h_ea_p + h_ea_m) > 0.0)
  EA_HETGS = ea_meg_mrg + ea_heg_mrg
end

; DONE! (except for tables and plots)
if KEYWORD_SET(USERSUPPORT) then begin
  ; Use fewer points for PIMMS files...
  npts = 512 ; for user support pimms use
  ; but remove all the clipped values
  to_clip = where(egrid GE emax, nwhere)
  if nwhere GE 2 then begin
    select = INDGEN(min(to_clip)+1)
  end else begin
    select = INDGEN(n_elements(egrid))
  end
  ord_str = STRCOMPRESS(STRING(mord,FORMAT='(I2)'),/REMOVE_ALL)
  if KEYWORD_SET(ACISI) then begin
    det_str = 'acis-i'
  end else begin
    det_str = 'acis-s'
  end
  output_tables,egrid(select),ea_hetgs(select),npts,!USTbl_DIR+ $
	'axaf_hetg-'+det_str+'_heg'+ord_str+'meg'+ord_str+'.area',$
	/regular
  output_tables,egrid(select),ea_heg_mrg(select),npts,!USTbl_DIR+ $
	'axaf_hetg-'+det_str+'_heg'+ord_str+'.area',$
	/regular
  output_tables,egrid(select),ea_meg_mrg(select),npts,!USTbl_DIR+ $
	'axaf_hetg-'+det_str+'_meg'+ord_str+'.area',$
	/regular
  ; since filter may change...
  ftrans = interpol_sort(t_filt,filt_e,egrid, /sil)
  output_tables,egrid(select),ftrans(select),npts,!USTbl_DIR+ $
	'axaf_'+det_str+'_obf.tbl',/regular
end else begin
  npts = n_elements(egrid) ; do 'em all! (for detailed use)
  ord_str = STRCOMPRESS(STRING(mord,FORMAT='(I2)'),/REMOVE_ALL)
  if KEYWORD_SET(ACISI) then begin
    det_str = 'AI'
  end else begin
    det_str = 'AS'
  end
  output_tables,egrid,ea_hetgs,npts,!Tbl_DIR+'EA_HETG'+ord_str+'-'+det_str+'.rdb',$
	/regular,/rdb
  output_tables,egrid,ea_heg_mrg,npts,!Tbl_DIR+'EA_HEG'+ord_str+'-'+det_str+'.rdb',$
	/regular,/rdb
  output_tables,egrid,ea_meg_mrg,npts,!Tbl_DIR+'EA_MEG'+ord_str+'-'+det_str+'.rdb',$
	/regular,/rdb
  ; since filter may change...
  ftrans = interpol_sort(t_filt,filt_e,egrid, /sil)
  output_tables,egrid,ftrans,npts,!Tbl_DIR+'ACIS_'+det_str+'_OBF.rdb',/regular,/rdb
end

RETURN
END