function eis_bin_windata, wd, xbin=xbin, ybin=ybin, offset=offset, $
                          xstart=xstart, ystart=ystart, yshift_wvl=yshift_wvl 


;+
; NAME:
;     EIS_BIN_WINDATA()
;
; PURPOSE:
;     Takes an EIS window data structure created by eis_getwindata, and
;     performs binning in solar-X and/or solar-Y. The error array is
;     correctly calculated, including missing pixels. E.g., if binning
;     is performed over a 3x3 block, and two of the pixels in one
;     block are missing, the error array is computed from the
;     remaining seven pixels. 
;
;     Note that, before averaging is done in the intensity array,
;     missing pixels are flagged as missing in the intensity array
;     (eis_getwindata replaces missing pixels in the intensity array
;     with a median of neighbouring pixels, which would mess up the
;     averaging process).
;
;     It is generally recommended that either XBIN is not used, or set to
;     small values (say 2 or 3). This is especially true for data-sets that
;     have missing exposures.
;
;     Since September 2010 the windata structure has contained the tags
;     WAVE_CORR, WAVE_CORR_TILT and WAVE_CORR_T which specify the
;     spectrum drift and slit tilt. These are automatically re-binned
;     along with the the other arrays.
;
; CATEGORY:
;     Hinode/EIS; windata
;
; CALLING SEQUENCE:
;     Result = EIS_BIN_WINDATA( Windata )
;
; INPUTS:
;     Wd:   An EIS window data structure created by eis_getwindata.
;
; OPTIONAL INPUTS:
;     Xbin: An integer specifying the binning in X. E.g., xbin=3 implies 3
;           spatial pixels in X are averaged. The intensity array will then be
;           reduced in size by a factor 3 in the X direction following the
;           averaging.
;
;     Ybin: An integer specifying the binning in Y. E.g., ybin=3 implies 3
;           spatial pixels in Y are averaged. The intensity array will then be
;           reduced in size by a factor 3 in the Y direction following the
;           averaging.
;
;     Xstart: By default the arrays will be re-binned starting from
;            X=0. The starting pixel can be changed by specifying
;            XSTART. E.g., if XBIN=10 then the default behavior is to
;            bin pixels 0:9, 10:19, etc. By specifying XSTART=5, the
;            binning would be 5:14, 15:24, etc.
;
;     Ystart:By default the arrays will be re-binned starting from
;            Y=0. The starting pixel can be changed by specifying
;            YSTART. E.g., if YBIN=10 then the default behavior is to
;            bin pixels 0:9, 10:19, etc. By specifying YSTART=5, the
;            binning would be 5:14, 15:24, etc.
;
;     Yshift_Wvl:  This is useful if you want to match the macropixels
;             with another wavelength window. By specifying the
;             wavelength of the reference window, the YSTART value
;             will be automatically computed to line up the
;             macropixels between the two windows. Note that this
;             overwrites the YSTART value (if set).
;
;     Offset: This keyword is now obsolete. The OFFSET array is now stored
;             in WINDATA.
;
; OPTIONAL INPUT/OUTPUT:
;
;     OFFSET  A 2-dimensional array of size (WD.NX,WD.NY) that
;             contains the combined spectrum drift/slit tilt
;             array. This is re-binned and returned to the user, *so
;             the input array is modified*. Note that OFFSET= is
;             intended for cases where the user is not using the
;             WAVE_CORR of WD in their analysis (see the documentation
;             about the EIS_AUTO_FIT suite of routines for more
;             details).
;
; OUTPUTS:
;     A structure with the same tags as WD is returned, but the
;     intensity and error arrays will be smaller following the
;     binning. If XBIN has been specified, then the .TIME,
;     .EXPOSURE_TIME and .TIME_CCSDS tags will also be modified. The
;     WAVE_CORR, WAVE_CORR_T and WAVE_CORR_TILT tags are also affected
;     by the re-binning.
;
;     In addition a *new* tag called START_PIX is added to the
;     structure. This is a two element array containing
;     [XSTART,YSTART]. This tag is non-standard and care should be
;     taken in using it in code.
;
; EXAMPLE:
;     wd=eis_getwindata(filename,195.12)
;     wd2=eis_bin_windata(wd,xbin=3,ybin=3,ystart=15)
;     wd2=eis_bin_windata(wd,xbin=3,ybin=3,yshift_wvl=284.16)
;
; MODIFICATION HISTORY:
;     Ver.1, 5-Nov-2008, Peter Young
;     Ver.2, 6-Nov-2008, Peter Young
;       modified SCALE tag based on binning factors
;     Ver.3, 11-Dec-2009, Peter Young
;       added OFFSET, XSTART and YSTART optional inputs; fixed bug
;       which occurs when all pixels are missing in a binned pixel
;       (spotted by Tongjiang Wang).
;     Ver.4, 18-Dec-2009, Peter Young
;       I have introduced a new tag called START_PIX to contain
;       [xstart,ystart]. This is non-standard so care has to be taken
;       in using this tag.
;     Ver.5, 11-Apr-2011, Peter Young
;       The WAVE_CORR, WAVE_CORR_TILT and WAVE_CORR_T tags are now
;       re-binned; OFFSET will replace WD.WAVE_CORR if it is
;       specified.
;     Ver.6, 11-May-2011, Peter Young
;       Now re-computes XCEN and YCEN.
;     Ver.7, 26-Mar-2014, Peter Young
;       Added check on which instrument the data are from.
;     Ver.8, 19-May-2020, Peter Young
;       Added YSHIFT_WVL input; updated header.
;     Ver.9, 07-Oct-2021, Peter Young
;       Slot data doesn't have the wave_corr tags so I've modified the
;       code to check for wave_corr.
;-


