PRO EIS_FIT_Y_SHIFT, INT, ERR, OFFSET, OUTINT, OUTERR, MISSING=MISSING

;+
; NAME
;
;     EIS_FIT_Y_SHIFT
;
; PROJECT
;
;     Hinode/EIS
;
; EXPLANATION
;
;     Takes 2D intensity and error arrays and works out new arrays
;     assuming a fractional spatial offset.
;
;     If a bin has a missing value, then each of the re-computed bins
;     that use this value will be set to missing. Note that missing
;     data are flagged through the missing values in ERR, not in INT. 
;
;     The offset must be sub-pixel.
;
;     This routine was used to perform the sub-pixel shifts to Fe XII
;     and Fe XIII density diagnostic lines in Young et al. (2009, A&A,
;     495, 587) - see Appendix C of this work.
;
; INPUTS
;
;     INT      A 2D array assumed to contain emission line intensities
;              over a spatial area.
;
;     ERR      A 2D array assumed to contain errors on emission line
;              intensities over a spatial area. Missing data points
;              should be flagged with the value MISSING (see below).
;
;     OFFSET   A positive value indicates that the intensity image is
;              moved up on the detector. A negative value indicates that
;              the intensity image is moved down on the detector.
;
; OPTIONAL INPUTS
;
;     MISSING  The value for missing data in the ERR array. If not
;              set, then a value of -100 is assumed.
;
; OUTPUTS
;
;     OUTINT   The output intensity array with corrected intensity
;              values. Missing data are set to the value MISSING.
;
;     OUTERR   The output error array with corrected error
;              values. Missing data are set to the value MISSING.
;
; HISTORY
;
;     Ver.1, 16-Dec-2009, Peter Young
;       - renamed from eis_fit_offset
;     Ver.2, 24-Feb-2014, Peter Young
;       - updated header; no change to code.
;-

siz=size(int)

nx=siz[1]
ny=siz[2]

IF n_elements(missing) EQ 0 THEN missing=-100.

outint=int-int+missing
outerr=err-err+missing

IF abs(offset) GT 1. THEN BEGIN
  print,'%EIS_FIT_Y_SHIFT:  OFFSET must have a value ABS(offset) LE 1. Returning...'
  return
ENDIF 

IF offset LT 0. THEN BEGIN
  alpha=abs(offset)
  FOR i=0,nx-1 DO BEGIN
    y=reform(int[i,0:ny-2])
    y_plus=reform(int[i,1:*])
   ;
    e=reform(err[i,0:ny-2])
    e_plus=reform(err[i,1:*])
   ;
    i_miss=where(e EQ missing OR e_plus EQ missing,n_miss)
   ;
    new_y=(1-alpha)*y + alpha*y_plus
    IF n_miss NE 0 THEN new_y[i_miss]=missing
   ;
    new_e=sqrt( (1.-alpha)^2*e^2 + (alpha)^2*e_plus^2 )
    IF n_miss NE 0 THEN new_e[i_miss]=missing
   ;
    outint[i,0:ny-2]=new_y
    outerr[i,0:ny-2]=new_e
  ENDFOR
ENDIF ELSE BEGIN
  alpha=abs(offset)
  FOR i=0,nx-1 DO BEGIN
    y=reform(int[i,1:*])
    y_minus=reform(int[i,0:ny-2])
   ;
    e=reform(err[i,1:*])
    e_minus=reform(err[i,0:ny-2])
   ;
    i_miss=where(e EQ missing OR e_minus EQ missing,n_miss)
   ;
    new_y=(1-alpha)*y + alpha*y_minus
    IF n_miss NE 0 THEN new_y[i_miss]=missing
   ;
    new_e=sqrt( (1.-alpha)^2*e^2 + (alpha)^2*e_minus^2 )
    IF n_miss NE 0 THEN new_e[i_miss]=missing
   ;
    outint[i,1:*]=new_y
    outerr[i,1:*]=new_e
  ENDFOR
ENDELSE

END