document.write(pagetitle)

Home Download Tutorial 1 Tutorial 2 Tutorial 3 Tutorial 4 Tutorial 5
```;+
; NAME:              XREBIN
;
; PURPOSE:           Rebin a vector
;
; USAGE:
;         yout = xrebin(xin, yin, xout [,/SPLINE] [,/LSQUADRATIC]
; ARGUMENTS:
;         XIN         (N+1) points of the borders of each input bin
;         YIN         values (N) of input bins (1 or 2D)
;         XOUT        points of the borders of output bins
; KEYWORDS:
;         SPLINE      Spline interpolation (default: linear)
;         SPLINF      Other version of spline interpolation
;                     3-4 times faster than SPLINE, the results are
;                     slightly different, but seemingly better.
;         CUBIC       Use the function INTERPOLATE with CUBIC=-05
;                     Cubic convolution is an interpolation method that
;                     closely approximates the theoretically optimum sinc
;                     interpolation function using cubic polynomials.
;         NAN         To treat the missing values in YIN as
;                     zeros when interpolating it. Use with caution.
; OUTPUT
;         YOUT        output 1 or 2D array, rebined at the borders
;                     requested by XOUT.
;                     YOUT is in double precision (it shall be converted
;                     in float if appropriate.
;
; DESCRIPTION:
;    This function rebins a vector or a n-dimensional array into
;    different bins along the first axis.
;
;    The rebining is made as integrating the signal, interpolating and
;    differentiating the result. The integrated signal is preserved.
;
;    Different interpolating algorithms are proposed, and only one of
;    If not of them is specified, the routine INTERPOL will make a linear
;    interpolation. This latter routin also handle the SPLINE LSQUADRATIC
;    QUADRATIC cases. The difference between SPLINF and SPLINE is
;    that the spline representation is computed on the whole vector
;    instead of on 4-points blocks (this is faster and often our preferred
;    option). CUBIC is using the function INTERPOLATE with CUBIC=-0.5 (sinc).
;
;    No check is performed to prevent extrapolation (beware to give proper input).
;
;    The computations are systematically made in double precision. The
;    reason is that the algorithm results in loss of relative precision equal
;    to the number of pixels. In float, this is often a unacceptable.
;
; AUTHOR: Mina Koleva (2007/03/01)
;-
; CATEGORY:    ULY_UTIL
;------------------------------------------------------------------------------

compile_opt idl2
on_error, 2

s = size(yin)
if s[0] lt 1  then message,'Yin should be a vector or array'

if s[1]+1 ne (size(xin))[1] then message,'Xin must have one element more than the 1st dim of Yin '

if s[0] eq 1 then begin    ; case of a 1D array

integr = [0, total(double(yin), 1, /CUMULATIVE, NaN=NaN)]  ; integrate yin
;   Note that we force the computation in double precision
;   The integration-differentiation process results in a loss of precision
;   which is often critical in float (and acceptable in double)
;   The loss of precision is equal to the number of pixels, so for 1e6
;   pixels, if the computation were made in float, the resulting precision
;   would not be better than a few 10%...

;interpolate yout for xout
case (1) of
keyword_set(splinf): begin
;           if spl_init generates an underflow error we clear it
y2 = spl_init(xin, integr)
integr_interp = spl_interp(xin, integr, y2, xout)
;           If any of the 4 input vectors is double precision, the
;           result is also in double precision
if status eq 0 then status = check_math(MASK=32)
end

keyword_set(cubic): begin
iout = interpol(lindgen(n_elements(xin)), xin, xout)
integr_interp = interpolate(integr, iout, /GRID, CUBIC=-0.5)
end

else : \$
integr_interp = \$
interpol(integr, xin, xout, \$

endcase
;differentiate
return, (shift(integr_interp,-1)-integr_interp)[0:n_elements(integr_interp)-2]

endif else begin   ; case of a nD array
dim = size(yin, /DIM)
integr = [dblarr([1,dim[1:*]]), total(double(yin), 1, /CUMULATIVE, NAN=NaN)]  ; integrate yin
integr = reform(integr, dim[0]+1, n_elements(integr)/(dim[0]+1), /OVER)
sr = size(integr)

case (1) of
keyword_set(splinf): begin
integr_interp = make_array(n_elements(xout), sr[2], /DOUB, /NOZERO)
for i=0, sr[2]-1 do begin
y2 = spl_init(xin, integr[*,i])
integr_interp[*,i] = spl_interp(xin, integr[*,i], y2, xout)
endfor
if status eq 0 then status = check_math(MASK=32)
end

keyword_set(cubic): begin
iout = interpol(lindgen(n_elements(xin)), xin, xout)
yout = lindgen(sr[2])
integr_interp = interpolate(integr, iout, yout, /GR, CUB=-0.5)
end

else: begin
; note: we can certainly improve this significantly if we avoid the
;  VALUE_LOCATE in interpol
integr_interp = make_array(n_elements(xout), sr[2], /DOUB, /NOZERO)
for i=0, sr[2]-1 do \$
integr_interp[*,i] = \$
interpol(integr[*,i], xin, xout, \$
end
endcase
;differentiate
s1=(size(integr_interp))[1]
return, reform((shift(integr_interp,-1)-integr_interp)[0:s1-2, *], \$
[s1-1, dim[1:*]])
endelse

; note: We can also use the function INTERPOLATE (with /GRID,CUB=-0.5)
;  We would not need the for loop (in the 2D case)
;

end

;------- TESTING XREBIN----------
pro test_xrebin

;test 1
print, 'TEST 1...',FORMAT="(A,\$)"
xin = [1,2,3,4]
yin = [1,2,1]
xout= [1,2,3,4]
yout = xrebin(xin,yin,xout)
a = where(yout - [1,2,1],n)
if n eq 0 then \$
print, 'OK' else \$
print, 'FAILED'

;test 2
print, 'TEST 2...',FORMAT="(A,\$)"
xin = [1,2,3,4]
yin = [1.,2.,1.]
xout = [1,1.5,2,2.5,3,3.5,4]
yout = xrebin(xin,yin,xout)
a = where(yout - [.5,.5,1,1,.5,.5],n)
if n eq 0 then \$
print, 'OK' else \$
print, 'FAILED'

;test 3
print, 'TEST 3...',FORMAT="(A,\$)"
xin = [1,1.5,2,2.5,3,3.5,4]
yin = [0.5,0.5,1.,1.,.5,.5]
xout = [1,2,3,4]
yout = xrebin(xin,yin,xout)
a = where(yout - [1,2,1],n)
if n eq 0 then \$
print, 'OK' else \$
print, 'FAILED'

end

;--- end ---

```
Part of ULySS package