IF n_params() LT 1 THEN BEGIN
  print,'Use:  IDL> wdout=eis_bin_windata( wd [, xbin=, ybin=, xstart=, ystart=, yshift_wvl= ])'
  return,-1
ENDIF 

;
; 'instr' takes the values
;    0 - Hinode/EIS
;    1 - SOHO/CDS
;    2 - IRIS
;
instr=0
IF strlowcase(trim(wd.hdr.telescop)) EQ 'iris' THEN instr=2

;
; Check if we have EIS slot data.
;
IF wd.hdr.slit_ind EQ 1 OR wd.hdr.slit_ind EQ 3 THEN slot=1 ELSE slot=0

;
; Check if the wave_corr tags exist.
;
wc_swtch=tag_exist(wd,'wave_corr')


IF n_elements(xbin) EQ 0 AND n_elements(ybin) EQ 0 THEN BEGIN
  print,'% EIS_BIN_WINDATA: please specify at least one of XBIN or YBIN. Returning...'
  return,-1
ENDIF

wd2=wd

;
; This only applies to EIS narrow slit data.
; 
IF wc_swtch EQ 0 AND instr EQ 0 AND slot EQ 0 THEN BEGIN
  print,'% EIS_BIN_WINDATA: Your WINDATA structure is out-of-date and does not contain the WAVE_CORR tag.'
  print,'                   Please re-generate the structure with EIS_GETWINDATA before calling EIS_BIN_WINDATA.'
  print,'                   Returning...'
  return,-1
ENDIF 

IF n_elements(offset) NE 0 THEN BEGIN
  siz=size(offset)
  IF siz[1] NE wd.nx OR siz[2] NE wd.ny THEN BEGIN
    print,'% EIS_BIN_WINDATA: the input OFFSET= has been specified, but the dimensions do not match those of WINDATA.'
    print,'                   Please check your data.  Returning...'
    return,-1
  ENDIF 
ENDIF 


;
; I have to deal with bad data here that weren't dealt with by eis_prep or
; eis_getwindata. I look for any exposure times that are zero and then set
; the bad data to missing where applicable.
;
k=where(wd.exposure_time EQ 0.,nk)
IF nk GT 0 THEN BEGIN
  wd.solar_x[k]=-1e6
  IF wc_swtch THEN BEGIN 
    wd.wave_corr[k,*]=wd.missing
    wd.wave_corr_t[k]=wd.missing
  ENDIF 
ENDIF 



