;+
;NAME: INTERPOL_EIS_EA
;
;PURPOSE:
;   Interpolate fit EIS effective areas to any date and wavelength
;
;INPUTS:
;   date - timestamp of the observation of interest (ISO format preferred)
;   wavelength - Wavelength of interest (in units of Angstroms). The EA curve
;       will be interpolated to the requested date and then the EA at
;       the given wavelength will be extracted (or interpolated further).
;
;       NB: If /short or /long are set, wavelength will instead become an OUTPUT
;       array containing the wavelength values for the selected waveband.
;       In this case, any input wave will be either cause an error or be ignored
;
;OPTIONAL INPUTS:
;   /short - If set, will output the EA in only the short wavelength band
;   /long - If set, will output the EA in only the long wavelength band
;   /radcal - If set, will output radiometric calibration factors for converting
;            DN to intensity
;   ea_file - filename and path containing compatible arrays of reference
;             effective areas and their timestamps. MUST include a 1D array
;             with string dates (UTC) named ALL_DATE_OBS and two 2D arrays
;             with the short and long waveband EA named, respectively,
;             ALL_FIT_EASW and ALL_FIT_EALW.
;
;OUTPUTS:
;   ea_out - 1D array or single value with the effective area in either the
;            selected waveband (default, short band) or singular wavelength.
;            If /radcal is set, will output calibration factors instead.
;
;OPTIONAL OUTPUTS:
;   None (at the moment)
;
;HISTORY:
;   Name----------Date----------Description
;   M Weberg    2022-May-02     Initial coding
;   M Weberg    2022-May-05     Checked math of /radcal method
;   H Warren    2022-May-16     Improved file search (added support for IDL<8.7)
;   M Weberg    2023-Oct-19     Enabled computing multiple wavelengths at once
;   M Weberg    2023-Dec-15     Added warning for dates outside fit dates
;
;TO-DO / LIMITATIONS:
;   - NONE (for now)
;-

FUNCTION INTERPOL_EIS_EA, date, wavelength, short=short, long=long, $
                          radcal=radcal, ea_file=ea_file, quiet=quiet

;@input#########################################################################
; Validate input values
;###############################################################################
IF N_ELEMENTS(date) NE 1 THEN MESSAGE, 'ERROR: please input a single date'

in_tai = ANYTIM2TAI(date)

IF in_tai LT ANYTIM2TAI('2006-10-20T10:20:00.000') THEN BEGIN
    ;Note: EIS was launched on '2006-09-22T21:36
    PRINT, ' WARNING: Selected date is before the start of normal '$
          +'EIS science operations. Output values may be inaccurate.'
ENDIF

IF NOT KEYWORD_SET(short) AND NOT KEYWORD_SET(long) THEN BEGIN
    n_input_wave = N_ELEMENTS(wavelength)
    loc_short = where(wavelength GE 165 and wavelength LE 213, n_short)
    loc_long = where(wavelength GE 245 and wavelength LE 292, n_long)
    IF (n_short+n_long LT n_input_wave) OR (n_short GT 0 AND n_long GT 0) THEN BEGIN
        MESSAGE, 'ERROR: Invalid wavelength(s). Please only select values in either'$
                +' the short (165 - 213) or long (245 - 292) wavelength bands.'
    ENDIF
ENDIF

IF KEYWORD_SET(short) THEN BEGIN
    wavelength = 1
ENDIF ELSE IF KEYWORD_SET(long) THEN BEGIN
    wavelength = 1000
ENDIF

; IF NOT KEYWORD_SET(slitsize) THEN slitsize = 1
;@load##########################################################################
; Load in the effective area arrays
;###############################################################################

IF N_ELEMENTS(ea_file) NE 1 THEN BEGIN
    ; ea_file = CONCAT_DIR(ROUTINE_DIR(), 'fit_eis_ea_2022-05-11.sav')
    file = find_with_def('interpol_eis_ea.pro', !path)
    break_file, file, disk, dir, name, ext
    files = FILE_SEARCH(dir, 'fit_eis_ea_*.sav')
    ea_file = files[-1]
ENDIF

RESTORE, ea_file
IF NOT KEYWORD_SET(quiet) THEN BEGIN
    PRINT, 'Loaded fit EA file: '
    PRINT, '   ', ea_file
