pro cg_get_UVIS_2015_calibration,u_wavelength,u_calibration,u_calibration_error,slit_width,input_spectrum,window_definition,$
bin_definition,sctime,evil_pixel_constant,noff=noff,keep_evil=keep_evil,unresolved=unresolved,channel=channel_, wdir



if n_params() le 0 then begin
Print,'Pro Get_UVIS_2015_calibration,u_wavelength,u_calibration,u_calibration_error,slit_width,input_spectrum,window_definition,bin_definition,sctime,noff=noff'
Print,'1999 Laboratory Calibration including adjustment to fit H2 model electron lamp spectrum'
Print,'Adjusted to correct for changes over time, as modeled by the observed change in the spectrum of Spica'
Print,' '
Print,'u_wavelength is the wavelength vector'
Print,' '
Print,'u_calibration is the data MULTIPLIER that converts counts/sec to kilorayleighs/pixel (or kilorayleighs/Angstrom)'
Print,' '
Print,'u_calibration_error is the 1-sigma error in u_calibration'
Print,' '
print,'slit_width selects ocultation (0), low_resolution (1), or high_resolution (2)'
print,' '
print,'input_spectrum selects calibration for continuous spectra (1) or for discrete lines (2) '
print,' '
print,'window_definition is a four element array that defines the four corners of the detector sub-array containing valid data'
print,'default is full spectral-spatial window (0,0,1023,63)'
print,' '
print,'bin_definition is a two element array that defines the spectralxspatial format'
print,'for the wavelength scale and the calibration vector'
print,' '
print,'sctime is a double-precision floats or longword integer representing the spacecraft time'
print,' '
print,'noff is a flag, nonzero to NOT apply the flat field (i.e. default is to apply flat field)'
print,' there is a preburn and postburn flat field, and sctime will determine which to use based on the time of the starburn event'
print, ' '
print,'keep_evil is a flag, nonzero to preserve the inclusion of evil pixels within the calibration'
print,'  this also prevents the application of the flat-field corrector' 
print,' '
print, 'required data files: '
print, '  EUV_1999_Lab_Cal.dat'
print, '  FLATFIELD_EUV_POSTBURN.txt'
print, '  FLATFIELD_EUV_PREBURN.txt'
print, '  FUV_1999_Lab_Cal.dat'
print, '  FLATFIELD_FUV_POSTBURN.txt'
print, '  FLATFIELD_FUV_PREBURN.txt'
print, '  spica_variability4_data.sav'
print, '  <flat-field modifier data files>
print, ' '
print, 'required routines:'
print, '  read_spica_ff_data.pro'
print, '  get_uvis_lab_sensitivity.pro'
print, '  get_uvis_ff_modifier.pro'
print, '  get_uvis_cal_modifier.pro'
return
end


channel = strupcase(channel_)

if keyword_set(wdir) eq 0 then begin
  ;
  ; wdir is the path containing all supporting data files, including: 
  ;   EUV_1999_Lab_Cal.dat, FUV_1999_Lab_Cal.dat, 
  ;   FLATFIELD_EUV_POSTBURN.txt, FLATFIELD_EUV_PREBURN.txt, FLATFIELD_FUV_POSTBURN.txt, FLATFIELD_FUV_PREBURN.txt
  ;   spica_variability4_data.sav, and flat-field modifiers (example, FUV2012_358_08_02_13_UVIS_177IC_ALPVIR001_PRIME_ff_modifier_1734944167.dat)
  ;
  wdir = '/users/Holsclaw/Cassini/UVIS/calibration_data/'
endif

if keyword_set(sctime) eq 0 then begin
  print, '*************************************'
  print, '** ERROR: you must provide sctime  **'
  print, '*************************************'
  return
endif

