PRO main_letg_m, mord, READ_DATA = read_data, OLD_FORMAT=old_format, $
	USERSUPPORT=usersupport, ACISI = acisi
;%=====================================================================
; 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
 ; 10/7/97 dd Copied from main_hetg_m
;-----------------------------------------------------------------------
; 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 LETG 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_letg
;       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_LETG-AS.tbl - total LETG-ACIS-S effective area
; ACIS_OBF.tbl   - acis optical blocking filter transmission (in case
; it changes, you can fix the EA tables easily).


@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
; As dph noted, the acis geom's have the angles included,
; change them for LETG angle
ave_angle = 0.5*(ABS(!meg_angle) + ABS(!heg_angle))
rgaps = rgaps*cos(ave_angle)/cos(!leg_angle)
lgaps = lgaps*cos(ave_angle)/cos(!leg_angle)

; 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...
  pix = 0.2                     ; mm/pix for acis
  rmin = -1.*!x_tg*(!hc/0.07)/!leg_period
  rmax = !x_tg*(!hc/0.07)/!leg_period
  npix = fix( (rmax-rmin) / pix)
  rr = dindgen(npix)/(npix-1) * (rmax-rmin)+rmin
  emax = 10.0                     ; truncate at 10 keV
  eleg = 1 * !hc * !x_tg / !leg_period / abs(rr)
  eleg = eleg < emax
end else begin
  emax = 10.0                     ; truncate at 10 keV
  eleg = mord * !hc * !x_tg / !leg_period / abs(rr)
  eleg = eleg < 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 leg, 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
   
   LEG_1_file = Data_dir + 'LETG_shell1_effic.rdb' ; dan's and john's membrane effic.
   LEG_3_file = Data_dir + 'LETG_shell3_effic.rdb' ;   by shell, includes vign.
   LEG_4_file = Data_dir + 'LETG_shell4_effic.rdb'
   LEG_6_file = Data_dir + 'LETG_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

; LETG tables have energy, then effic from -25 to +25 orders; 
; same for rdb tables
   maxord = 25
   
   message, 'Begin reading LETG effic files...', /info
foo = systime(1)
   leg1 = rdb_read(leg_1_file)
   l1_e = leg1.energy
   leg_eff1 = leg1(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', leg_1_file

   message, '... continue reading LETG effic files...', /info
foo = systime(1)
   leg3 = rdb_read(leg_3_file)
   l3_e = leg3.energy
   leg_eff3 = leg3(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', leg_3_file

   message, '... continue reading LETG effic files...', /info
foo = systime(1)
   leg4 = rdb_read(leg_4_file)
   l4_e = leg4.energy
   leg_eff4 = leg4(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', leg_4_file

   message, '... continue reading LETG effic files...', /info
foo = systime(1)
   leg6 = rdb_read(leg_6_file)
   l6_e = leg6.energy
   leg_eff6 = leg6(*).(maxord+1+mord)
print, systime(1)-foo, ' Seconds to read ', leg_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, l1_e, leg_eff1, l3_e, leg_eff3, l4_e, leg_eff4, l6_e, leg_eff6, $
    file=!Sav_DIR+'LETGS_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+'LETGS_eff'+'_m'+strcompress(string(mord),/rem)+'.sav'
ENDIF


; interpolate to common grid.

; "lm_" => Common Leg
;
message, 'Interpolating to LEG grid...', /info

lm_dqe_fi = interpol_sort(dqe_fi, det_e, eleg, /sil)
lm_dqe_bi = interpol_sort(dqe_bi, det_e, eleg, /sil)
if KEYWORD_SET(ACISI) then begin
  lm_t_filt = interpol_sort(ti_filt, filt_e, eleg, /sil)
end else begin
  lm_t_filt = interpol_sort(t_filt, filt_e, eleg, /sil)
end
lm_a1 = interpol_sort(a1, !hc/hrma_wav, eleg, /sil)
lm_a3 = interpol_sort(a3, !hc/hrma_wav, eleg, /sil)
lm_a4 = interpol_sort(a4, !hc/hrma_wav, eleg, /sil)
lm_a6 = interpol_sort(a6, !hc/hrma_wav, eleg, /sil)
lm_leg_eff1 = interpol_sort(leg_eff1, l1_e, eleg, /sil)
lm_leg_eff3 = interpol_sort(leg_eff3, l3_e, eleg, /sil)
lm_leg_eff4 = interpol_sort(leg_eff4, l4_e, eleg, /sil)
lm_leg_eff6 = interpol_sort(leg_eff6, l6_e, eleg, /sil)

; determine where backside (BI) chips are in rr, eleg 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
  lm_dqe_hybrid = lm_dqe_fi
end else begin
  ; ACIS-S: merge the detector QE arrays:
  lm_dqe_hybrid = lm_dqe_fi
  lm_dqe_hybrid[ll_BI] = lm_dqe_bi[ll_BI]
end

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

EA_LEG = (lm_a1*lm_leg_eff1 + lm_a3*lm_leg_eff3 + $
	lm_a4*lm_leg_eff4 + lm_a6*lm_leg_eff6 )*!hrma_vign * !LETG_vign $
	* lm_dqe_hybrid * lm_t_filt * vis

; need to "fold" in half about rr=0.  Define a final grid, min([eleg]) to
; max([eleg]). 
; Interpolate left and right eleg 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 LEG 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 / !leg_period / rgrid) <  emax )  
end

EA_leg_plus = EA_LEG(ll_ge0)
EA_leg_minus = EA_LEG(ll_lt0)

el_plus = eleg(ll_ge0)
el_minus = eleg(ll_lt0)

; interpolate and add
message,'Interpolate and add',/info
l_ea_p = interpol_sort(ea_leg_plus,el_plus,egrid, /sil)
l_ea_m = interpol_sort(ea_leg_minus,el_minus,egrid, /sil)

; "fix" leg areas - interpol extrapolates, need to zero off-chip
; areas.
ll = where(egrid LT min(el_plus),nll)
if nll GE 1 then begin
  l_ea_p(ll) = 0
end
ll = where(egrid LT min(el_minus),nll)
if nll GE 1 then begin
  l_ea_m(ll) = 0
end

; Merge and set negative values to zero
if mord EQ 0 then begin
  EA_LEG_mrg = ((l_ea_p) > 0.0)
  EA_LETGS = ea_leg_mrg 
end else begin
  EA_LEG_mrg = ((l_ea_p + l_ea_m) > 0.0)
  EA_LETGS = ea_leg_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_letgs(select),npts,!USTbl_DIR+ $
	'axaf_letg-'+det_str+'_order'+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_letgs,npts,!Tbl_DIR+'EA_LETG'+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