PRO resample, x, y, x1, x2, dx

; Resample a function y(x) onto a uniform grid
; of stepsize dx between x1 and x2.

colors = GetColor(/Load)


; ** WARNING **
;
; Where the two wavelength scales are very similar, this
; procedure used to introduce significant numerical noise.
; This is a double precision version and seems to work better.

; ** WARNING ** 14.06.2017
;
; This routine is WRONG when resampling into binb sizes much
; larger than original. Compare with 'bin' (this folder) whicj
; gives much better results. 
; 
; TODO -- resolve source of differences and correct both 


  x=x*1.d0
  y=y*1.d0
  nx = n_elements(x)

; define new grid
  nxx = (x2-x1)/dx + 1L
  if (floor(nxx) mod 2 eq 1) then nxx = nxx + 1.0
  xx = Dindgen(nxx)*dx + x1
  yy = DblArr(nxx)

  c1 = 0L
  dx2 = dx/2.0D0
  FOR ix = 0L,nxx-2 DO BEGIN ; loop over new grid

    b1 = xx(ix)-dx2  &  b2 = xx(ix)+dx2
    c2 = 0L 

    IF xx(ix) GT x(0) AND xx(ix) LT x(nx-1) THEN BEGIN

      j1 = c1
      FOR jx = j1,nx-2 DO BEGIN       ; locate data
        xg = x(jx) 
        dg2 = 0.5D0*(x(jx+1)-x(jx))
        dg1 = dg2
        IF jx GT 0 THEN dg1 = 0.5D0*(x(jx)-x(jx-1))
        IF xg+dg2 LE b1 THEN c1 = min([jx+1,nx-2])
	IF xg-dg1 LT b2 THEN c2 = max([c1,jx])
        IF xg-dg1 GE b2 THEN GOTO, next_stage
      ENDFOR

  next_stage:

      IF c2 EQ c1 THEN BEGIN
        s = y(c1)

      ENDIF ELSE BEGIN
 
        IF c2 GT 1.01*c1 THEN BEGIN                     ; bin up 
          s = y(c1)*(0.5*(x(c1)+x(c1+1)) - b1)
          IF c2 GT c1+1 THEN BEGIN
            FOR jx = c1+1,c2-1 DO BEGIN
              s = s + y(jx)*0.5D0*(x(jx+1)-x(jx-1))
            ENDFOR
          ENDIF
          s = s + y(c2)*(b2-0.5D0*(x(c2-1)+x(c2)))
          s = s / (b2-b1) 
	  IF ( NOT FINITE(s) ) THEN PRINT,'WARNING - resample: s=NaN',c2,x(c2),x(c2+1),x(c2+1)-x(c2),b2,b1,b2-b1

  
          ENDIF ELSE BEGIN                          ; interpolate
	  s = ( y(c1) + y(c1+1) ) / 2.0d0
	  IF ( (x(c1+1)-x(c1)) GT 0.d0 ) THEN s = y(c1) + (y(c1+1)-y(c1))*(xx(ix)-x(c1))/(x(c1+1)-x(c1))
	  IF ( NOT FINITE(s) ) THEN PRINT,'WARNING - resample: s=NaN',c1,x(c1),x(c1+1),x(c1+1)-x(c1)

        ENDELSE
      ENDELSE


    ENDIF ELSE BEGIN

      IF xx(ix) LE x(0) THEN s = y(0)
      IF xx(ix) GE x(nx-1) THEN s = y(nx-1)

    ENDELSE

    yy(ix) = s
  ENDFOR

  
  x = xx(0:nxx-2)
  y = yy(0:nxx-2)

END