ENDIF

IF TYPENAME(all_date_obs) EQ 'LIST' THEN BEGIN
    ; Loading old format of files
    fit_dates = all_date_obs.toarray()
    fit_easw = all_recal_easw4.toarray()
    fit_ealw = all_recal_ealw4.toarray()
    sw_wave = wsw
    lw_wave = wlw
ENDIF ELSE BEGIN
    ; Newer, cleaner files
    fit_dates = fit_ea.date_obs
    fit_easw = fit_ea.sw_ea
    fit_ealw = fit_ea.lw_ea
    sw_wave = fit_ea.sw_wave
    lw_wave = fit_ea.lw_wave
ENDELSE

ref_tai = ANYTIM2TAI(fit_dates)

IF in_tai LT ref_tai[0] THEN BEGIN
    PRINT, ' WARNING: Selected date is before the first calibrated date on '$
          +fit_ea.date_obs[0]+ '. Returning first fit calibration'
    in_tai = ref_tai[0]
ENDIF

IF in_tai GT ref_tai[-1] THEN BEGIN
    PRINT, ' WARNING: Selected date is after the last calibrated date on '$
          +fit_ea.date_obs[-1]+ '. Returning last fit calibration'
    in_tai = ref_tai[-1]
ENDIF

; Select out the desired waveband
; IF KEYWORD_SET(short) OR wavelength LT 220 THEN BEGIN
IF KEYWORD_SET(short) OR MAX(wavelength) LT 220 THEN BEGIN
    ; Short wavelength band of 170 - 211 Angstrom
    ref_ea = fit_easw
    ref_wave = sw_wave
ENDIF ELSE BEGIN
    ; Long wavelength band of 245 - 292 Angstrom
    ref_ea = fit_ealw
    ref_wave = lw_wave
ENDELSE

;@interpol######################################################################
; Interpolate the EA curve to the selected date and wavelength value
;###############################################################################

n_ref_waves = N_ELEMENTS(ref_wave)
new_ea = FLTARR(n_ref_waves)
FOR w=0, (n_ref_waves-1) DO BEGIN
    new_ea[w] = INTERPOL(ref_ea[*,w], ref_tai, in_tai)
ENDFOR

IF NOT KEYWORD_SET(short) AND NOT KEYWORD_SET(long) THEN BEGIN
    out_ea = INTERPOL(new_ea, ref_wave, wavelength, /spline)
ENDIF ELSE BEGIN
    wavelength = ref_wave
    out_ea = new_ea
ENDELSE

;@radcal########################################################################
; Compute the radcal curve or value (if requested)
;###############################################################################
IF KEYWORD_SET(radcal) THEN BEGIN
    ; Define unit conversion factors
    ;sr_pix = (!PI/(180.0*60.0^2))^2 / 1.0^2 ; 2.35044e-11 NB: asusming 1" slit!

    sr_factor = (725.0/1.496e8)^2
    ergs_to_photons = 6.626e-27*2.998e10*1.e8
    gain = 6.3 ;Hiro 6.45 8/2/2007 changed to 6.3 by Hiro & Louisa 15/4/2007
    phot_to_elec = 12398.5/3.65
    tau_sensitivity = 1894.0 ; updated from 2300.0 on 19/04/2011

    PRINT, 'Returning radcal values for converting [DN/s] to [ergs/(sr cm^2 s)]'
    PRINT, '   Note: You may still need to adjust for exposure time and slitsize.'

    ; Convert from [DN/s] to [photons/(arcsec^2 * cm^2 * s)]
    radcal = (wavelength*gain) / (out_ea*phot_to_elec)

    ; Convert to units of [ergs/(sr cm^2 s)]
    radcal = radcal * ergs_to_photons / wavelength / sr_factor ;/ slitsize

    ; radcal = radcal * 1.9866e-8 / wavelength / (!PI/(180.0*60.0^2))^2 ;GDZ code

    ; photon to ergs calculation (from eis_cal__define, which leaves out the EA)
    ; self.ergs_to_photons/self.sr_factor/lam*photons[*,*,irast] $
    ;                             /exp[irast]/abs(slitsize)

    out_ea = radcal
ENDIF

RETURN, out_ea
END