@iris_calib_spectrum.pro
;+
; NAME:
;   IRIS_FLUX_RADCAL
;
; PURPOSE:
;   Wrapper for the absolute radiometric calibration of an IRIS spectrum. 
;   It returns either the calibration factor (/onlyfactor) or the calibration factor 
;   and the datas calibrated (/factor_and_data). Several options for the units output 
;   are given. Some other useful variables are given (e.g. wavelengths for each 
;   calibrated data)
;
; CATEGORY:
;   Data analysis
;
; CALLING SEQUENCE:
;  output_structure = iris_flux_radcal(initial_time, final_time)
;  (See EXAMPLES for other options)
;
; INPUTS:
;       t0: Initial time for IRIS data 
;       t1: Ending time for IRIS data
;
; OPTIONAL INPUTS as a keyword:
;       local: local filename. If giving, t0 and t1 values are ignored.
;
; KEYWORD PARAMETERS
;       units: four options are considered. For the sake of simplicity, 
;              we refer to them as:
;              0) lambda (cgs): erg cm^-2 s^-1 sr^-1 Å^-1    <<< default option 
;              1) lambda (SI): W m^-2 sr^-1 nm^-1
;              2) Hz (cgs): erg cm^-2 s^-1 sr^-1 Hz^-1
;              3) Hz (SI): W m^-2 sr^-1 Hz^-1
;       onlyfactor: it returns only the wavelength-depending calibration factor    <<< default option 
;       factor_and_data: it return the wavelength-depending calibration factor and the data properly calibrated.
;       mosaic: it looks for the closest level2 fits file to the mosaic date (YYYYMMDD), then, calculates the calibration from that file info. 
;       verbose: it prompts about units and/or factor/data options if units and/or factor/data keywords are not given.
;
; DESCRIPTION:
;       The wrapper is making use of IRIS_GET_RESPONSE developed by JP. Wuelser (LMSAL),
;       and IRIS_CALIB_SPECTRUM developed G. Vissers (University of Stockholm). 
;        
; EXAMPLES:
;         t0 = '2015-08-18 00:00:00' & t1 = '2015-08-19 00:00:00'
;         - Interactived mode to calibrate selected IRIS data between times t0 and t1. Selection of data, units and 
;           output mode are requested by prompt:
;              iris_datacal = iris_flux_radcal(t0,t1)     
;         - Interactived mode to calibrate selected IRIS data between times t0 and t1, units passed by keyword:
;              iris_datacal = iris_flux_radcal(t0,t1, units = 2)     
;         - Calibration of local file. Initial and final time values are ignored. Units and output mode are requested
;           by the prompt:
;              iris_datacal = iris_flux_radcal('', '', local='/my_iris_l2_files/iris_l2_fits_file.fits')
;         - Automatic calibration of local file. Initial and final time values are ignored:
;              iris_datacal = iris_flux_radcal('', '', local='/my_iris_l2_files/iris_l2_fits_file.fits', $
;                                              units= 2, /factor_and_data)
;         - Example of the message showing the output structure variables when units = 2 and /factor_and_data
;              ********************************************************************************
;              RESULTS
;
;              Factor, WL, and calibrated data are avalaible as output_name.variable_name,
;              being variable_name:
;         
;              FILENAME
;              VER_RF2
;              BLD_VERS
;              Flux_units 
;              Info_source_file
;              data_Hz_2814
;              data_Hz_2832
;              data_Hz_C_II_1336
;              data_Hz_Mg_II_k_2796
;              data_Hz_O_I_1356
;              data_Hz_Si_IV_1394
;              data_Hz_Si_IV_1403
;              factor_Hz_2814
;              factor_Hz_2832
;              factor_Hz_C_II_1336
;              factor_Hz_Mg_II_k_2796
;              factor_Hz_O_I_1356
;              factor_Hz_Si_IV_1394
;              factor_Hz_Si_IV_1403
;              wl_2814
;              wl_2832
;              wl_C_II_1336
;              wl_Mg_II_k_2796
;              wl_O_I_1356
;              wl_Si_IV_1394
;              wl_Si_IV_1403
;              ********************************************************************************
;
; MODIFICATION HISTORY:
;
; 2017 Oct 23: A. Sainz Dalda (asainz.solarphysics@gmail.com) (First version)
; 2018 March:  A. Sainz Dalda: verbose and mosaic options are now available.
;                              Flux_units is given as a variable in the output structure.
;                              Therefore, factor_ and data_ are now given in Flux_units,
;                              and there is no distinction in the variable name. 
;
;   $Id$
;-
;

