PRO eae_flux_anal, eae
; Go through an eae array and calculate values for:
;  - flux method
;  - flux and flux_err
;  - focal_method
;  - focal_qe
;  - focal_line, focal_line_plus, focal_line_minus
;  - effective_method
;  - effective, 'err_plus, 'err_minus, 'syst_plus, 'syst_minus
;
; 11/19/97 dd Use the detector QE curves to have flux and
; focal_line in photons instead of counts.
; 6/29/98 dd Use new fpc qe fits
; 4/6/99 dd Use CALDBaxafcal/xrcf/cip files...

; FPC QEs
; Tags: Energy (keV) and QE
fpc_x2_qe = rdb_read(!DDXRCFCAL+ $
		'/cip/fpc_x2.fitQe_N19980902.rdb')
fpc_5_qe = rdb_read(!DDXRCFCAL+ $
		'/cip/fpc_5.fitQe_N19980902.rdb')
fpc_hn_qe = rdb_read(!DDXRCFCAL+ $
		'/cip/fpc_hn.fitQe_N19980902.rdb')

; SSD QEs
; Tags: 990208 energy is in eV and QE
ssd_5_qe = rdb_read(!DDXRCFCAL+ $
	'/cip/ssd_5_qe_N19990208.rdb')
ssd_5_qe.energy = ssd_5_qe.energy/1000.0
ssd_x_qe = rdb_read(!DDXRCFCAL+ $
	'/cip/ssd_x_qe_N19990208.rdb')
ssd_x_qe.energy = ssd_x_qe.energy/1000.0

; Arrays for the fluxes
hn_flux = fltarr(n_elements(eae))
f5_flux = fltarr(n_elements(eae))
s5_flux = fltarr(n_elements(eae))
hn_flux_err = fltarr(n_elements(eae))
f5_flux_err = fltarr(n_elements(eae))
s5_flux_err = fltarr(n_elements(eae))
 
for ir=0,n_elements(eae)-1 do begin

  ; Calculate the fluxes first (counts/s cm^2)

  ; FPC_HN
  area = eae_bnd_hrma_area('FPC_HN',eae(ir).fpc_hn_aperture)
  if area LT 0.0 OR eae(ir).fpc_hn_time LT 0.1 then begin
    hn_flux(ir) = -1.0  ; not valid      
  end else begin
    hn_flux(ir) = eae(ir).fpc_hn_counts/( $
	eae(ir).fpc_hn_time * area)
    hn_flux_err(ir) = eae(ir).fpc_hn_err/( $
	eae(ir).fpc_hn_time * area)
  end
  ; FPC_5
  area = eae_bnd_hrma_area('FPC_5',eae(ir).fpc_5_aperture)
  if area LT 0.0 OR eae(ir).fpc_5_time LT 0.1 then begin
    f5_flux(ir) = -1.0  ; not valid      
  end else begin
    f5_flux(ir) = eae(ir).fpc_5_counts/( $
	eae(ir).fpc_5_time * area)
    f5_flux_err(ir) = eae(ir).fpc_5_err/( $
	eae(ir).fpc_5_time * area)
  end

  ; SSD_5
  area = eae_bnd_hrma_area('SSD_5',eae(ir).ssd_5_aperture)
  if area LT 0.0 OR eae(ir).ssd_5_time LT 0.1 then begin
    s5_flux(ir) = -1.0  ; not valid      
  end else begin
    s5_flux(ir) = eae(ir).ssd_5_counts/( $
	eae(ir).ssd_5_time * area)
    s5_flux_err(ir) = eae(ir).ssd_5_err/( $
	eae(ir).ssd_5_time * area)
  end

