FUNCTION PROSHIFT,cube,ds,reb
;
;Calculation of line profile shift by means of a two-slit method.
;Two "slits" with a constant distance DS between them are shifted
;along the line profile(s) (last dimension of CUBE) in order to
;achieve equal intensities on both slits. Then the position of the
;"slits" shows the relative shift(s) of the line profile(s).
;
;INPUTS: CUBE - 1D-3D array, where the spectral dimension is last
;	    DS - distance between the slits (pixels)
;	    REB - rebinning coefficient to make the spectral sampling finer
;OUTPUT: relative shift of the profile in pixels (float)

on_error,1
si=size(cube)

if si(0) eq 1 then begin	 ;a single spectral profile (lambda)
  cub=rebin(cube,si(1)*reb)	 ;rebin to refine spectral resolution
  cub=cub(0:(si(1)-1)*reb)	 ;cut away the constant tail
  difc=cub-shift(cub,ds*reb) 	 ;compare shifted spectra
  cub=0
  difc=abs(difc(ds*reb:*))	 ;abs difference in the useful range
  mi=min(difc,pos)		 ;at which position is the minimum of difference?
  RETURN,float(pos)/reb		 ; a scalar
endif

if si(0) eq 2 then begin		;1D spectral image (x,lambda)
  cub=rebin(cube,si(1),si(2)*reb)	;rebin to refine spectral resolution
  cub=cub(*,0:(si(2)-1)*reb)	 	;cut away the constant tail
  difc=cub-shift(cub,0,ds*reb)		;compare shifted spectra
  cub=0
  difc=abs(difc(*,ds*reb:*))		;abs difference in the useful range
  prosh=fltarr(si(1))
				;at which position is the minimum of difference?
  for i=0,si(1)-1 do begin
    mi=min(difc(i,*),pos)
    prosh(i)=float(pos)
  endfor
  RETURN,prosh/reb		;a 1D float vector
endif

if si(0) eq 3 then begin			;2D spectral image (x,y,lambda)
  cub=rebin(cube,si(1), si(2),si(3)*reb) 	;rebin to refine spectral resolution
  cub=cub(*,*,0:(si(3)-1)*reb)	 		;cut away the constant tail
  difc=cub-shift(cub,0,0,ds*reb)		;compare shifted spectra
  cub=0
  difc=abs(difc(*,*,ds*reb:*))			;abs difference in the useful range
  prosh=fltarr(si(1),si(2))
				;at which position is the minimum of difference?
  for i=0,si(1)-1 do begin
  for j=0,si(2)-1 do begin
    mi=min(difc(i,j,*),pos)
    prosh(i,j)=float(pos)
  endfor
  endfor
  RETURN,prosh/reb		;a 2D float array
endif

message,'Input array must have 1-3 dimensions'
END