function iris_flux_radcal, t0, t1, local=local, units=units, $ 
         onlyfactor=onlyfactor, factor_and_data=factor_and_data, mosaic=mosaic, mgii_only=mgii_only, verbose=verbose

if keyword_set(local) eq 0  then begin
   ;;; Get/show information of IRIS data available between t0 and t1 data time 

   tl = iris_time2timeline(t0,t1)
   info = get_infox(tl, tag_names(tl))
   dim_info = size(info)
   for j=0, dim_info[1]-1 do print, j,' --- ', info[j]
   
   
   if keyword_set(mosaic) eq 0 then begin 
      ;;; Selection of an IRIS Level data set 
      read, ans, prompt='Select an observation data set... '
      print, 'Selected data ', ans, ' ---> ',  info[ans]
      print
      info_sel = strsplit(info[ans], " ", /extract)
      ds = 'level2'
      l2_fits = iris_ssw_l2_query(info_sel[0], info_sel[1], /fits, ds=ds)
      if l2_fits ne -1 then begin
         l2_fits = l2_fits[where(strpos(l2_fits, '_raster_') ne -1)] ;;; We just care abour raster fits files.
         dim_l2 = size(l2_fits)
         for j=0, dim_l2[1]-1 do print, j,' --- ', l2_fits[j]
         ;;; Selection of a file 
         print
         read, ans, prompt='Select an observation data file... '
         print
      endif  	 
   endif else begin 
      where_Mg = (where(strpos(tl.DESCRIPTION,'Mg') ne -1, count))
      if where_Mg[0] ne -1 then begin 
         l2_fits = - 1
	 j = 0
         while l2_fits[0] eq -1 and j le count do begin
	       ans = where_Mg[j] 
               info_sel = strsplit(info[ans], " ", /extract)
               ds = 'level2'
               l2_fits = iris_ssw_l2_query(info_sel[0], info_sel[1], /fits, ds=ds)
               l2_fits = l2_fits[where(strpos(l2_fits, '_raster_') ne -1)] ;;; We just care abour raster fits files.
               dim_l2 = size(l2_fits)
	       j = j + 1       
         endwhile 
	 ans = 0
      endif 	 
   endelse

   ans = fix(ans)
   if ans gt dim_l2[1]-1 or ans lt 0 then return, -1
   print, 'Getting data ', ans, ' ---> ',  l2_fits[ans]
   url = l2_fits[ans]
   data_home = './'
   home = data_home
   iris_get_l2files, url, ds = ds, data_home = data_home
   
   ;;; Getting information from the selected data filename
   ;;; Copied & pasted from iris_get_l2files.pro
   tmp = str_sep(url,ds)
   tmp = tmp[1]
   tmp = str_sep(tmp, '/')
   nt = n_elements(tmp)
   year = tmp[nt-5]
   month = tmp[nt-4]
   day = tmp[nt-3]
   obs_folder = tmp[nt-2]
   file = tmp[nt-1]
   outdir = home + year + '/' + month+ '/' + day + '/' + obs_folder + '/'
   ;;; Copied & pasted from iris_get_l2files.pro
   local_file = outdir+file
   download_flag = 'File info obtained from file web/remotely downloaded after a time based search: '+local_file
endif else begin
   local_file = file_search(local)
   filename = local_file
   download_flag = 'File info obtained from file loaded from local directory: '+local_file
   print
   print, 'Working on local file... '+local_file
   print
endelse 

;;; Reading data locally 
d = iris_load(local_file)
d->show_lines
nwin = d->getnwin()
radcal = 0
radcal = dictionary('filename', local_file)
read_iris_l2, local_file, index, data, wave = name_win


;;; How to calibrate the data
ans_factor = 0 ;;; lambda (cgs) is the default option
if keyword_set(units) eq 0  and keyword_set(verbose) eq 1 then begin
   print
   ans_factor = -1
   print
   while ans_factor lt 0 or ans_factor gt 3 do begin 
        ;read, ans_factor, prompt = 'In what units do you want your data to be calibrated? (0 = lambda (cgs), 1 = lambda (SI), 2 = Hz (cgs), 3 = Hz (SI)): '
        print
        print, 'In what units do you want your data to be calibrated?' 
        print
	print, '0) lambda (cgs): erg cm^-2 s^-1 sr^-1 Å^-1'  
        print, '1) lambda (SI): W m^-2 sr^-1 nm^-1'
        print, '2) Hz (cgs): erg cm^-2 s^-1 sr^-1 Hz^-1'
        print, '3) Hz (SI): W m^-2 sr^-1 Hz^-1'
        print
        read, ans_factor, prompt = 'Select an option: '
   endwhile	
   print
