FUNCTION eis_fix_windata, windata, thresh=thresh, quiet=quiet, night_thresh=night_thresh

;+
; NAME:
;     EIS_FIX_WINDATA
;
; PURPOSE:
;     This routine sets bad data within the windata structure to be
;     missing. The "bad data" includes:
;       - anomalously bright exposures (probably 2048 DN data)
;       - exposures with zero exposure time
;       - exposures during spacecraft night (eclipse season)
;
;     UPDATE: the 2048 data and zero exposure time data should be
;     automatically flagged by eis_prep as of Dec 2019.
;
; CATEGORY:
;     Hinode/EIS; windata.
;
; CALLING SEQUENCE:
;     Result = EIS_FIX_WINDATA( Windata )
;
; INPUTS:
;     Windata: An EIS data structure in the format returned by the
;              routine EIS_GETWINDATA.PRO.
;
; OPTIONAL INPUTS:
;     Thresh:  The threshold for identifying anomalously bright data
;              columns is set to thresh_factor*median value of
;              data. By setting THRESH=, the value of thresh_factor
;              can be modified. The current default is 4.0.
;     Night_Thresh:  This is threshold for detecting low signal data
;              during spacecraft night (eclipse season). By default,
;              the routine sets the threshold to be 0.3*median value
;              data. 
;	
; KEYWORD PARAMETERS:
;     QUIET:   If set, then the routine does not display the data in a
;              graphics window.
;
; OUTPUTS:
;     Returns the input WINDATA structure, but with any "bad" data set
;     to missing.
;
; EXAMPLE:
;     IDL> new_wd=eis_fix_windata(wd)
;
; MODIFICATION HISTORY:
;     Ver.1, 18-May-2011, Peter Young
;     Ver.2, 25-Oct-2011, Peter Young
;        added THRESH= optional input, and /QUIET keyword; changed
;        default threshold value to 4.
;     Ver.3, 14-Aug-2012, Peter Young
;        corrected bug whereby exposures above threshold value were
;        not being set to missing when /quiet was set.
;     Ver.4, 5-Nov-2012, Peter Young
;        now exits if the median value is found to be negative (this
;        happened for a spectral atlas).
;     Ver.5, 29-Jan-2019, Peter Young
;        If the median is negative (see above), now artificially force
;        the data to be positive.
;     Ver.6, 10-May-2020, Peter Young
;        Updated to automatically flag night time exposures; updated
;        header.
;     Ver.7, 18-Dec-2020, Peter Young
;        Added a check on nexp_prp as the routine doesn't work
;        for multiple exposure data.
;-



IF n_params() LT 1 THEN BEGIN
  print,'Use:  IDL> wdout=eis_fix_windata(wd [, /quiet, thresh=, night_thresh= ] )'
  print,'             (Default:   thresh=4)'
  return,windata
ENDIF 

IF windata.hdr.nexp_prp GT 1 THEN BEGIN
   print,'% EIS_FIX_WINDATA: Sorry, this routine does not work for nexp_prp > 1 data. Returning...'
   return,windata
ENDIF 

IF n_elements(thresh) EQ 0 THEN thresh_factor=4.0 ELSE thresh_factor=float(thresh)

wdout=windata

k=where(windata.exposure_time EQ 0.,nk)

IF nk NE 0 THEN BEGIN
  IF NOT keyword_set(quiet) THEN BEGIN
    print,'% EIS_FIX_WINDATA: '+trim(nk)+' exposures have zero exposure time. They have been set to missing.'
  ENDIF 
  wdout.err[*,k,*]=windata.missing
ENDIF 

int=windata.int
k=where(windata.err EQ windata.missing,nk)
IF nk NE 0 THEN int[k]=windata.missing

im_chck=average(int,1,missing=windata.missing)
exposure_av=average(im_chck,2,missing=windata.missing)
k=where(exposure_av NE windata.missing)
val=median(exposure_av[k])

IF n_elements(night_thresh) EQ 0 THEN night_thresh=0.3*val


;
; The method doesn't work if VAL is less than zero, so increase
; the data values to force them to be positive.
;
IF val LE 0. THEN BEGIN
  IF NOT keyword_set(quiet) THEN BEGIN
    print,'%EIS_FIX_WINDATA: median of data was below zero, so data forced to be > 0.'
  ENDIF 
  min_val=min(exposure_av[k])
  exposure_av[k]=exposure_av[k]-min_val
  val=median(exposure_av[k])
ENDIF


IF NOT keyword_set(quiet) THEN BEGIN
  !p.multi=[0,2,1]
  plot_image,im_chck, $
             xtit='Exposure number', $
             ytit='Solar-Y pixels', $
             charsiz=1.5

  plot,exposure_av>0,psym=10,/xsty, $
       xtit='Exposure number', $
     ytit='Average intensity', $
             charsiz=1.5

  IF night_thresh NE windata.missing THEN BEGIN
    oplot,th=2,linestyle=2,[0,1e5],night_thresh*[1,1]
    xyouts,0,night_thresh,'Night threshold',charsiz=1.5
  ENDIF 
  
  oplot,th=2,[0,1e5],val*thresh_factor*[1,1]
  !p.multi=0
ENDIF 



;
; Flag pixels for which the mean exposure value is greater than
; thresh_factor times the median.
;
k=where(exposure_av GE thresh_factor*val,nk)

IF night_thresh NE windata.missing THEN BEGIN
  j=where(exposure_av LE night_thresh,nj)
ENDIF ELSE BEGIN
  nj=0
ENDELSE 

IF nk GT 0 THEN BEGIN 
  IF NOT keyword_set(quiet) THEN BEGIN
    print,'% EIS_FIX_WINDATA: '+trim(nk)+' exposures are anomalously bright. They have been set to missing.'
  ENDIF 
  wdout.int[*,k,*]=windata.missing
ENDIF 

IF nj GT 0 THEN BEGIN 
  IF NOT keyword_set(quiet) THEN BEGIN
    print,'% EIS_FIX_WINDATA: '+trim(nj)+' exposures have been flagged as night-time exposure. They have been set to missing.'
    print,'                   Threshold value: '+trim(string(format='(f8.1)',night_thresh))
    print,'                   (Use NIGHT_THRESH= to manually change this value.)'
  ENDIF 
  wdout.int[*,j,*]=windata.missing
ENDIF 

return,wdout

END