print,' '
print,'**************************************************************************************************************'
Print,'** 1999 FUV Channel Laboratory Calibration including adjustment to fit H2 model electron lamp spectrum      **'
Print,'**   Adjusted to correct for changes over time, as modeled by the observed change in the spectrum of Spica  **'
;*********************************************
;set up calibration matrix parameters and format
; slit width
;   This is the ratio of the low-res slit width divided by that of the chosen slit.
;   This fractional slit area will be used to scale the low-res sensitivity.
;
case channel of
'FUV': begin
  if n_params() lt 4 then slit_width=1 ; default is low resolution slit
  if (slit_width ne 1) and (slit_width ne 2) and (slit_width ne 0) then read,'*******  Select occultation (0), low resolution (1), or high resolution entrance slit ',slit_width
  if slit_width eq 1 then print,'** Calibration values for low_resolution (0.15mm wide) spectrograph slit                                    **'
  if slit_width eq 2 then print,'** Calibration values for hi_resolution (0.075mm wide) spectrograph slit                                    **'
  if slit_width eq 0 then print,'** Calibration values for occultation (0.80mm wide) spectrograph slit                                    **'
  if slit_width eq 0 then slit_width=1.5/8
end
'EUV': begin
  if n_params() lt 4 then slit_width=1 ; default is low resolution slit
  if (slit_width ne 1) and (slit_width ne 2) and (slit_width ne 0) then read,'*******  Select occultation (0), low resolution (1), or high resolution entrance slit ',slit_width
  if slit_width eq 1 then print,'** Calibration values for low_resolution (0.20mm wide) spectrograph slit **'
  if slit_width eq 2 then print,'** Calibration values for hi_resolution (0.10mm wide) spectrograph slit **'
  if slit_width eq 0 then print,'** Calibration values for occultation (0.80mm wide) spectrograph slit                                  **'
  if slit_width eq 0 then slit_width=1./4
end
endcase

; input spectrum
if n_params() lt 5 then input_spectrum=1 ; default is continuous
if (input_spectrum ne 1) and (input_spectrum ne 2) then read,'*******  Select either a continus input spectrum (1) or discrete line (2) ',input_spectrum
if input_spectrum eq 2 then Print,'** Multiply the Data (counts/pixel/sec) by u_calibration to obtain kiloRayleighs/pixel                      **'
if input_spectrum eq 1 then Print,'** Multiply the Data (counts/pixel/sec) by u_calibration to obtain kiloRayleighs/Angstrom                   **'
;
; window definition
;
if n_params() lt 6 then window_definition=[0,0,1023,63]  ; default is full detector area
if n_elements(window_definition) ne 4 then begin
    window_definition=intarr(4)
    read,'*******  input the column and row values of the upper left corner of the window ',temp1,temp2
    window_definition(0)=temp1
    window_definition(1)=temp2
    read,'*******  input the column and row values of the lower right corner of the window ',temp1,temp2
    window_definition(2)=temp1
    window_definition(3)=temp2
end
windef:
if (window_definition(2) lt window_definition(0)) or (window_definition(2) lt window_definition(0))  then begin
print,'******** faulty window definition try again'
window_definition=intarr(4)
    read,'*******  input the column and row values of the upper left corner of the window ',temp1,temp2
    window_definition(0)=temp1
    window_definition(1)=temp2
    read,'*******  input the column and row values of the lower right corner of the window ',temp1,temp2
    window_definition(2)=temp1
    window_definition(3)=temp2
end
if (window_definition(2) lt window_definition(0)) or (window_definition(2) lt window_definition(0)) then goto,windef
ul_w=[window_definition(0),window_definition(1)]
lr_w=[window_definition(2),window_definition(3)]
format='(I4)'
Print,'** Upper Left Window Corner is   ',string(ul_w(0),format=format),' x ',string(ul_w(1),format=format),'                                                                **'
Print,'** Lower Right Window Corner is  ',string(lr_w(0),format=format),' x ',string(lr_w(1),format=format),'                                                                **'
;
; bin_definition
;
if n_params() lt 7 then bin_definition=[1,1]  ;default is spectral & spatial binning=1
if n_elements(bin_definition) ne 2 then begin
    bin_definition=intarr(2)
    read,'*******  input spectral binning parameter (1=no binning) ' ,temp
    bin_definition(0)=temp
    read,'*******  input spatial binning parameter ((1=no binning) ' ,temp
    bin_definition(1)=temp