IF n_elements(xstart) NE 0 THEN BEGIN
  IF xstart GE wd.nx-1 THEN BEGIN
    print,'% EIS_BIN_WINDATA: XSTART must be smaller than size of array. Returning...'
    return,-1
  ENDIF
ENDIF ELSE BEGIN
  xstart=0
ENDELSE 

IF n_elements(ystart) NE 0 THEN BEGIN
  IF ystart GE wd.ny-1 THEN BEGIN
    print,'% EIS_BIN_WINDATA: YSTART must be smaller than size of array. Returning...'
    return,-1
  ENDIF
ENDIF ELSE BEGIN
  ystart=0
ENDELSE

;
; This overwrites the ystart value (if set).
;
IF n_elements(yshift_wvl) NE 0 THEN BEGIN
  y1=eis_ccd_offset(mean(wd.wvl))
  y2=eis_ccd_offset(yshift_wvl)
 ;
  ypix=round(y1[0]-y2[0])
  IF ypix GE 0 THEN BEGIN
    ystart=ypix
  ENDIF ELSE BEGIN
    a=-ypix/float(ybin)
    n=ceil(a)
    ystart=ypix+n*ybin
  ENDELSE
  print,'% EIS_BIN_WINDATA: yshift_wvl value gives ystart='+trim(ystart)+'.'
ENDIF 


;
; I need to create new XCEN and YCEN values. First I find the origin
; of the original array, and then compute the new origin.
; Note that I'm using a simplistic method of computing the
; origin that assumes an X-pixel has a width SCALE[0]. This is
; consistent, however, with the usage of the PLOT_IMAGE routine. 
;
scale=wd.scale
IF wd.hdr.slit_ind EQ 2 THEN slit_wid=2.0 ELSE slit_wid=1.0
orig=fltarr(2)
orig[0]=wd.xcen - float(wd.nx)/2.*scale[0]
orig[1]=wd.ycen - float(wd.ny)/2.*scale[1]
orig_new=orig+[xstart,ystart]*wd.scale


nl=wd.nl
nx=wd.nx-xstart
ny=wd.ny-ystart
time=wd.time[xstart:*]
time_ccsds=wd.time_ccsds[xstart:*]
exposure_time=wd.exposure_time[xstart:*]
solar_x=wd.solar_x[xstart:*]
solar_y=wd.solar_y[ystart:*]
data_quality=wd.data_quality[xstart:*]
scale=wd.scale
IF wc_swtch THEN BEGIN 
  wave_corr=wd.wave_corr[xstart:*,ystart:*]
  wave_corr_tilt=wd.wave_corr_tilt[ystart:*]
  wave_corr_t=wd.wave_corr_t[xstart:*]
ENDIF 

;
; These are the 'start' arrays. Specifying XSTART or YSTART will
; modify them.
;
int0=wd.int[*,xstart:*,ystart:*]
err0=wd.err[*,xstart:*,ystart:*]

;
; If OFFSET has been specified, then re-define
;
IF n_elements(offset) NE 0 THEN offset=offset[xstart:*,ystart:*]

;
; Compute new array sizes based on XBIN and YBIN inputs.
; SWTCH indicates what binning is being performed:
;   0   - no binning
;   1   - binning only in X
;   10  - binning only in Y
;   11  - binning in X and Y
;
swtch=0
IF n_elements(xbin) NE 0 THEN BEGIN
  mx=nx/xbin
  swtch=swtch+1
  scale[0]=scale[0]*xbin
ENDIF ELSE BEGIN
  mx=nx
  xbin=1
ENDELSE
;
IF n_elements(ybin) NE 0 THEN BEGIN
  my=ny/ybin 
  swtch=swtch+10
  scale[1]=scale[1]*ybin
ENDIF ELSE BEGIN
  my=ny
  ybin=1
endelse

;
; INDEX is used to pick out the index of the last pixel in each group
; of binned X-pixels. This is used, e.g., for extracting the exposure
; start times.
;
IF n_elements(xbin) NE 0 THEN index=indgen(mx)*xbin+xbin-1 ELSE index=indgen(nx)


