FUNCTION iris_get_calib, wvl, time_obs, photons=photons, units=units, $
                         exp_time=exp_time, ybin=ybin, perang=perang, $
                         eff_area=eff_area
;+
; NAME:
;      IRIS_GET_CALIB
;
; PURPOSE:
;      This routine returns an array containing the multiplicative
;      calibration factor corresponding to the input wavelength
;      array.
;
; CATEGORY:
;      IRIS; radiometric calibration.
;
; CALLING SEQUENCE:
;      Result = IRIS_GET_CALIB( WVL, TIME_OBS )
;
; INPUTS:
;      Wvl:    A 1D wavelength array in angstroms.
;      Time_Obs:   A string giving a time in a standard IDL format. 
;
; OPTIONAL INPUTS:
;      Exp_Time:  The exposure time in seconds. If not set, then 1
;                 second is assumed.
;      Ybin:      Binning of pixels in Y-direction. If not set, then 1
;                 is assumed.
;	
; KEYWORD PARAMETERS:
;      PHOTONS:  The output units will be given in photons rather than
;                ergs.
;      PERANG:   If set, then the intensity is returned in "per
;                angstrom" units rather than "per pixel" units.
;
; OUTPUTS:
;      An array of same size as WVL giving the factor by which an
;      intensity array in DN is multiplied to give the intensity in
;      the units indicated by UNITS (see below).
;
;      If a problem is found, then the value -1 is returned.
;
; OPTIONAL OUTPUTS:
;      Units:   A string giving the units of the output array.
;      Eff_Area: An array of same size as WVL containing the effective
;                area in cm^2.
;
; EXAMPLE:
;      IDL> wd=iris_getwindata(file,1402.77)
;      IDL> cal=iris_get_calib(wd.wvl,wd.hdr.date_obs)
;
; MODIFICATION HISTORY:
;      Ver.1, 31-Jan-2018, Peter Young
;      Ver.2, 6-Mar-2018, Peter Young
;          Added /perang keyword.
;      Ver.3, 27-May-2019, Peter Young
;          Added eff_area optional output.
;-


IF n_params() LT 2 THEN BEGIN
  print,'Use:  IDL> cal = iris_get_calib( wvl, time_obs [, /photons, units=, exp_time=, ybin='
  print,'                                 /perang ])'
  return,-1
ENDIF 

;
; Read the IRIS calibration file.
;
iresp=iris_get_response(time_obs)
IF NOT keyword_set(quiet) THEN print,'% IRIS_GET_CALIB: instrument response version is '+iresp.version


IF n_elements(exp_time) EQ 0 THEN exp_time=1.0
IF n_elements(ybin) EQ 0 THEN ybin=1

;
; Convert pixel size to solid angle. The slit is set as 0.33",
; and a Y-pixel is 0.16635". 
;
pix_size=(2*!pi/(360.*3600.))^2 * 0.33 * 0.16635 * float(ybin)

;
; Work out whether spectrum is NUV or FUV.
;
midwvl=mean(wvl)
IF midwvl LE 1500. THEN j=0 ELSE j=1

n=n_elements(wvl)

;
; Wavelengths are given in nm in the iresp structure. I have to do
; some trickery to deal with wavelengths where the eff. area is zero. 
;
yi=fltarr(n)
x=iresp.lambda*10.   ; convert from nm to angstroms
y=iresp.area_sg[*,j]   ; effective area in cm^2

;
; Restrict the effective area to the wavelength range of the input
; spectrum.
;
k=where(x LT min(wvl))
getmax=max(x[k],imax)
i0=k[imax]
;
k=where(x GT max(wvl))
getmin=min(x[k],imin)
i1=k[imin]
;
x=x[i0:i1]
y=y[i0:i1]

;
; Further filter out wavelengths with zero eff. area
;
k=where(y NE 0,nk)
IF nk GT 0 THEN BEGIN
  x=x[k]
  y=y[k]
ENDIF

;
; Now to do spline fit to the remaining good eff. area values.
;
y2=spl_init(x,y)
;
k=where(wvl GE min(x) AND wvl LE max(x))
yi[k]=spl_interp(x,y,y2,wvl[k])

dn2phot_sg=iresp.dn2phot_sg[j]


;
; The array scale_factor contains the multiplicative factors that
; convert from DN to the output intensity units. The effective area
; (eff_area) is an optional output.
;
scale_factor=fltarr(n)
eff_area=fltarr(n)
eff_area[k]=yi[k]

;
; Compute conversion factor to phot cm^-2 s^-1 sr^-1 units.
; Note that only pixels defined by k will have an intensity; all
; others will be missing.
;
k=where(yi NE 0.)
scale_factor[k]=dn2phot_sg/yi[k]/float(exp_time)/pix_size

;
; Note: the /perang keyword may generate funny results if an unusual
; wavelength array is specified. It is meant to be used if WVL
; corresponds to the wavelength of an IRIS data window.
;
IF keyword_set(perang) THEN BEGIN
  ang=fltarr(n)
  ang[1:n-1]=wvl[1:n-1]-wvl[0:n-2]
  ang[0]=wvl[1]-wvl[0]
  ang_units='Angstrom^-1'
ENDIF ELSE BEGIN
  ang=fltarr(n)+1.0
  ang_units='pixel^-1'
ENDELSE 

;
; Convert to erg if necessary.
;
IF NOT keyword_set(photons) THEN BEGIN
  scale_factor[k]=scale_factor[k]*1.986e-8/wvl[k]/ang
  units='erg cm^-2 s^-1 sr^-1 '+ang_units
ENDIF ELSE BEGIN
  units='photon cm^-2 s^-1 sr^-1 '+ang_units
ENDELSE

return,scale_factor

END