end
format='(I4)'
Print,'** Spectral x Spatial Binning is ', string(bin_definition(0),format=format),' x ',string(bin_definition(1),format=format),'                                                                **'
;
; read the sensitivity in units of (counts/second) / (kilorayleigh)
;   this is the full-slit, low-resolution, monochromatic extended source sensitivity measured in the laboratory in 1997, and updated in 1999
;
 
 
 cg_get_uvis_lab_sensitivity_temp, wcal, ucal, ucalerror, channel=channel, wdir=wdir

;
; apply slit width correction
;   The sensitivity corresponds to the low-resolution slit of both channels.
;   'slit_width' is a the width of the chosen slit (high, low, occ) divided by the low-res slit width
;
ucal=ucal/slit_width
ucalerror=ucalerror/slit_width
;
; retreive wavelength scale
;

case channel of
  'EUV': cg_e_flight_wavelength,wf
  'FUV': cg_f_flight_wavelength,wf 
endcase 
;
; interpolate the discrete-wavelength sensitivity to the full spectral range of the detector
;

ucal=interpol(ucal,wcal,wf)/60.
ucalerror=interpol(ucalerror,wcal,wf)/60.

;
; define the average pixel bandpass in angstroms (dispersion x pixel width)
;
case channel of 
'FUV': pixel_bandpass = 0.78
'EUV': pixel_bandpass = 0.605
endcase
;
; apply factor for continuous input spectrum
; sensitivity units now (counts/second) / (kilorayleigh) / (angstrom)
;
if input_spectrum eq 1 then begin
  ucal=ucal*pixel_bandpass
  ucalerror=ucalerror*pixel_bandpass
end
;
; Apply a time-dependent adjustment to the lab-measured sensitivity based on repeated observations of Spica
;
 cg_get_uvis_cal_modifier, sctime, channel, specmod, wdir
ucal = ucal * specmod
;
; construct a two-dimensional array and populate with the 1-D sensitivity
;

u_cal=fltarr(1024,64)
u_cal_error=fltarr(1024,64)
for k=2,61 do u_cal(*,k)=ucal            ; calibration per spatialpixel at the detector center row
for k=2,61 do u_cal_error(*,k)=ucalerror  ; error per spatial pixelat the detector center row
;

; apply flat field correction (default)
;
sctimeburn = 1402021717d  ; approximate spacecraft time of the starburn even (June 6, 2002)
if (keyword_set(noff) eq 0) and (keyword_set(keep_evil) eq 0) and (keyword_set(unresolved) eq 0) then begin
  ; adjust sensitivity to account for elimination of evil pixels in original calibration
  if channel eq 'FUV' then begin
    ucal = ucal / 0.91
    ucalerror = ucalerror / 0.91
  endif
  
  ;
  ; read postburn ff
  if sctime gt sctimeburn then begin
    file = wdir+'FLATFIELD_'+channel+'_POSTBURN.txt'
    cg_read_spica_ff_data,file,ff
    print, '** POSTBURN flat field applied                                                                              **'
  endif
  ; read preburn ff
  if sctime lt sctimeburn then begin
    file = wdir+'FLATFIELD_'+channel+'_PREBURN.txt'
    cg_read_spica_ff_data,file,ff
    print, '** PREBURN flat field applied                                                                               **'
  endif
  ndxnan = where( finite(ff) eq 0 )
  ;
  ; adjust flat field normalization to account for asymmetry in histogram distribution
  ;
  if channel eq 'FUV' then begin
    ff = ff * 1.05
  endif
  ;
  ; the flat field must be multiplied by the data, or equivalently, divided into the _sensitivity_ (calibration=1/sensitivity)
  ;
  case channel of 
    'FUV': begin ; GMH June 27, 2010: row 2 and row 61 in the FUV flat-field corrector appear erroneous.  For now eliminate the corrector by setting to 1
    ff[*,2] = 1.
    ff[*,61] = 1.
  end
  else:
  endcase
  ;
  ; GMH Sept 20, 2010:
  ; Set the evil pixel sensitivity to some constant fraction of its neighboring good pixels
  ;
  ;evil_pixel_constant = 0.5
  ;evil_pixel_constant = 0.
  ;u_cal[ndxnan] = u_cal[ndxnan] * evil_pixel_constant
  ;ff[ndxnan] = 1.
  ;
  ; apply the flat-field corrector 
  ;
 
  u_cal = u_cal / ff
  u_cal_error = u_cal_error / ff
  ;
  ; apply the flat-field modifier
  ;

  cg_get_uvis_ff_modifier, sctime, channel, arrmod, wdir
  arrmod[ndxnan] = 1.
  u_cal = u_cal * arrmod
  u_cal_error = u_cal_error * arrmod