;
; Create a new windata structure with the correct dimensions. This
; involves removing the old tags, and then adding them again with the
; new dimensions.
;
wd3=rem_tag(wd2,['int','err','time','time_ccsds','exposure_time', $
                 'solar_x','solar_y','data_quality','wave_corr', $
                 'wave_corr_t','wave_corr_tilt'])
wd2=add_tag(wd3,fltarr(nl,mx,my),'int')
wd2=add_tag(wd2,fltarr(nl,mx,my),'err')
wd2=add_tag(wd2,dblarr(mx),'time')
wd2=add_tag(wd2,strarr(mx),'time_ccsds')
wd2=add_tag(wd2,fltarr(mx),'exposure_time')
wd2=add_tag(wd2,fltarr(mx),'solar_x')
wd2=add_tag(wd2,fltarr(my),'solar_y')
wd2=add_tag(wd2,fltarr(mx),'data_quality')
IF wc_swtch THEN BEGIN 
  wd2=add_tag(wd2,fltarr(mx,my),'wave_corr')
  wd2=add_tag(wd2,fltarr(my),'wave_corr_tilt')
  wd2=add_tag(wd2,fltarr(mx),'wave_corr_t')
ENDIF 
wd2=add_tag(wd2,intarr(2),'start_pix')   ; new tag
wd2.nx=mx
wd2.ny=my
wd2.scale=scale
;
wd3=0

wd2.start_pix=[xstart,ystart]

;
; Need to re-define time arrays if X-binning
;  - new exposure time is XBIN*(old exposure time)
;  - I'm not sure what to do with data_quality (!)
;
IF n_elements(xbin) NE 0 THEN BEGIN
  wd2.time=time[index]
  wd2.time_ccsds=time_ccsds[index]
  FOR i=0,mx-1 DO wd2.exposure_time[i]=total(exposure_time[i*xbin:(i+1)*xbin-1])
  wd2.data_quality=data_quality[index]
  FOR i=0,mx-1 DO wd2.solar_x[i]=average(solar_x[i*xbin:(i+1)*xbin-1],missing=-1e6)
ENDIF ELSE BEGIN
  wd2.time=time
  wd2.time_ccsds=time_ccsds
  wd2.exposure_time=exposure_time
  wd2.data_quality=data_quality
  wd2.solar_x=solar_x
ENDELSE

;
; Work out new XCEN and YCEN
;
wd2.xcen=orig_new[0] + float(wd2.nx)/2.*wd2.scale[0]
wd2.ycen=orig_new[1] + float(wd2.ny)/2.*wd2.scale[1]


;
; Only solar_y tag needs re-defining if Y-binning
;  - when averaging solar_y there should be a 'missing' value,
;    but I don't know what Viggo assigns for missing pointing values.
;
IF n_elements(ybin) NE 0 THEN BEGIN
  FOR i=0,my-1 DO wd2.solar_y[i]=average(solar_y[i*ybin:(i+1)*ybin-1])
ENDIF ELSE BEGIN
  wd2.solar_y=solar_y
ENDELSE 


