function rebinw,f,x,y,perbin=perbin,nbin=nbin,wrange=wrange,slowok=slowok,\$ numbins=numbins,verbose=verbose, _extra=e ;+ ;function rebinw ; rebins the input array on a new grid and returns the output array. ; this routine is preferred over simple INTERPOL or SPLINE because ; this is guaranteed to e.g., conserve flux while rebinning spectra. ; ;syntax ; ff=rebinw(f,x,y,/perbin,numbins=numbins,nbin=nbin,wrange=wrange,\$ ; /slowok,verbose=verbose) ; ;parameters ; f [INPUT; required] array values ; x [INPUT; required] absissae for F ; 1: if size matches that of F, assumed to be mid-bin values ; 2: if size is size(F)+1, assumed to be bin boundaries ; (i.e., all the bin-beginners and a final bin-ender) ; 3: if neither of the above, this is assumed to be the ; desired *output* grid, and the absissae of the input ; are assumed to span the linear range 1..N(F) ; * Y is ignored on input, but will be overwritten on output ; y [I/O] the new absissa values, taken from NUMBINS if it is set. ; 1: if scalar or 1-element vector on input, assumed to be ; number of bins in output ; * linear binning if +ve, log binning if -ve ; 2: if >1-element vector on input, assumed to be the desired ; output grid of bin boundaries, with the last element ; defining the final bin-ending value. ; 3: if not defined, takes value from NBIN ; NOTE: F extrapolates as zeros. ; ;keywords ; perbin [INPUT] if set, assumes that units of F(X) are [.../bin] ; nbin [I/O] number of bins in output ; * overwritten if defined via X or Y ; * linear binning if +ve, log binning if -ve ; * default: 1 ; numbins [INPUT] same as NBINS, except that if set, it ; overrides NBIN _and_ Y ; * so, just to reiterate: ; if NUMBINS is _not_ defined, then ; -- if Y is a scalar, say = NN, then ; output will be rebinned to have NN elements, and ; Y will be overwritten to contain the new grid, and ; NBIN will be reset to NN ; -- if Y is a vector, then ; output will be rebinned using Y as bin boundaries, and ; NBIN will be reset to N(Y)-1 ; -- if Y is not defined, but NBIN is defined, then ; output will be rebinned to have NBIN elements, and ; Y will on output contain the new grid ; if NUMBINS is defined, then ; -- output will be rebinned to have NUMBINS elements, and ; Y will contain the new grid on output, and ; NBIN will be reset to NUMBINS ; wrange [INPUT] output grid range ; * ignored if defined via array Y (or X -- case 3) ; * overrides if determined via X and CASE 1 of Y ; * default: [1.,500.] ; slowok [INPUT] if set, skips call to FINDEX and does it the ; slow IDL way ; verbose [INPUT] controls chatter ; _extra [JUNK] here only to prevent crashing the program ; ;restrictions ; beware that if Y is not sorted in ascending order, on output ; it will be! ; requires external procedure FINDEX.PRO ; ;history ; vinay kashyap (Jan98) ; forced x to be required inputs (Jun98) ; changed call to INTERPOL to combo call to INTERPOLATE/FINDEX (VK;99Jul) ; corrected long-standing "feature" of NBIN being useless, changed ; default behavior of Y=UNDEFINED, NBIN, and X=/=F (VK; MarMM) ; added keyword SLOWOK (VK; Jul01) ; added keyword NUMBINS (VK; Dec'02) ; added keyword VERBOSE; now uses /CUMULATIVE in TOTAL() (VK; Feb'03) ; beware of /CUMULATIVE in TOTAL(), which returns an array of same ; size as input, and does _not_ begin with 0 (VK; Apr'03) ; added warning if Y gets resorted (VK; Mar04) ; cleaned up behavior in the case of non-unique/non-monotonic inputs ; (VK; Mar05) ;- forward_function findex ; usage ok='ok' & np=n_params() & nf=n_elements(f) & nx=n_elements(x) ny=n_elements(y) if np lt 2 then ok='Insufficient parameters' else \$ if nf eq 0 then ok='Input function undefined' else \$ if nx eq 0 then ok='bin boundaries undefined' else \$ if nf eq 1 then ok='Cannot make a spectrum out of one element' if ok ne 'ok' then begin print,'Usage: g=rebinw(f,x,y,/perbin,numbins=numbins,nbin=nbin,\$' print,' wrange=wrange,/slowok,verbose=verbose)' print,' rebins F(X) into G(Y) on a different grid while conserving norm' if np ne 0 then message,ok,/info return,-1L endif ; verbosity vv=0L & if keyword_set(verbose) then vv=long(verbose[0])>1 ; initialize via keywords forcekey=0 if keyword_set(numbins) then begin forcekey=1 & nbin=long(numbins(0)) & ny=0 if vv gt 1 then message,/informational,\$ 'Constructing output grid from NUMBINS' endif if n_elements(nbin) eq 0 then nbin=1L else nbin=long(nbin(0)) if nbin eq 0 then nbin=1L ;NBIN cannot be 0 if n_elements(wrange) ge 2 then wr=[wrange(0),wrange(1)] else wr=[1.,500.] ; compatibility checks xx=x & wrr=[min(xx),max(xx)] ;case 2: default if nx eq nf then begin ;(case 1: mid-bin values if vv gt 2 then message,/informational,\$ 'Converting input from mid-bin values to bin boundaries' dx=x(1:*)-x & dx=0.5*dx & xx=[x-dx,x(nx-1)+dx(nx-2)*[-1,1]] wrr=[xx(0),xx(nx)] endif else if nx ne (nf+1) then begin ;)(case 3: X is the output grid if vv gt 2 then message,/informational,\$ 'Assuming input grid to be array indices' xx=lindgen(nf+1)+1 if nx gt 0 then yy=x else yy=1L wrr=[xx(0),xx(nf)] y=yy ;ignore Y on input if not keyword_set(forcekey) then ny=n_elements(yy) endif ;handled X) nx=n_elements(xx) if n_elements(wrange) eq 2 then wrr=wr ; figure out default output grid if nbin gt 0 then begin mbin=nbin if n_elements(wrr) eq 0 then wrr=[min(xx),max(xx)] wmin=min(wrr,max=wmax) dw=float(wmax-wmin)/float(mbin) yydef=[wmin+findgen(mbin)*dw,wmax] endif else begin mbin=-nbin & wmin=min(alog10(abs(wrr)>(1e-10)),max=wmax) dw=(wmax-wmin)/mbin yydef=[wmin+findgen(mbin)*dw,wmax] & yydef=10.^(yydef) endelse ; now figure out which output grid to use case ny of 0: begin ;(case 3 -- not case 3 of X yy=yydef & ny=n_elements(yy) & nbin=ny-1L end ;NY=0) 1: begin ;(case 1 nbin=long(y(0)) if nbin eq 0 then nbin=1L ;NBIN cannot be 0 if nbin gt 0 then begin mbin=nbin & wmin=min(wrr,max=wmax) dw=float(wmax-wmin)/float(mbin) yy=[wmin+findgen(mbin)*dw,wmax] endif else begin mbin=-nbin & wmin=min(alog10(abs(wrr)),max=wmax) dw=(wmax-wmin)/mbin yy=[wmin+findgen(mbin)*dw,wmax] & yy=10.^(yy) endelse ny=n_elements(yy) & nbin=ny-1L ;if nbin eq 0 then begin ;(case 3 above ; ny=n_elements(yy) ; if ny eq 1 then nbin=yy(0) else nbin=ny-1 ;endif else yy=y ;NBIN=0) end ;NY=1) else: begin ;(case 2 -- also case 3 of X if not keyword_set(forcekey) then begin yy=y & nbin=ny-1 & wrr=[yy(0),yy(ny-1)] endif else yy=yydef ;if keyword_set(wrange) then message,'ignoring specified WRANGE',/info end ;NY>1) endcase wr=wrr ;if ny eq 1 then begin ;(make regular grid ; if nbin lt 0 then begin ;(log grid ; mbin=-nbin & wmin=min(alog10(abs(wr)),max=wmax) ; dw=(wmax-wmin)/mbin ; yy=[wmin+findgen(mbin)*dw,wmax] & yy=10.^(yy) ; endif else begin ;)(linear grid ; mbin=nbin & wmin=min(wr,max=wmax) ; dw=float(wmax-wmin)/float(mbin) ; yy=[wmin+findgen(mbin)*dw,wmax] ; endelse ;NBIN>0) ;endif ;NY=1) ; get bin widths ff=f[*] & dx=xx(1:*)-xx & dy=yy(1:*)-yy ; make sure everything's in the right order o1=where(dx lt 0,mo1) if mo1 gt 0 then begin oz=sort(xx) xmid=0.5*(xx(1:*)+xx) xx=xx(oz) oz=sort(xmid) & ff=ff(oz) & dx=xx(1:*)-xx endif o1=where(dy lt 0,mo1) if mo1 gt 0 then begin oz=sort(yy) & yy=yy(oz) & dy=yy(1:*)-yy message,'WARNING: On output Y will be sorted in increasing order',\$ /informational endif ; and correct for stupid user tricks o1=where(dx eq 0,mo1) if mo1 gt 0 then begin o2=where(dx gt 0,mo2) if mo2 eq 0 then begin message,"Give us this day a sliced array",/info return,total(ff)+0*yy endif message,'forcing non-uniqueness on the absissae',/info if vv gt 100 then stop,'Halting. Type .CON to continue' mdx=min(dx(o2)) dx(o1)=((spline(o2,dx(o2),o1,0.1)<(mdx/10.))>(mdx/1e6)) ;dx(o1)=min(dx(o2))/10. for i=1L,mo1-1L do xx(o1(i))=xx(o1(i)-1)+dx(o1(i)-1) ; resort to make sure it is monotonic os=sort(xx) xmid=0.5*(xx(1:*)+xx) xx=xx[os] os=sort(xmid) & ff=ff[os] endif if not keyword_set(perbin) then ff=ff*dx ;(units correction ; make cumulative function if float(!version.RELEASE) lt 5.3 then begin cf=0*xx+0*ff(0) for i=0L,nf-1L do cf(i+1)=cf(i)+ff(i) endif else cf=[0.,total(ff,/cumulative)] base=min(cf,max=top) ; interpolate if not keyword_set(slowok) then cg=interpolate(cf,findex(xx,yy)) else \$ cg=interpol(cf,xx,yy) ;(slower) cg = ((cg > base) < top) ;extrapolate ; undo accumulation gg=cg(1:*)-cg y=yy ;output if not keyword_set(perbin) then gg=gg/dy ;units correction) return,gg end