endif else begin
   if keyword_set(units) eq 1 then begin
      ans_factor =  units	   
      print
      while ans_factor lt 0 or ans_factor gt 3 do begin 
        ; read, ans_factor, prompt = 'In what units do you want your data to be calibrated? (0 = lambda (cgs), 1 = lambda (SI), 2 = Hz (cgs), 3 = Hz (SI)): '
        print
        print, 'In what units do you want your data to be calibrated?' 
        print
	print, '0) lambda (cgs): erg cm^-2 s^-1 sr^-1 Å^-1'  
        print, '1) lambda (SI): W m^-2 sr^-1 nm^-1'
        print, '2) Hz (cgs): erg cm^-2 s^-1 sr^-1 Hz^-1'
        print, '3) Hz (SI): W m^-2 sr^-1 Hz^-1'
        print
        read, ans_factor, prompt = 'Select an option: '
     endwhile	
   endif 
endelse

;;; Return a dictionary variable including calibration factor and WL or calibration, WL and data calibrated
;;; as an IDL structure
ans_data = 1 ;;; onlyfactor is default option 
if (keyword_set(onlyfactor) or keyword_set(factor_and_data)) eq 0  and keyword_set(verbose) eq 1 then begin
   print
   read, ans_data, prompt = 'Do you want to calibrate data? (1 = No, just the factor array, 2 = YES) '
   print
endif else begin 
   if keyword_set(onlyfactor) then ans_data = 1
   if keyword_set(factor_and_data) then ans_data = 2
endelse 

;;; Calculating calibrated data for each window available in the selected file
radcal = orderedhash('FILENAME',local_file)
radcal['VER_RF2'] = index[0].VER_RF2
radcal['BLD_VERS']= index[0].BLD_VERS
radcal['Info_source_file']= download_flag
ini_win = 0
fin_win = nwin-1
if keyword_set(mgii_only) then begin
      nwin = 1 
      iwindx = d->getwindx(2796.20)
      ini_win = iwindx[0]
      fin_win = iwindx[0]