CASE swtch OF 

 ; --- Binning in X only ---
  1: BEGIN
    FOR i=0,mx-1 DO BEGIN
      k0=i*xbin
      k1=(i+1)*xbin-1
     ;
      int=int0[*,k0:k1,*]
      err=err0[*,k0:k1,*]
     ;
      imiss=where(err EQ wd.missing,nmiss)
      IF nmiss NE 0 THEN int[imiss]=wd.missing
     ;
      wd2.int[*,i,*]=average(int,2,missing=wd.missing)
     ;
      err2=err^2
      IF nmiss NE 0 THEN err2[imiss]=wd.missing
     ;
      av_err2=average(err2,2,missing=wd.missing)
      i_notmiss=where(av_err2 NE wd.missing,n_notmiss)
     ;
      IF n_notmiss NE 0 THEN av_err2[i_notmiss]=sqrt(av_err2[i_notmiss])/sqrt(xbin)
      wd2.err[*,i,*]=av_err2
    ENDFOR
  END

 ; --- Binning in Y only ---
  10: BEGIN
    FOR i=0,my-1 DO BEGIN
      k0=i*ybin
      k1=(i+1)*ybin-1
     ;
      int=int0[*,*,k0:k1]
      err=err0[*,*,k0:k1]
     ;
      imiss=where(err EQ wd.missing,nmiss)
      IF nmiss NE 0 THEN int[imiss]=wd.missing
     ;
      wd2.int[*,*,i]=average(int,3,missing=wd.missing)
     ;
      err2=err^2
      IF nmiss NE 0 THEN err2[imiss]=wd.missing
     ;
      av_err2=average(err2,3,missing=wd.missing)
      i_notmiss=where(av_err2 NE wd.missing,n_notmiss)
     ;
      IF n_notmiss NE 0 THEN av_err2[i_notmiss]=sqrt(av_err2[i_notmiss])/sqrt(ybin)
      wd2.err[*,*,i]=av_err2
    ENDFOR
  END

 ; --- Binning in X and Y ---
  11: BEGIN
    FOR i=0,mx-1 DO BEGIN
      FOR j=0,my-1 DO BEGIN
        k0=i*xbin
        k1=(i+1)*xbin-1
       ;
        l0=j*ybin
        l1=(j+1)*ybin-1
       ;
        int=int0[*,k0:k1,l0:l1]
        err=err0[*,k0:k1,l0:l1]
       ;
        imiss=where(err EQ wd.missing,nmiss)
        IF nmiss NE 0 THEN int[imiss]=wd.missing
       ;
        wd2.int[*,i,j]=average(int,[2,3],missing=wd.missing)
       ;
        err2=err^2
        IF nmiss NE 0 THEN err2[imiss]=wd.missing
       ;
        av_err2=average(err2,[2,3],missing=wd.missing)
        i_notmiss=where(av_err2 NE wd.missing,n_notmiss)
       ;
        IF n_notmiss NE 0 THEN av_err2[i_notmiss]=sqrt(av_err2[i_notmiss])/sqrt(xbin*ybin)
        wd2.err[*,i,j]=av_err2
      ENDFOR
    ENDFOR
  END

ENDCASE


;
; Update the three WAVE_CORR arrays in WINDATA
;
IF wc_swtch THEN BEGIN 
  wave_corr=wd.wave_corr
  wave_corr_tilt=wd.wave_corr_tilt
  wave_corr_t=wd.wave_corr_t
;
  FOR i=0,mx-1 DO BEGIN
    FOR j=0,my-1 DO BEGIN
      k0=i*xbin
      k1=(i+1)*xbin-1
   ;
      l0=j*ybin
      l1=(j+1)*ybin-1
   ;
      wd2.wave_corr[i,j]=average(wave_corr[k0:k1,l0:l1],missing=-100.)
    ENDFOR
  ENDFOR
;
  FOR j=0,my-1 DO BEGIN
    l0=j*ybin
    l1=(j+1)*ybin-1
 ;
    wd2.wave_corr_tilt[j]=average(wave_corr_tilt[l0:l1])
  ENDFOR
;
  new_wave_corr_t=fltarr(mx)
  FOR i=0,mx-1 DO BEGIN
    k0=i*xbin
    k1=(i+1)*xbin-1
 ;
    wd2.wave_corr_t[i]=average(wave_corr_t[k0:k1],missing=-100.)
  ENDFOR
ENDIF 

;
; If OFFSET been specified, then overwrite the WAVE_CORR tag with the
; re-binned OFFSET array.
;
; NOTE: ideally I should also re-calculate the WAVE_CORR_TILT and
; WAVE_CORR_T tags, too. These aren't important for eis_auto_fit so
; I'm going to ignore them for now.
;
IF n_elements(offset) ne 0 THEN BEGIN 
  new_offset=fltarr(mx,my)    ; re-binned array
  FOR i=0,mx-1 DO BEGIN
    FOR j=0,my-1 DO BEGIN
      k0=i*xbin
      k1=(i+1)*xbin-1
     ;
      l0=j*ybin
      l1=(j+1)*ybin-1
     ;
      new_offset[i,j]=average(offset[k0:k1,l0:l1])
    ENDFOR
  ENDFOR
  offset=new_offset
ENDIF 

return,wd2

END