endif

;
; do not apply flat-field corrector
;
if (keyword_set(noff) eq 1) and (keyword_set(keep_evil) eq 0) and (keyword_set(unresolved) eq 0) then begin
  ; adjust sensitivity to account for elimination of evil pixels in original calibration
  if channel eq 'FUV' then begin
    ucal = ucal / 0.91
    ucalerror = ucalerror / 0.91
  endif
  ;
  ; if no flat field is applied, set all evil pixels to NaN as the calibration now only applies to good pixels
  file = wdir+'FLATFIELD_'+channel+'_POSTBURN.txt'
  cg_read_spica_ff_data,file,ff
  ndxnan = where( finite(ff) eq 0 )
  u_cal[ndxnan] = !values.f_nan
  u_cal_error[ndxnan] = !values.f_nan
endif
;
; Unresolved source
;
if (keyword_set(unresolved) eq 1) and (keyword_set(keep_evil) eq 0) then begin
  ; adjust sensitivity to account for elimination of evil pixels in original calibration
  if channel eq 'FUV' then begin
    ucal = ucal / 0.91
    ucalerror = ucalerror / 0.91
  endif
  ;
  file = wdir+'FLATFIELD_'+channel+'_PREBURN.txt'
  cg_read_spica_ff_data,file,ff_pre
  file = wdir+'FLATFIELD_'+channel+'_POSTBURN.txt'
  cg_read_spica_ff_data,file,ff_post
  
  sens_burn = 1.
  if sctime gt sctimeburn then sens_burn = ff_pre / ff_post
  
  u_cal = u_cal * sens_burn
  u_cal_error = u_cal_error * sens_burn
  ;
  ; apply the flat-field modifier
  ;
 
  cg_get_uvis_ff_modifier, sctime, channel, arrmod, wdir
  ;arrmod[ndxnan] = 1.
  u_cal = u_cal * arrmod
  u_cal_error = u_cal_error * arrmod
endif

;
; if we keep the evil pixels, then don't force the calibration only apply to good pixels
; 
if (keyword_set(keep_evil) eq 1) then begin
endif

print,'**************************************************************************************************************'
print,' '
;
; extract the correct window padded out to integral binning
;
spectral_bin=bin_definition(0)
spatial_bin=bin_definition(1)
;get size of unbinned window
width=lr_w(0)-ul_w(0)+1
height=lr_w(1)-ul_w(1)+1
; pad it out to integral multiples of the binning parameters
width_pad=spectral_bin*ceil(1.*width/spectral_bin)
height_pad=spatial_bin*ceil(1.*height/spatial_bin)
; get the size of the final returned array partial bins at the ends are included
width_final=ceil(1.*width/spectral_bin)
height_final=ceil(1.*height/spatial_bin)
;
; create binned and padded versions of the wavelength vector, sensitivity, and sensitivity uncertainty
;

  ;12/31/2017; Code added by Todd Bradley to interpolate calibration factor values so that there are no NaN values
  for zzz = 0,63 do begin
    u_cal_tempp = u_cal(*,zzz)
    d=where(finite(u_cal_tempp),complement = c)
    u_cal_interpol = interpol(u_cal_tempp(d),d,findgen(1024))
    ;if total(u_cal_interpol) gt 0. then stop
    u_cal_interpol(c) = u_cal_interpol(c)*.1
    u_cal(*,zzz) = u_cal_interpol
    ;if total(u_cal_interpol) gt 0. then stop
  endfor