endif
for j = ini_win, fin_win do begin
   if (ans_factor ge 0 and ans_factor le 3) and (ans_data eq 1 or ans_data eq 2)  then begin 
      print, 'Recovering data for window... ', d->getline_id(j)
      lambda =  d->getlam(j)
      dim_lambda =  size(lambda)
      name_win = d->getline_id(j)
      exptime = d->getexp(iwin=j)
      ;print, name_win+'a'
      aux_name_win  = strcompress(name_win)
      if strpos(aux_name_win," ", /reverse_search) eq strlen(aux_name_win)-1 then  aux_name_win = strmid(aux_name_win, 0, strlen(aux_name_win)-1)
      name_win_ok = (strrep(strcompress(aux_name_win), ' ', '_'))[0]
      name_wl_ok = 'wl_'+name_win_ok 
      ;radcal[name_wl_ok]  = lambda
      radcal[name_wl_ok]  = lambda
      ; name_key_factor = ['factor_l_', 'factor_l_SI_', 'factor_Hz_', 'factor_Hz_SI_'] + name_win_ok
      name_key_factor = ['factor_'] + name_win_ok
      name_key_dn2photon = ['dn2photon_'] + name_win_ok
      name_key_Ephoton= ['Ephoton_'] + name_win_ok
      ; name_key_texp_used = ['texp_used_'] + name_win_ok
      ; name_key_texp_mean = ['texp_mean_'] + name_win_ok
      name_key_texp_arr = ['texp_arr_'] + name_win_ok
      ; name_key_texp_stddev = ['texp_stddev_'] + name_win_ok
      ; name_key_texp_minmax = ['texp_minmax_'] + name_win_ok
      units_info = ['lambda (cgs): erg cm^-2 s^-1 sr^-1 Å^-1', 'lambda (SI): W m^-2 sr^-1 nm^-1', $
                    'Hz (cgs): erg cm^-2 s^-1 sr^-1 Hz^-1', 'Hz (SI): W m^-2 sr^-1 Hz^-1']

      case ans_factor of
           0: begin & i_nu = 0 & si = 0 & end
           1: begin & i_nu = 0 & si = 1 & end
           2: begin & i_nu = 1 & si = 0 & end
           3: begin & i_nu = 1 & si = 1 & end
      end
      print, 'Calculating calibration values with...'
      print, 'obs. time = ',(index.date_obs)[0]
      print, 'binxy     = ',(index.sumspatn)[0]
      print, 'binl      = ',(index.sumsptrn)[0]
      w_t = where(exptime ne 0.0) 
      exptime_ok = mean((exptime)[w_t])
      minmax_exptime_ok = minmax((exptime)[w_t])
      stddev_exptime_ok = stddev((exptime)[w_t])
      exptime_used = exptime_ok
      exptime_used = 1.0
      print, 'exp. time = ', exptime_used ;(index.exptimen)[0]
      ; Equals to 1 to recover the actual value through the corresponding keyword
      factor = 1
      dn2photon = 1
      Ephoton = 1
      ic = iris_calib_spectrum(indgen(dim_lambda[1]), lambda, time = (index.date_obs)[0], $
                               binxy = (index.sumspatn)[0], binl = (index.sumsptrn)[0], t_exp=exptime_used, $
                               factor = factor, dn2photon = dn2photon, Ephoton = Ephoton, i_nu = i_nu, si = si, verbose=1)
		       
      radcal[name_key_factor[0]]= factor
      radcal[name_key_dn2photon[0]]  = dn2photon
      radcal[name_key_Ephoton[0]]  = Ephoton
      ; radcal[name_key_texp_used[0]]  = exptime_used
      ; radcal[name_key_texp_mean[0]]  = exptime_ok
      radcal[name_key_texp_arr[0]]  = exptime 
      ; radcal[name_key_texp_minmax[0]]  = minmax_exptime_ok
      ; radcal[name_key_texp_stddev[0]]  = stddev_exptime_ok

      ;;; Exposure time correction
      ;apply_indvidual_texp = 0
      ;if (minmax_exptime_ok)[1] - (minmax_exptime_ok)[0] ge 1.0 then  apply_indvidual_texp = 1.
      ;apply_indvidual_texp = 1

      if (fix(ans_data) eq 2) then begin 
         read_iris_l2, local_file, index, data, wave = name_win
         dimdata = size(data)
         data_new  =  data*0
         ; name_key_data = ['data_l_', 'data_l_SI_', 'data_Hz_', 'data_Hz_SI_'] + name_win_ok
         name_key_data = ['data_'] + name_win_ok
         for s = 0, dimdata[3] - 1 do begin  
          for k = 0, dimdata[2] - 1 do begin  
              ;if apply_indvidual_texp eq 0 then begin		  
                 ;data_new[*,k,s] = data[*,k,s]*radcal[name_key_factor[ans_factor]]*(exptime_used/exptime_ok)
              ;endif else begin 	 
	         data_new[*,k,s] = data[*,k,s]*radcal[name_key_factor[0]]*(exptime_used/exptime[s])
	      ;endelse  	 
          endfor
         endfor
         radcal[name_key_data[0]] = data_new
      endif
   endif
endfor
radcal['Flux_units'] =  units_info[ans_factor]
out = -1 

;;; Assigning some extra variables and keys. Printing all of them
if (ans_factor ge 0 and ans_factor le 3) and (ans_data eq 1 or ans_data eq 2)  then begin 
   keys = radcal.Keys()
   keys = keys.ToArray()
   keys = keys[sort(keys)]
   dimk = size(keys)
   keys_ok = keys
   keys_ok[0] = 'FILENAME'
   keys_ok[1] = 'VER_RF2'
   keys_ok[2] = 'BLD_VERS'
   keys_ok[3] = 'Flux_units'
   keys_ok[4:dimk[1]-1] = keys[where(keys ne 'FILENAME' and keys ne 'VER_RF2' and keys ne 'BLD_VERS' and keys ne 'Flux_units')]
   keys = keys_ok
   print
   print,'********************************************************************************'
   print, 'RESULTS'
   print
   if ans_data eq 1 then begin  
         print, 'Factor and WL are avalaible as output_name.variable_name,'
         print, 'being variable_name:'
         print
         for k = 0, dimk[1]-1 do print, keys[k] 
   endif else begin 
         print, 'Factor, WL, and calibrated data are avalaible as output_name.variable_name,'
         print, 'being variable_name:'
         print
         for k = 0, dimk[1]-1 do print, keys[k] 
   endelse
   print
   print, 'Factor will convert original DNs to Flux in '+radcal['Flux_units']
   print
   print,'********************************************************************************'
   print
   out = radcal[keys]
   out = out.tostruct()
endif    

return, out  

end