FUNCTION iris_calib_spec, spec, photons=photons, quiet=quiet, perang=perang

;+
; NAME:
;	IRIS_CALIB_SPEC()
;
; PURPOSE:
;       This routine takes a 1D IRIS spectrum in the format output by
;       iris_sum_spec.pro and calibrates it into erg cm^-2 s^-1 sr^-1
;       pixel^-1 units.
;
; CATEGORY:
;	IRIS; calibration.
;
; CALLING SEQUENCE:
;	Result = IRIS_CALIB_SPEC( Spec )
;
; INPUTS:
;	Spec:	A structure in the format output by IRIS_SUM_SPEC.
;
; KEYWORD PARAMETERS:
;	PHOTONS: If set, then the output units are photon cm^-2 s^-1
;	         sr^-1 pixel^-1.
;       QUIET:   If set, then no information messages are printed.
;       PERANG:  If set, then the spectrum is given in "per-Angstrom"
;                units rather than "per-pixel".
;
; OUTPUTS:
;       A spectrum structure with the same format as the input
;       structure, but the following tags are changed:
;
;         .int   Converted to intensity units.
;         .err   Converted to intensity units.
;         .int_units  Updated to reflect new units.
;         .time_stamp Updated to new spectrum formation time.
;
;       If a problem is found, then -1 is returned.
;
; EXAMPLE:
;       IDL> spec=iris_sum_spec(file,xpix=10,ypix=200,wavel=1402)
;       IDL> calspec=iris_calib_spec(spec)
;
; MODIFICATION HISTORY:
;       Ver.1, 26-Aug-2016, Peter Young
;       Ver.2, 21-Sep-2016, Peter Young
;         I now handle the wavelengths where the effective area is
;         zero by just setting the output intensity to missing (I was
;         getting some strange spline points previously).
;       Ver.3, 5-Dec-2016, Peter Young
;         The calibration version number is no longer hard-coded.
;       Ver.4, 8-Mar-2018, Peter Young
;         Added /perang keyword.
;-


IF n_params() EQ 0 THEN BEGIN
  print,'Use:  IDL> iris_calib_spec, spec [, /photon, /perang, /quiet]'
  return,-1
ENDIF 

;
; Copy spec to outspec, and reset the intensity and error arrays to
; missing. 
;
outspec=spec
outspec.int=spec.missing
outspec.err=spec.missing


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

;
; Get exposure time and pixel size from input spectrum.
;
texp=spec.exp_time
pix_size=(2*!pi/(360.*3600.))^2 * spec.xscale * spec.yscale

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

n=n_elements(spec.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(spec.wvl))
getmax=max(x[k],imax)
i0=k[imax]
;
k=where(x GT max(spec.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(spec.wvl GE min(x) AND spec.wvl LE max(x))
yi[k]=spl_interp(x,y,y2,spec.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.
;
n=n_elements(spec.int)
scale_factor=fltarr(n)

;
; If /perang is set, then create the array dwvl containing the size of
; each wavelength bin. 
;
IF keyword_set(perang) THEN BEGIN
  dwvl=fltarr(n)
  dwvl[1:*]=spec.wvl[1:*]-spec.wvl[0:n-2]
  dwvl[0]=spec.wvl[1]-spec.wvl[0]
  wunits='Angstrom^-1'
ENDIF ELSE BEGIN
  dwvl=fltarr(n)+1.0
  wunits='pixel^-1'
ENDELSE 


;
; 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(spec.int NE spec.missing AND yi NE 0.)
scale_factor[k]=dn2phot_sg/yi[k]/pix_size/texp

;
; Convert to erg if necessary.
;
IF NOT keyword_set(photons) THEN BEGIN
  scale_factor[k]=scale_factor[k]*1.986e-8/spec.wvl[k]
  outspec.int[k]=scale_factor[k]*spec.int[k]/dwvl[k]
  outspec.err[k]=scale_factor[k]*spec.err[k]/dwvl[k]
  outspec.int_units='erg cm^-2 s^-1 sr^-1 '+wunits
ENDIF ELSE BEGIN
  outspec.int[k]=scale_factor[k]*spec.int[k]/dwvl[k]
  outspec.err[k]=scale_factor[k]*spec.err[k]/dwvl[k]
  outspec.int_units='photon cm^-2 s^-1 sr^-1 '+wunits
ENDELSE 

outspec.time_stamp=systime()

return,outspec

END