;  - flux method, flux, and flux_err
; Choose the flux method as appropriate and
; convert flux to photons/s cm^2  by dividing by the QE.
; Default is not flux set:
  eae(ir).flux_method = 'No flux available (eae_flux_anal.pro)'
  eae(ir).flux = -1.0
  eae(ir).flux_err = 0.0
  ; if SSD or ACIS or 2C prefer to use ssd_5 flux
  if STRPOS(eae(ir).fp_detector,'SSD') GE 0 OR $
	STRPOS(eae(ir).fp_detector,'ACIS') GE 0 OR $
	STRPOS(eae(ir).fp_detector,'2C') GE 0 then begin
    eae(ir).flux_method = 'SSD_5 only (eae_flux_anal.pro)'
    ; Get SSD_5 QE at this energy
    above = where(ssd_5_qe.energy GE eae(ir).energy)
    eae(ir).flux = s5_flux(ir)/ssd_5_qe(above(0)).qe
    eae(ir).flux_err = s5_flux_err(ir)/ssd_5_qe(above(0)).qe
  end else begin
    ; FPC, HSI, HRC in focal plane
    ; Use HN if available
    if hn_flux(ir) GT 0.0 then begin
      eae(ir).flux_method = 'FPC_HN only (eae_flux_anal.pro)'
      ; Get fpc_hn QE at this energy
      above = where(fpc_hn_qe.energy GE eae(ir).energy)
      eae(ir).flux = hn_flux(ir)/fpc_hn_qe(above(0)).qe
      eae(ir).flux_err = hn_flux_err(ir)/fpc_hn_qe(above(0)).qe
    end else begin
      ; if no HN then use fpc_5
      if f5_flux(ir) GT 0.0 then begin
        eae(ir).flux_method = 'FPC_5 only (eae_flux_anal.pro)'
        ; Get fpc_5 QE at this energy
        above = where(fpc_5_qe.energy GE eae(ir).energy)
        eae(ir).flux = f5_flux(ir)/fpc_5_qe(above(0)).qe
        eae(ir).flux_err = f5_flux_err(ir)/fpc_5_qe(above(0)).qe
      end else begin
        if s5_flux(ir) GT 0.0 then begin
          ; Finally use SSD_5
          eae(ir).flux_method = 'SSD_5 only (eae_flux_anal.pro)'
          ; Get ssd_5 QE at this energy
          above = where(ssd_5_qe.energy GE eae(ir).energy)
          eae(ir).flux = s5_flux(ir)/ssd_5_qe(above(0)).qe
          eae(ir).flux_err = s5_flux_err(ir)/ssd_5_qe(above(0)).qe
        end
      end
    end
  end

;  - focal_method
  eae(ir).focal_method = 'photons/s (eae_flux_anal.pro)'
;  - focal_qe
  ; Find focal QE based on focal detector
  CASE eae(ir).fp_detector OF
    'FPC2': begin
      above = where(fpc_x2_qe.energy GE eae(ir).energy)
      fp_qe = fpc_x2_qe(above(0)).qe
    end
    'SSD': begin
      above = where(ssd_x_qe.energy GE eae(ir).energy)
      fp_qe = ssd_x_qe(above(0)).qe
    end
    'HSI': begin
      fp_qe = 1.0
    end
    '2C1': begin
      fp_qe = 1.0
    end
    '2C0': begin
      fp_qe = 1.0
    end
    'HRC,I': begin
      fp_qe = 1.0
    end
    'HRC,S': begin
      fp_qe = 1.0
    end
    'ACIS,S3': begin
      fp_qe = 1.0
    end
    ELSE: begin
      print, 'eae_flux_anal.pro : unrecognized fp_detector:'+ $
	eae(ir).fp_detector
      fp_qe = 1.0
    end
  ENDCASE
  eae(ir).focal_qe = fp_qe

;  - focal_line, focal_line_plus, focal_line_minus
  qe_fact = 1.0/eae(ir).focal_qe
  eae(ir).focal_line = qe_fact*eae(ir).focal_counts / $
			eae(ir).focal_time
  eae(ir).focal_line_plus = qe_fact*eae(ir).focal_err_plus / $
			eae(ir).focal_time
  eae(ir).focal_line_minus = qe_fact*eae(ir).focal_err_minus / $
			eae(ir).focal_time
;  - effective_method
  eae(ir).effective_method = 'BNDREF (eae_flux_anal.pro)'
;  - effective, 'err_plus, 'err_minus, 'syst_plus, 'syst_minus
  eae(ir).effective = eae(ir).focal_line / eae(ir).flux
  eae(ir).effective_err_plus = eae(ir).focal_line_plus / eae(ir).flux
  eae(ir).effective_err_minus = eae(ir).focal_line_minus / eae(ir).flux
; For this analysis diffracted orders can be underestimated
  if eae(ir).grating EQ 'NONE' then begin
    eae(ir).effective_syst_plus = eae(ir).effective * 0.05
    eae(ir).effective_syst_minus = eae(ir).effective * 0.05
  end else begin
    if eae(ir).order EQ 0.0 then begin
      eae(ir).effective_syst_plus = eae(ir).effective * 0.10
      eae(ir).effective_syst_minus = eae(ir).effective * 0.10
    end else begin
      eae(ir).effective_syst_plus = eae(ir).effective * 0.40
      eae(ir).effective_syst_minus = eae(ir).effective * 0.10
    end
  end
end

RETURN
END