wave_temp=fltarr(width_pad)
u_cal_temp=fltarr(width_pad,height_pad)
u_err_temp=fltarr(width_pad,height_pad)
wave_temp(0:width-1)=wf(ul_w(0):lr_w(0))
u_cal_temp(0:width-1,0:height-1)=u_cal(ul_w(0):lr_w(0),ul_w(1):lr_w(1))
u_err_temp(0:width-1,0:height-1)=u_cal_error(ul_w(0):lr_w(0),ul_w(1):lr_w(1))
u_wavelength=fltarr(width_final)
u_calibration=fltarr(width_final,height_final)
u_calibration_error=fltarr(width_final,height_final)


 

;
; average wavelength values within each bin
;
for k=0,width_final-1 do u_wavelength(k)=mean(wave_temp(k*spectral_bin:(k+1)*spectral_bin-1))
;
; sum the sensitivity
; sum the sensitivity uncertainty in quadrature
;
;the lines of code in the k loop below that are commented out was Todd's attempt to deal with the nan issue

for k=0,height_final-1 do for l=0,width_final-1 do begin
    ;toddtemp = u_cal_temp(l*spectral_bin:(l+1)*spectral_bin-1,k*spatial_bin:(k+1)*spatial_bin-1); commented out on 6/3/2017
    ;d=finite(toddtemp,/nan); commented out 6/3/2017
    ;dd=where(d eq 1); commented out 6/3/2017
    ;ddd=where(dd eq 0); commented out 6/3/2017
    ;dddd = where(dd eq spectral_bin-1); commented out 6/3/2017
    ;if total(dd) ne -1 then begin; commented out 6/3/2017
    ;  toddtemp(dd) = (toddtemp(dd-1) + toddtemp(dd+1))/2.; commented out 6/3/2017
    ;  stop; commented out 6/3/2017
    ;  if ddd ne -1 then toddtemp(0) = toddtemp(1); commented out 6/3/2017
    ;  if dddd ne -1 then toddtemp(spectral_bin-1) = toddtemp(spectral_bin-2); commented out 6/3/2017
    ;  u_cal_temp(l*spectral_bin:(l+1)*spectral_bin-1,k*spatial_bin:(k+1)*spatial_bin-1) = toddtemp ; commented out on 6/3/2017
    ;end  ; commented out 6/3/2017
    ;if k eq 54 then stop; commented out on 6/3/2017
    u_calibration(l,k)=total(u_cal_temp(l*spectral_bin:(l+1)*spectral_bin-1,k*spatial_bin:(k+1)*spatial_bin-1))
    u_calibration_error(l,k)=total(u_err_temp(l*spectral_bin:(l+1)*spectral_bin-1,k*spatial_bin:(k+1)*spatial_bin-1))/sqrt(spatial_bin*spectral_bin)
;if k eq 10 then stop
end



cal_temp=(u_calibration eq 0)*1.e10+u_calibration
err_temp=(u_calibration eq 0)*1.e10+u_calibration_error


;
; invert the sensitivity to obtain the 'calibration' factor
;   units will be (kR) / (counts/second) for input_spectrum=0
;   units will be (kR/angstrom) / (counts/second) for input_spectrum=1
;
u_calibration=1./cal_temp
u_calibration_error=err_temp/cal_temp/cal_temp
u_calibration=(u_calibration gt 1.e-9)*u_calibration
u_calibration_error=(u_calibration_error gt 1.e-9)*u_calibration_error

end
return
end