;+
; NAME:
;       eis_prep
;
; PURPOSE:
;
; CATEGORY:
;       Hansteen/Wikstol Data analysis SW
;
; CALLING SEQUENCE:
;
; INPUTS:
;
;       input:Name of FITS file with EIS data to be calibrated. If the file 
;             is already open a data object can be sent instead of a filename.
;
; KEYWORD PARAMETERS:
;
;                               -- Dark current data --
;       default: Use the default file containing dark current data.
;
;       dc_file: Optional name of the dark current file.
;
;                                   -- Saving --
;       save: Save the data and the error to files. The data is saved
;             to "eis_l1_YYYYMMDD_HHMMSS.fits" The error is saved to
;             "eis_er_YYYYMMDD_HHMMSS.fits". The default is to save to
;             the current directory.
;
;       outdir: Optional directory for storing the fits files.
;
;       datfile: Optional name of the data fits file.
;
;       errfile: Optional name of the error fits file.
;
;                              -- Calibration options --
;       retain: Retain negative pixel values after dark current subtraction.
;
;       nocr: Don't remove cosmic rays (despiking).
;
;       nohp: Don't remove hot pixels.
;
;       nowp: Don't remove warm pixels.
;
;       nodp: Don't remove dusty pixels.
;
;       noabs: Don't perform absolute calibration.
;
;       photons: Return absolutely calibrated data in units of photons rather
;                than ergs.
;
;       refill: Fill missing mixels and errors with average values of the 
;               neighbouring valid pixels according to the method of P.Young.
;               Only up to three consecutive pixels in the same column are
;               filled. The method is intended to recover scientific 
;               information mostly lost because of the warm pixels.
;
;                                 -- Despiking --
;       neighbours: Mark the neighbours of bad pixels as bad in EIS_DESPIKE.
;
;       n_iter:
;
;       _extra:
;
;                                  -- General --
;       quiet: Don't issue xwindow messages.
;
;       cleanup: Cleanup data if set.
;
; OUTPUTS:
;
;       If save is set, a data file and an error file are written: their names 
;       are "eis_l1_YYYYMMDD_HHMMSS.fits" and "eis_er_YYYYMMDD_HHMMSS.fits" by
;       default. These names can be superseded by setting the datfile and 
;       errfile keywords.`
;
; CALLS:
;
;
; COMMON BLOCKS:
;
;
; PROCEDURE:
;
;
; RESTRICTIONS:
;
;
; MODIFICATION HISTORY:
;       22-Jun-2006  Oivind Wikstol - Version 1.0
;       12-Feb-2007  Viggo Hansteen - Slight change in structure, added 
;                                     despiking. Still to do are error 
;                                     estimates, bad pixel marking
;                                     and quality control (ie debugging)!
;       28-Mar-2007  Harry P Warren - New despiking keywords, new output 
;                                     keywords, fix to DC subtraction
;       02-Apr-2007  Viggo Hansteen - Force conversion to float even if only 
;                                     doing DC subtraction
;       10-Oct-2007  A. Gardini     - Modified the naxis1 computation;
;                                     set nohp=nohp and noff=noff
;       10-Oct-2007  Harry P Warren - fixed no_neighbours despiking keyword
;       19-Oct-2007  Harry P Warren - added delvarx to help with memory 
;                                     problems
;       23-Oct-2007  A. Gardini     - Marking of the dusty pixels
;       20-Nov-2007  A. Gardini     - CR removal performed after hot/dusty
;                                     pixels correction
;       28-Nov-2007  A. Gardini     - Added a check on the number of
;                                     missing pixels during CR removal
;       30-Nov-2007  Viggo Hansteen - Inserted getexp method in place of
;                                     difference between t1 and t2
;        7-Jan-2008  A. Gardini     - Check if exposure time = 0 in abs cal.
;        8-Jan-2008  A. Gardini     - Modified the CAL_ keywords in l1 files.
;       10-Jan-2008  A. Gardini     - Corrections in the HP calibration
;        1-Feb-2008  Mike Marsh	    - Added photons keyword, errors 
;                                     corrected based on photon stats.
;        2-Feb-2008  Mike Marsh     - Set photon calibration status.
;        7-Feb-2008  Peter Young    - Changed parameter calling eis_despike
;       12-Feb-2008  A. Gardini     - Checked that the data file level is 0, 
;                                     and that ABS did not precede CR.
;       19-Feb-2008  A. Gardini     - Added warm pixels calibration.
;                                     Uploaded on 1-Apr-2004.
;        4-Apr-2008  A. Gardini     - Moved CR removal before HP. Set xbox=1.
;        6-May-2008  A. Gardini     - Removed the texp=0 check in abs cal.
;       16-May-2008  Peter Young    - Previously eis_prep set any pixels with
;                                     photon numbers < 1 after dc subtraction
;                                     to be missing as the photon errors on 
;                                     such pixels did not make sense. 
;                                     For these pixels the error is now set i
;                                     to the read noise and they are no 
;                                     longer flagged as missing.
;       19-May-2008  A. Gardini     - Implemented the option above according
;                                     to the new keyword retain.
;        9-Jul-2008  A. Gardini     - Set the new values of the read out 
;                                     noise according to the quadrant.
;       20-Feb-2009  A. Gardini     - Moved the error computation in DN from 
;                                     the CR removal to the ABS calibration.
;        2-Apr-2009  A. Gardini     - Fixed datfile and errfile. 
;        6-Apr-2009  A. Gardini     - Added the refill option from P. Young. 
;       28-Apr-2009  A. Gardini     - Fixed an 1 pixel error in the  pixel-
;                                     wvl conversion for the CCD-A (LW) in
;                                     the absolute calibration.
;       30-Apr-2009  A. Gardini     - Set keywords to eis_cal.
;        5-May-2009  A. Gardini     - Modified the refill option to deal 
;                                     with retained negative pixels.
;        6-May-2009  A. Gardini     - Fixed the infinite loops in hp/wp.
;       24-Jun-2009  A. Gardini     - Fixed the photon unit tag.
;       26-Jun-2009  A. Gardini     - Roll back of the previous change.
;       16-Apr-2010  T. Fredvik     - Added keyword hkwavecorr: if set correct
;                                     for orbital variations of the line
;                                     centre due to temperature variations,
;                                     using Kamio-san's code.
;       30-Jun-2010  V. Hansteen    - Added keyword
;                                     "/correct_sensitivity" which
;                                     should correct for CCD
;                                     degradation as a funciton of time. 
;       11-Aug-2010  T. Fredvik     - Added check for missing data when
;                                     calculating error arrays
;       19-Apr-2011  V. Hansteen    - Mark ZERO DN as "missing" data
;       14-Jan-2014  V. Hansteen/H. Warren - New (correct) error
;                                            estimate for non-absolute
;                                            calibrated data.
;       22-Jan-2015  H. Warren      - Removed stop
;       29-Jan-2015  H. Warren      - corrected "correct_sensistivity" bug on line 615
;       17-Oct-2019  H. Warren      - corrected i_not_missing bug
;       17-Oct-2019  P. Young       - flag 2048 DN points as missing;
;                                     modified input behavior so that
;                                     if input is a string then it is
;                                     not returned as an object.
;       21-Oct-2019  P. Young       - modified 2048DN check to require
;                                     entire column have this value;
;                                     also pixel values are set to 0.
;       06-Dec-2019   P. Young      - fixed 2048 check for nexp_prp>1;
;                                     renamed as eis_prep_original
;                                     (from eis_prep).
;
; $Id: eis_prep.pro 748 2011-04-19 10:03:15Z viggoh $
;
;-
pro eis_prep_original,input, $
;             default=default, dc_file=dc_file, ff_file=ff_file,$
             default=default, dc_file=dc_file,$
             save=save,outdir=outdir,datfile=datfile,errfile=errfile,$
             retain=retain,$
             refill=refill,$
             nocr=nocr,nohp=nohp,nowp=nowp,nodp=nodp,noabs=noabs,$
             photons=photons,$
             correct_sensitivity=correct_sensitivity,$
;            noff=noff,$
             quiet=quiet,cleanup=cleanup,$
             neighbours=neighbours,n_iter=n_iter, hkwavecorr=hkwavecorr,$
             _extra=ex
;
  ver_no='$Revision: 748 $'
  if n_params() ne 1 then begin
    message,'eis_prep,data [data,file], /default,/save, $',/info
    message,'    /nohp,/nodp,/nocr,/noabs,/photons,/correct_sensitivity, $',/info
    message,'    /nowp, $',/info
;   message,'    /noff, $',/info
    message,'    dc_file=dc_file, $',/info
;   message,'    ff_file=ff_file, $',/info
    message,'    outdir=outdir,datfile=datfile,errfile=errfile, $',/info
    message,'    _extra=ex,/quiet',/info
    return
  endif
  if n_elements(default) eq 0 then default=0
  if n_elements(save) eq 0 then save=0
  if n_elements(retain) eq 0 then retain=0
  if n_elements(nocr) eq 0 then nocr=0
  if n_elements(nohp) eq 0 then nohp=0
  if n_elements(nowp) eq 0 then nowp=0
  if n_elements(nodp) eq 0 then nodp=0
; if n_elements(noff) eq 0 then noff=0
  if n_elements(noabs) eq 0 then noabs=0
  if n_elements(photons) eq 0 then photons=0
  if n_elements(correct_sensitivity) eq 0 then correct_sensitivity=0
  if n_elements(quiet) eq 0 then quiet=0
  IF n_elements(hkwavecorr) EQ 0 THEN hkwavecorr = 0
  if default and (n_elements(dc_file) ne 0) then begin
    message,'WARNING: /default overrides dc_file',/info
    dc_file=''
  endif
; if default and (n_elements(dc_file) ne 0 or n_elements(ff_file) ne 0) then begin
;   message,'WARNING: /default overrides dc_file and/or ff_file',/info
;   dc_file=''
;   ff_file=''
; endif
  if n_elements(dc_file) eq 0 then dc_file=''
; if n_elements(ff_file) eq 0 then ff_file=''

; open data file if necessary
;
; PRY 17-Oct-2019: introduced "file_type" below and renamed "data" to
;                  "input". 
;  - allows routine to remember if input was given as a string or
;    object (see end of routine).
;
  case datatype(input) of
    'STR': begin
           file_type=0b
           data=obj_new('eis_data',input,datasource='fits')
           fitslev = data-> getfitslev()
           if fitslev ne 0 then begin
             message,'eis_prep runs on level 0 files only',/info
             return
             end 
           end
    'OBJ': BEGIN 
      file_type=1b
      data=input
    END 
    else : begin
      message,'data parameter must be either eis_data object or filename',/info
      return
           end
  endcase

  fitslev=1

  data-> setfits_reformat, 'ver_rf1', 'EIS_prep '+ver_no
  get_utc, date, /ccsds
  home_inst = data-> gethome_inst()
  data-> setfitslev,fitslev

  if fitslev eq 1 then begin
    data-> setfits_reformat, 'ver_rf1', 'EIS_prep Ver '+ver_no
    data-> setfits_reformat, 'date_rf1', date
    data-> setfits_reformat, 'orig_rf1', home_inst
  endif

; read in DC data: either manually pick, given in dc_file, or default set 
; in eis_cal object

  if not default and dc_file eq '' then begin
    dc_dir=getenv('EIS_CAL_DATA')+path_sep()+'dc'+path_sep()
    title = 'Pick Dark Current file'
    dc_file = dialog_pickfile(path = dc_dir, filter = '*.idlsave', title = title)
  endif

;; read in Flat Field (FF) data (FF file must be manually picked)

; if not default and ff_file eq '' then begin
;   ff_dir=getenv('EIS_CAL_DATA')+path_sep()+'ff'+path_sep()
;   title = 'Pick Flat Field'
;   ff_file = dialog_pickfile(path = ff_dir, filter = '*.idlsave', title = title)
; endif

;; define calibration object with dc_file and ff_file
;  cal = obj_new('eis_cal', dc_file = dc_file, ff_file = ff_file, /gen_aeff)
; define calibration object with dc_file
  time_obs=(data->gethdr())->getdate_obs()
  slit_ind=data->getinfo('slit_ind')
  cal = obj_new('eis_cal', dc_file = dc_file, /gen_aeff,$
                time_obs = time_obs, slit_ind = slit_ind)
  calstat=''
;
; CHECK for SATURATED DATA and for data set to ZERO
;
; PRY 18-Oct-2019: added check for 2048 DN points
;
  missing=-100.
  if not (data->getcalstat()).dc then begin
    for iwin=0,data->getnwin()-1 do begin
      wd=float(data->getvar(iwin))
      (data->getcal())->setupvar,wd,iwin ; set up error array
      err=wd*0.0
     ;
      i=where(wd eq 16383,ni)
      if ni ne 0 then err[i]=missing
     ;
      i=where(wd eq 0,ni)
      if ni ne 0 then err[i]=missing
     ;
     ; PRY, 21-Oct-2019
     ; The following handles 2048 DN pixels. It assumes that an entire
     ; column has 2048 values, so I check the min and max of each column,
     ; and if they're both 2048 then I set the column to
     ; missing. Note that I replace the WD 
     ; value with 0, which gets rid of the bright column usually seen
     ; for these pixels.
     ;
      nexp=data->getnexp()
      FOR i=0,nexp-1 DO BEGIN
        wdmin=min(reform(wd[*,*,i]),dim=2)
        wdmax=max(reform(wd[*,*,i]),dim=2)
        k=where(wdmin EQ 2048 AND wdmax EQ 2048,nk)
        IF nk GT 0 THEN BEGIN
          err[k,*,i]=-50.
        ENDIF 
      ENDFOR
      k=where(err EQ -50.,nk)
      IF nk GT 0 THEN BEGIN
        print,'% EIS_PREP: window '+trim(iwin)+' has '+trim(nk)+' DN2048 pixels. Set to missing.'
        err[k]=missing
        wd[k]=0
      ENDIF 

      data->setvar,wd,iwin
      data->seterr,err,iwin
      delvarx,wd,err
    endfor
  ENDIF

;
; DARK CURRENT subtraction
;
  if not (data->getcalstat()).dc then begin
    if retain then $
      message = ['Dark current subtraction...     ','Retained values less than 0.'] $
    else message = ['Dark current subtraction...     ']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    for iwin=0,data->getnwin()-1 do begin
      wd=float(cal->dc(data,iwin))
      if not retain then wd = wd > 0
      data->setvar, wd, iwin
      err=data->geterr(iwin)
      if not retain then begin
        i=where(wd lt 1,ni)
        if ni ne 0 then err[i]=missing
      endif
      data->seterr,err,iwin
      delvarx,wd,err
    endfor
    data->setcalstat, 'dc', 1  ; set calibration status to 'true' for dc
    if retain then data->setcalstat, 'retain', 1  
  endif else begin
    message,'Dark current subtraction done previously',/info
  endelse
  if (data->getcalstat()).dc then calstat='DC' else begin
    message,'Dark current subtraction must be performed in order ot proceed',/info
    return
  end
;
; COSMIC RAY removal 
;
  if not (data->getcalstat()).cr and not (data->getcalstat()).abs and $
    not nocr then begin
    message = ['Running new_spike...','...this will take QUITE a while']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    for iwin=0,data->getnwin()-1 do begin
      nexp_prp=data->getnexp_prp()
      wd=data->getvar(iwin)
      if keyword_set(neighbours) then no_neighbours = 0 else no_neighbours = 1
      eis_despike,wd,wdds,nexp_prp,swtch=swtch,no_neighbours=no_neighbours,n_iter=n_iter,_extra=ex
      data->setvar, wdds, iwin
      err=data->geterr(iwin)
      i_miss=where(swtch gt 0.5,ni)
      if ni ne 0 then err[i_miss]=missing
; ERROR estimation moved inside the ABSOLUTE CALIBRATION
;      i_not_miss=where(err gt missing,ni)
;      if ni ne 0 then err[i_not_miss]=sqrt(wd[i_not_miss])
      data->seterr,err,iwin

      delvarx,wd,err,wdds,swtch
      ;;print,"  MEMORY = "+trim(double(memory(/current))/1.0E+6)
    endfor
    data->setcalstat,'cr', 1 ; set calibration status to 'true' for cr
  endif else begin
    if (data->getcalstat()).cr then message,'Cosmic ray removal done previously',/info else if (data->getcalstat()).abs then message,'Cosmic ray removal cannot be performed if the absolute calibration has been already done',/info
  endelse
  if (data->getcalstat()).cr then calstat=calstat+' CR'
;
; MARK HOT PIXELS and FILL in DATA with MEDIAN
;
  if not (data->getcalstat()).hp and not nohp then begin
    message = ['Hot pixel correction...']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    xbox=1 & ybox=3
    for iwin=0,data->getnwin()-1 do begin
      wd=data->getvar(iwin)
      err=data->geterr(iwin)
      hp=cal->hot_pixels(data,iwin) ; stores in hp the part of the hp matrix coincident with the window
      for iexp=0,data->getnexp()-1 do begin ;loop over exposures
        i=where(hp eq 1,ni) ; ni == count
        if ni ne 0 then begin
          wdum=wd[*,*,iexp]
          edum=err[*,*,iexp]
          edum[i]=missing
          err[*,*,iexp]=edum
; fill in hot pixels with median (remembering where they are in err)
          wdum[i]=missing
          ct = ni
          repeat begin
; remember that other pixels, not just hp, are already marked as missing AG
            ct_old=ct
            newmed = fmedian(wdum,xbox,ybox,missing=missing,/only_missing)
            wdum[i] = newmed[i]
            newmiss = where(newmed[i] eq missing,ct)
          end until (ct eq 0) or (ct ge ct_old)
          if ((ct ge ct_old) and not quiet) then message,'WARNING: '+string(ct)+' hot pixels impossible to refill',/info
          wd[*,*,iexp]=wdum
        endif
      endfor
      data->seterr,err,iwin
      data->setvar,wd,iwin
    endfor
    data->setcalstat, 'hp', 1  ; set calibration status to 'true' for hp
  endif else begin
    if (data->getcalstat()).hp then message,'Hot pixels marked previously',/info
  endelse
  if (data->getcalstat()).hp then calstat=calstat+' HP'
;
; MARK WARM PIXELS and FILL in DATA with MEDIAN
; Notice that it is now considered one warm pixel map for all the exposures 
; of a given window, because usually they have the same duration. 
; In some cases however, a sit-and-stare observations consists
; of exposures of different duration. 
; In these cases, different warm pixels maps should possibly have been 
; considered.
;
  if not (data->getcalstat()).wp and not nowp then begin
    message = ['Warm pixel correction...']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    xbox=1 & ybox=3
    for iwin=0,data->getnwin()-1 do begin
      wd=data->getvar(iwin)
      err=data->geterr(iwin)
; stores in wp the part of the wp matrix coincident with the window
      wp=cal->warm_pixels(data,iwin) 
      if wp[0] ne -1 then begin
        for iexp=0,data->getnexp()-1 do begin ;loop over exposures
          i=where(wp eq 1,ni) ; ni == count
          if ni ne 0 then begin
            wdum=wd[*,*,iexp]
            edum=err[*,*,iexp]
            edum[i]=missing
            err[*,*,iexp]=edum
; fill in warm pixels with median (remembering where they are in err)
;           wdum[i]=0. ; TEMP!
            wdum[i]=missing
            ct=ni
            repeat begin
              ct_old=ct
              newmed = fmedian(wdum,xbox,ybox,missing=missing,/only_missing)
              wdum[i] = newmed[i]
              newmiss = where(newmed[i] eq missing,ct)
            end until (ct eq 0) or (ct ge ct_old)
            if ((ct ge ct_old) and not quiet) then message,'WARNING: '+string(ct)+' warm pixels impossible to refill',/info
            wd[*,*,iexp]=wdum
          endif
        endfor
        data->seterr,err,iwin
        data->setvar,wd,iwin
      end
    endfor
    data->setcalstat, 'wp', 1  ; set calibration status to 'true' for wp
  endif else begin
    if (data->getcalstat()).wp then message,'Warm pixels marked previously',/info
  endelse
  if (data->getcalstat()).wp then calstat=calstat+' WP'
;
; MARK DUSTY PIXELS
;
  if not (data->getcalstat()).dp and not nodp then begin
    message = ['Dusty pixel correction...']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    for iwin=0,data->getnwin()-1 do begin
      err=data->geterr(iwin) ; returns error window iwin
      dp=cal->dusty_pixels(data,iwin)
      for iexp=0,data->getnexp()-1 do begin
        i=where(dp eq 1,ni)
        if ni ne 0 then begin
          edum=err[*,*,iexp]
          edum[i]=missing
          err[*,*,iexp]=edum
        endif
      endfor
      data->seterr,err,iwin
    endfor
    data->setcalstat, 'dp', 1  ; set calibration status to 'true' for dp
  endif else begin
    if (data->getcalstat()).dp then message,'Dusty pixels marked previously',/info
  endelse
  if (data->getcalstat()).dp then calstat=calstat+' DP'
;
; FLAT FIELD correction
;
;  if not (data->getcalstat()).ff and not noff then begin
;    message = ['Flat field correction...']
;    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
;    for iwin=0,data->getnwin()-1 do begin
;      wd0=data->getvar(iwin) > 1.0
;      wd=cal->ff(data,iwin)
;      data->setvar, wd, iwin
;      err=data->geterr(iwin)
;      i=where(wd le 0.0,ni)
;      if ni ne 0 then err[i]=missing
;      i_not_miss=where(err ne missing)
;      err[i_not_miss]=wd[i_not_miss]*err[i_not_miss]/wd0[i_not_miss]
;      data->seterr,err,iwin
;    endfor
;    data->setcalstat, 'ff', 1  ; set calibration status to 'true' for ff
;  endif else begin
;    if (data->getcalstat()).ff then message,'Flat field done previously',/info
;  endelse
;  if (data->getcalstat()).ff then calstat=calstat+' FF'
;


;
; ABSOLUTE CALIBRATION and ERROR estimation
;
  if not (data->getcalstat()).abs and not noabs then begin
    message = ['Absolute calibration...']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    exp=data->getexp()
    t0 = where(exp eq 0.,ct0) 
    if (ct0 gt 0) then begin
      message,'Some exposure times are equal to 0. Set to 1 to complete calibration',/info
      exp[t0] = 1.
    endif 
    case data->getinfo('slit_ind') of
        0: begin
           slitsize=1.0*1.0
           if not photons then begin
             data->setunit,['erg/cm!e2!N/s/sr/'+string("305B),'','']
           endif else begin
             data->setunit,['Photon-Events']
           endelse
           end
        2: begin
           slitsize=1.0*2.0
           if not photons then begin
             data->setunit,['erg/cm!e2!N/s/sr/'+string("305B),'','']
           endif else begin
             data->setunit,['Photon-Events']
           endelse
           end
        else: begin
          slitsize=-1.0
           if not photons then begin
             data->setunit,['erg/cm!e2!N/s/sr','','']
           endif else begin
             data->setunit,['Photon-Events']
           endelse
          message,'This is slot data - using slitsize of 1.0 and dispersion of 1.0 to compute absolute calibration',/info
            end
    endcase

    for iwin=0,data->getnwin()-1 do begin
      data->getwin,iwin,dn,pos
      err=data->geterr(iwin)
      i_not_miss=where(err ne missing, count)

      if count eq 0 then begin
        intensity = err
        data->setvar, intensity, iwin
        print, ' ! no good data, skipping this window ' + trim(iwin)
        continue
      endif

      ;; --- read out noise definition
      quad=cal->quadrant(pos)
      case quad of
      0: begin
         read_noise = 2.26 ;; in DN
         end
      1: begin
         read_noise = 2.29 ;; in DN
         end
      2: begin
         read_noise = 2.37 ;; in DN
         end
      3: begin
         read_noise = 2.24 ;; in DN
         end
      endcase

      if pos[0] gt 2147 then begin
        detector='A'
        pos[0]=pos[0]-2148
        lambda=((cal->getlambda()).scale_a)[pos[0]:pos[0]+pos[1]-1]
      endif else begin
        detector='B'
        lambda=((cal->getlambda()).scale_b)[pos[0]:pos[0]+pos[1]-1]
      endelse
      aeff = cal->getaeff(detector)
      effarea=interpol(aeff.eff,aeff.lambda,lambda,/spline)
      sub=where(lambda lt aeff.lambda[0] or lambda gt aeff.lambda[n_elements(aeff.lambda)-1])
      if sub[0] ne -1 then effarea[sub]=1.e-20
; intensity in photons, for error calculation
      intensity_ph = cal->dn_to_photon(lambda,float(dn))	

      data->setcalstat,'sens',0
      if not photons then begin
        days_since_launch=0.0
        if correct_sensitivity then days_since_launch=data->days_since_launch()
        if days_since_launch ne 0.0 then begin
          if correct_sensitivity eq 1 then begin
             action=' Correcting ' 
             (data->getcal())->setcorrected_sensitivity,1
             data->setcalstat, 'sens', 1
          endif ; else  begin
;            action=' Uncorrecting '
;            days_since_launch=-1*days_since_launch
;             (data->getcal())->setcorrected_sensitivity,0
;          endelse          
          message,string(days_since_launch,format='(f7.1)')+' days since lauch.' $
                  +action+'sensitivity with tau = ' $
                 +string((data->getcal())->gettau_sensitivity(),format='(f7.1)')+' days.',/info
        endif
;        print,(data->getcal())->sensitivity_loss(days_since_launch)
        intensity = cal->photon_to_ergs(lambda*effarea,intensity_ph, $
                           exp=exp,slitsize=slitsize,days_since_launch=days_since_launch) 
      endif else intensity=intensity_ph ;leave intensity in photons if photons keyword set

      ;; --- add in the read noise
      rnoise_dn = fltarr(size(dn,/dimensions))+read_noise  ;sets read noise array in DN
      rnoise_ph = cal->dn_to_photon(lambda,float(rnoise_dn))  ;read noise array in photons
      ;; PRY 12-Feb-08 - modified line below by adding "> 0."
      err[i_not_miss] = sqrt(intensity_ph[i_not_miss] > 0. +rnoise_ph^2)  ;photon variance = photon counts

      data->setvar, intensity,iwin

      ;;  Changed method of calculating error using the photon_to_ergs method
      ;;  Old method follows below (commented out). PRY
      if not photons then begin
        i_miss=where(err EQ missing)   
        err = cal->photon_to_ergs(lambda*effarea,err, $
                     exp=exp,slitsize=slitsize,days_since_launch=days_since_launch)
        err[i_miss]=missing   
      endif
;     err[i_not_miss]=intensity[i_not_miss]*err[i_not_miss]/intensity_ph[i_not_miss]	;multiply by fractional error in photons

      data->seterr,err,iwin

      delvarx,err,intensity,i_not_miss
      ;print,"  MEMORY = "+trim(double(memory(/current))/1.0E+6)
    endfor
    data->setcalstat, 'abs', 1  ; set calibration status to 'true' for ff
    if photons then data->setcalstat, 'phot', 1  ;Mike Marsh 2/2/08 set photon calibration status to 'true' if keyword set
  endif else begin
    if (data->getcalstat()).abs then begin
      message,'Absolute calibration done previously',/info
      if correct_sensitivity eq -1 and (data->getcal())->getcorrected_sensitivity() then begin
        message,'Uncorrecting for sensitivity with tau = ' $
           +string((data->getcal())->gettau_sensitivity(),format='(f7.1)')+' days.',/info
        for iwin=0,data->getnwin()-1 do begin
          data->getwin,iwin,intensity,pos
          err=data->geterr(iwin)
;
          days_since_launch=data->days_since_launch()
          sens_loss=(data->getcal())->sensitivity_loss(abs(days_since_launch))
          intensity=intensity*sens_loss
          err=err*sens_loss
          data->setvar,intensity,iwin
          data->seterr,err,iwin
;
          (data->getcal())->setcorrected_sensitivity,0
          data->setcalstat, 'sens', 0
;
          delvarx,err,intensity
        endfor    
      endif
    endif else begin ; error computation in DN
      for iwin=0,data->getnwin()-1 do begin
        data->getwin,iwin,dn,pos
        err=data->geterr(iwin)
        i_not_miss=where(err ne missing, nmissing)
        ;; --- read out noise definition
        quad=cal->quadrant(pos)
        case quad of
          0: begin
            read_noise = 2.26 ;; in DN
          end
          1: begin
            read_noise = 2.29 ;; in DN
          end
          2: begin
            read_noise = 2.37 ;; in DN
          end
          3: begin
            read_noise = 2.24 ;; in DN
          end
        endcase
        ;; [Start HPW 13-JAN-2016]
        ;; get wavelength 
        if pos[0] gt 2147 then begin
          detector='A'
          pos[0]=pos[0]-2148
          lambda=((cal->getlambda()).scale_a)[pos[0]:pos[0]+pos[1]-1]
        endif else begin
          detector='B'
          lambda=((cal->getlambda()).scale_b)[pos[0]:pos[0]+pos[1]-1]
        endelse
        ;; sets read noise array in DN
        rnoise_dn = fltarr(size(dn,/dimensions))+read_noise
        IF nmissing GT 0 THEN begin
          err[i_not_miss] = sqrt(dn[i_not_miss] > 0. + rnoise_dn^2)*sqrt((cal->photon_to_dn(lambda,1))[i_not_miss])
        endif ELSE err[*] = -100
        data->seterr,err,iwin
        delvarx,err,dn,i_not_miss
      endfor
      ;; [Stop HPW 13-JAN-2016]      
    endelse
  endelse
  if (data->getcalstat()).abs then calstat=calstat+' ABS'
  if (data->getcalstat()).sens then calstat=calstat+' SENS'
  if (data->getcalstat()).phot then calstat=calstat+' PHOT'
  if (data->getcalstat()).retain then calstat=calstat+' RETAIN' 
;
; UPDATE data header
;
  hdr = data->gethdr()
  main_hdr = hdr->getmain_hdr()
  bte_hdr = hdr-> getbte_hdr()
  fits_reformat=data->getfits_reformat()

  fxaddpar, main_hdr, 'DATA_LEV', (data-> getfitslev()), $
                    'FITS Level (0, 1, 2)'
  fxaddpar, main_hdr, 'DATE_RF1', fits_reformat.date_rf1, $
                    'Date and time of Level 1 reformat'
  fxaddpar, main_hdr, 'ORIG_RF1', fits_reformat.orig_rf1, $
                    'Institution where Level 1 reformat was done'
  fxaddpar, main_hdr, 'VER_RF1', fits_reformat.ver_rf1, $
                    'Fits Level 1 reformatter version no.'
  fxaddpar, bte_hdr, 'DATA_LEV', (data-> getfitslev()), $
                    'FITS Level (0, 1, 2)'
  fxaddpar, bte_hdr, 'DATE_RF1', fits_reformat.date_rf1, $
                    'Date and time of Level 1 reformat'
  fxaddpar, bte_hdr, 'ORIG_RF1', fits_reformat.orig_rf1, $
                    'Institution where Level 1 reformat was done'
  fxaddpar, bte_hdr, 'VER_RF1', fits_reformat.ver_rf1, $
                    'Fits Level 1 reformatter version no.'

  calflag=data->getcalstat()
  fxaddpar,main_hdr,'CALSTAT',calstat,'Calibration status', after = 'CMIRR  '
  fxaddpar,bte_hdr,'CALSTAT', calstat,'Calibration status', after = 'CMIRR  '
  fxaddpar,main_hdr,'CAL_DC',calflag.dc,after='CALSTAT'
  fxaddpar,bte_hdr,'CAL_DC',calflag.dc,after='CALSTAT'
  fxaddpar,main_hdr,'CAL_HP',calflag.hp,after='CAL_DC'
  fxaddpar,bte_hdr,'CAL_HP',calflag.hp,after='CAL_DC'
  fxaddpar,main_hdr,'CAL_DP',calflag.dp,after='CAL_HP'
  fxaddpar,bte_hdr,'CAL_DP',calflag.dp,after='CAL_HP'
  fxaddpar,main_hdr,'CAL_WP',calflag.dp,after='CAL_DP'
  fxaddpar,bte_hdr,'CAL_WP',calflag.dp,after='CAL_DP'
; fxaddpar,main_hdr,'CAL_FF',calflag.ff
; fxaddpar,bte_hdr,'CAL_FF',calflag.ff
  fxaddpar,main_hdr,'CAL_CR',calflag.cr,after='CAL_WP'
  fxaddpar,bte_hdr,'CAL_CR',calflag.cr,after='CAL_WP'
; fxaddpar,main_hdr,'CAL_WVL',calflag.wvl
; fxaddpar,bte_hdr,'CAL_WVL',calflag.wvl
  fxaddpar,main_hdr,'CAL_ABS',calflag.abs,after='CAL_CR'
  fxaddpar,bte_hdr,'CAL_ABS',calflag.abs,after='CAL_CR'
  fxaddpar,main_hdr,'CAL_PHOT',calflag.phot,after='CAL_ABS'
  fxaddpar,bte_hdr,'CAL_PHOT',calflag.phot,after='CAL_ABS'
  fxaddpar,main_hdr,'CAL_SENS',calflag.sens,after='CAL_PHOT'
  fxaddpar,bte_hdr,'CAL_SENS',calflag.sens,after='CAL_PHOT'
  fxaddpar,main_hdr,'CAL_RETA',calflag.retain,after='CAL_SENS'
  fxaddpar,bte_hdr,'CAL_RETA',calflag.retain,after='CAL_SENS'
  sxdelpar, main_hdr, ['CAL_FF','CAL_WVL']
  sxdelpar, bte_hdr, ['CAL_FF','CAL_WVL']
;
  if calflag.sens then begin
    fxaddpar,main_hdr,'TAU_SENS',(data->gettau_sensitivity()),$
      'Sensitivity decrease e-folding time (days)',after='CAL_RETA'
    fxaddpar,bte_hdr,'TAU_SENS',(data->gettau_sensitivity()),$
      'Sensitivity decrease e-folding time (days)',after='CAL_RETA'
  endif
; max and min values in line windows have changed.
; TDMIN and TDMAX must be set again
  nwin = data-> getnwin()
  naxis1=0l
  for iwin=0, nwin-1 do begin
    data-> getwin, iwin, wd, pos
    wdmin = min(wd)  ; minimum value in line window
    wdmax = max(wd)  ; maximum value in line window
    param1 = 'TDMIN'+strtrim(string(iwin+1), 2)
    param2 = 'TDMAX'+strtrim(string(iwin+1), 2)
    fxaddpar, bte_hdr, param1, wdmin
    fxaddpar, bte_hdr, param2, wdmax
    param1 = 'TFORM'+strtrim(string(iwin+1), 2)
    winsize=ulong((data->getxw())[iwin])*ulong((data->getyw())[iwin])
    winstring=strtrim(string(winsize),2)
    case datatype(wd) of
      'FLO': begin
        naxis1=winsize*4l+naxis1
        fxaddpar,bte_hdr,param1,winstring+'E','Real*4 (floating point)'
             end
      'INT': begin
        naxis1=winsize*2l+naxis1
        fxaddpar,bte_hdr,param1,winstring+'I','Integer*2 (short integer)'
             end
      else : message,'Unexpected datatype for wd, all bets are off!',/info
    endcase
  endfor
; naxis1=naxis1+73l   ; NB this is a HACK!!!! should count number of bytes added in auxilary data
  naxis1=naxis1+65l
  if ptr_valid((data->getaux_data()).mhc_hz_t10) then naxis1=naxis1+8l
  if ptr_valid((data->getaux_data()).v_sat) then naxis1=naxis1+12l

                      ; by reading the bte_hdr as it exists.
  fxaddpar,main_hdr,'NAXIS1',naxis1,'Number of bytes per row'
  fxaddpar,bte_hdr,'NAXIS1',naxis1,'Number of bytes per row'
  hdr-> setmain_hdr, main_hdr
  hdr-> setbte_hdr, bte_hdr
;
; REFILLING missing pixels according to the method discussed by PRY in
; http://msslxr.mssl.ucl.ac.uk:8080/eiswiki/Wiki.jsp?page=WarmPixelStudy
;
  if keyword_set(refill) then begin
    message = ['Refilling missing pixel',$
               'and errors: up to three',$
               'pixels in a column are filled.']
    if not quiet then xmessage,message,wbase=wbase,font='helvetica'
    for iwin=0,data->getnwin()-1 do begin
      data->getwin,iwin,wd,pos
      err=data->geterr(iwin)
; computation of wvl and ea
      if pos[0] gt 2147 then begin
        detector='A'
        pos[0]=pos[0]-2148
        lambda=((cal->getlambda()).scale_a)[pos[0]:pos[0]+pos[1]-1]
      endif else begin
        detector='B'
        lambda=((cal->getlambda()).scale_b)[pos[0]:pos[0]+pos[1]-1]
      endelse
      aeff = cal->getaeff(detector)
      effarea=interpol(aeff.eff,aeff.lambda,lambda,/spline)
; the routine spl_interp returns slightly different results.
      y2 = spl_init(aeff.lambda,aeff.eff)
      effarea=spl_interp(aeff.lambda,aeff.eff,y2,lambda)
      sub=where(lambda lt aeff.lambda[0] or lambda gt aeff.lambda[n_elements(aeff.lambda)-1])
      if sub[0] ne -1 then effarea[sub]=1.e-20
;-----
      sz=size(wd)
      nl=sz[1]
      ny=sz[2]
      nt=sz[3]
      wvl=fltarr(nl,ny,nt)
      ea=fltarr(nl,ny,nt)
      for k=0,nl-1 do begin
        wvl[k,*,*]=lambda[k]
        ea[k,*,*]=effarea[k]
      endfor
;-----
      i=where(err eq missing,nmiss)
      if nmiss gt 0 then begin
        frac=[1.0,1.2,1.2,1.3,1.3]
;
        chck=bytarr(nl,ny,nt)
        miss_map=bytarr(nl,ny,nt)
        miss_map(i)=1b
        tofill_map=miss_map
        x=fltarr(nl,ny,nt)
        y=fltarr(nl,ny,nt)
        x=wd*ea
        y=err^2*ea^2*wvl
        j=where(err ne missing and wd gt 0.,nj)
        f=poly_fit(x[j],y[j],1)
;
        wdaux=fltarr(nl,ny,nt)
        wdout=fltarr(nl,ny,nt)
        wdout=wd
        errout=err
;       
        chck[*]=1b
        wdaux[*,0:ny-2,*]=wd[*,1:ny-1,*]
        chck[*,0:ny-2,*]=miss_map[*,1:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[4]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac_val
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,1:ny-1,*]=wd[*,0:ny-2,*]
        chck[*,1:ny-1,*]=miss_map[*,0:ny-2,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[4]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[4]
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,2:ny-3,*]=0.5*(wd[*,0:ny-5,*]+wd[*,4:ny-1,*])
        chck[*,2:ny-3,*]=miss_map[*,0:ny-5,*]+miss_map[*,4:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[3]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[3]
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,3:ny-2,*]=(wd[*,0:ny-5,*]+3.5*wd[*,4:ny-1,*])/4.5
        chck[*,3:ny-2,*]=miss_map[*,0:ny-5,*]+miss_map[*,4:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[2]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[2]
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,1:ny-4,*]=(3.5*wd[*,0:ny-5,*]+wd[*,4:ny-1,*])/4.5
        chck[*,1:ny-4,*]=miss_map[*,0:ny-5,*]+miss_map[*,4:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[2]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[2]
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,2:ny-2,*]=(wd[*,0:ny-4,*]+2.*wd[*,3:ny-1,*])/3.
        chck[*,2:ny-2,*]=miss_map[*,0:ny-4,*]+miss_map[*,3:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[1]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[1]
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,1:ny-3,*]=(2.*wd[*,0:ny-4,*]+wd[*,3:ny-1,*])/3.
        chck[*,1:ny-3,*]=miss_map[*,0:ny-4,*]+miss_map[*,3:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[1]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])*frac[1]
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
;
        chck[*]=1b
        wdaux[*,1:ny-2,*]=0.5*(wd[*,0:ny-3,*]+wd[*,2:ny-1,*])
        chck[*,1:ny-2,*]=miss_map[*,0:ny-3,*]+miss_map[*,2:ny-1,*]
        k=where(chck[i] eq 0b,nk)
        frac_val=frac[0]
        if nk ne 0 then begin
          wdout[i[k]]=wdaux[i[k]]
;         errout[i[k]]=sqrt(wdout[i[k]]*f[1]+f[0])
          errout[i[k]]=sqrt( 1.0/wvl[i[k]]/ea[i[k]]^2 * $
                       (f[1]*(wdout[i[k]]>0.)*ea[i[k]]+f[0]) )*frac_val
          tofill_map[i[k]]=0b
        endif
; For all remaining pixels, set intensity and error to be missing.
        k=where(tofill_map eq 1b,nk)
        if nk ne 0 then begin
          wdout[k]=missing
          errout[k]=missing
        endif
;
        data->setvar,wdout,iwin
        data->seterr,errout,iwin
      end
    endfor
 endif 

  IF hkwavecorr THEN data->sethkpixcorr

;
  xkill,wbase
  obj_destroy,cal
;
; SAVE data to fits file
;
  if keyword_set(save) then begin
    file = data->getfilename()
    break_file,file,disk,dir,name,ext
    if not keyword_set(datfile) then $
      datfile = str_replace(name,'l0','l1')+'.fits'
    if not keyword_set(errfile) then $
      errfile = str_replace(name,'l0','er')+'.fits'
    if keyword_set(outdir) then begin
      datfile = concat_dir(outdir,datfile)
      errfile = concat_dir(outdir,errfile)
    endif

    message,' saving data to  '+datfile,/info
    data->save,file=datfile
    message,' saving error to '+errfile,/info
    data->saveerr,file=errfile
    file = datfile
  endif

   ;
 ; PRY, 17-Oct-2019
 ;   Depending if the input FILE was a string (file_type=0) or an
 ;   object (file_type=1), then I do a clean-up of DATA.
 ;
  IF file_type EQ 1 THEN file=temporary(data) ELSE obj_destroy,data

  if n_elements(cleanup) ne 0 then obj_destroy,data
;
  return
end