; docformat = 'rst' ;+ ; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010 ;- ;+ ; :Description: ; Calculates the area under the y(x) curve provided. The function must be ordered in x, but it does not ; matter whether the order is increasing or decreasing. ; ; There are two possible ways that the input function is interpreted, determined by the keyword local. ; See the description of local below, since it may significantly alter the result. ; ; If the method selected requires interpolation, it can be linear, 3 or 4-point quadratic, or spline. This is determined ; by setting the quadratic, lsquadratic, or spline keywords (linear if none of these 3 keywords is set). ; ; :Returns: ; If xmin and xmax are not provided, returns the area under the curve over its whole extension. ; ; If xmin and xmax are scalars, returns the area between these x values. If they are arrays of n elements, ; returns the n areas, calculated starting at each xmin element, and ending at the corresponding xmax element. ; ; Every value in xmax is must be larger than or equal to the correspoinding value in xmin. ; ; :Params: ; x : in, required ; An array of locations where the function is sampled. Must be ordered (increasing or decreasing). ; y : in, required ; An array with the function values corresponding to the locations in x. ; :Keywords: ; xmin : in, out, optional ; A scalar or array with the start of the range(s) where the area is to be calculated. If not provided, ; the minimum location of the function is used. Must have the same number of elements as xmax, and every ; element in xmin must be smaller than or equal to the corresponding element in xmax. If any value of xmin ; is smaller than the beginning of the x range, it is clipped to the beginning. ; xmax : in, out, optional ; A scalar or array with the end of the range(s) where the area is to be calculated. If not provided, ; the maximum location of the function is used. Must have the same number of elements as xmin, and every ; element in xmax must be larger than or equal to the corresponding element in xmin. If any value of xmax ; is larger than the end of the x range, it is clipped to the end. ; cumulative : in, optional, default=0 ; If set, and xmin and xmax are not provided, the areas returned are the cumulative areas at the end of each x point, ; starting at the first x point. If xmax and xmin are provided, this keyword is ignored. ; xinc : out, optional ; Returns the x values, in increasing order. ; yinc : out, optional ; Returns the y values, in order of increasing x. ; xstart : out, optional ; If local is not set, returns the location where each input pixel starts. If local is set, this is ignored. ; xend : out, optional ; If local is not set, returns the location where each input pixel ends. If local is set, this is ignored. ; xbox : out, optional ; Used when not in local mode, to return the input x values in pairs of equal values, so that ; plotting xbox,ybox shows the rectangles that are how the input values were interpreted. ; ybox : out, optional ; Used when not in local mode, to return the input y values in pairs of equal values, so that ; plotting xbox,ybox shows the rectangles that are how the input values were interpreted. ; newton : in, optional, default=0 ; If local is set and newton is set, the function integrations are done with int_tabulated, which ; uses a 5 point Newton-Cotes formula. If not set, trapezoid integration is done. If local is not set, this has no effect. ; local : in, optional, default=0 ; Determines how the input function is interpreted, which determines how the areas are calculated. If set, the function is ; interpreted literally, that is, as in a simple mathematical function: the y values are the local evaluation y(x). ; ; If not set, the y values are interpreted as a measured sample, that is, as the average of the "flux" falling inside the region ; centered at the corresponding x values: x locations are interpreted as the centers of bins. Since this interpretation is equivalent ; to a literal interpretation of a function made of rectangular regions, in this method there are no interpolations (partial areas, if ; any, are just the corresponding fractions of the rectangles). ; lsquadratic : in, optional, default=0 ; If the method requires interpolation, this is passed on to interpol, to select 4 point quadratic interpolation. ; ; If none of lsquadratic, quadratic and spline are set, interpolation is linear. ; quadratic : in, optional, default=0 ; If the method requires interpolation, this is passed on to interpol, to select 3 point quadratic interpolation. ; ; If none of lsquadratic, quadratic and spline are set, interpolation is linear. ; spline : in, optional, default=0 ; If the method requires interpolation, this is passed on to interpol, to select spline interpolation. ; ; If none of lsquadratic, quadratic and spline are set, interpolation is linear. ; ; ; :Examples: ; ; Make up a simple constant function to show the difference between local and non-local modes:: ; ; x=[0d0,1d0] ; y=[2d0,2d0] ; print,pp_integral(x,y,xmin=xmin,xmax=xmax) ; ;4.0000000 ; print,xmin,xmax ; ;-0.50000000 1.5000000 ; ;(the x locations are considered the middle of the bins, so the bins extend beyond the range min(x),max(x)) ; print,pp_integral(x,y,xmin=lxmin,xmax=lxmax,/local) ; ;2.0000000 ; print,lxmin,lxmax ; ;0.0000000 1.0000000 ; ; Make up a well sampled function, so that the difference between local and non-local is small:: ; ; x=dindgen(10001)*!dpi/1d4 ; y=sin(x) ; a=pp_integral(x,y,xmin=[0d0,0d0],xmax=!dpi*[0.5d0,1d0]) ; b=pp_integral(x,y,xmin=[0d0,0d0],xmax=!dpi*[0.5d0,1d0],/local) ; print,a ; ;1.0003142 2.0000000 ; print,b ; ;1.0000000 2.0000000 ; print,a-b ; ;0.00031410992 -4.9348009e-08 ; ; ; ; :Author: Paulo Penteado (pp.penteado@gmail.com), Feb/2010 ;- function pp_integral,x,y,xmin=xmin,xmax=xmax,cumulative=cumulative,$ xinc=ix,yinc=iy,xstart=xstart,xend=xend,xbox=xbox,ybox=ybox,$ newton=newton,local=local,$ lsquadratic=lsquadratic,quadratic=quadratic,spline=spline compile_opt idl2 ;Check that the dimensions of x and y match and are acceptable nx=n_elements(x) & ny=n_elements(y) if (nx ne ny) || (nx eq 0) then message,'x and y provided must be arrays of the same length' ;Check that the dimensions of xmin and xmax match and are acceptable nxmin=n_elements(xmin) & nxmax=n_elements(xmax) if (nxmin ne nxmax) then message,'xmin and xmax provided must be arrays of the same length' ;Defaults local=n_elements(local) eq 1 ? local : 0 newton=local ? 0 : n_elements(newton) eq 1 ? newton : 0 cumulative=n_elements(cumulative) eq 1 ? cumulative : 0 ;Make ix and iy to work with, 1D and in increasing order (x is assumed sorted, but can be increasing or decreasing order) if (x[nx-1] lt x[0]) then begin ix=reverse(reform(x)) & iy=reverse(reform(y)) endif else begin ix=reform(x) & iy=reform(y) endelse if (~local) then begin ;Input function is histogram-like (rectangle rule) ;Calculate the x start and end points for the tabulated function xstart=[ix[0]-(ix[1]-ix[0])/2d0,(ix[0:nx-2]+ix[1:nx-1])/2d0] xend=[xstart[1:nx-1],ix[nx-1]+(ix[nx-1]-ix[nx-2])/2d0] ;Default interval xmin=n_elements(xmin) gt 0 ? xmin : cumulative ? replicate(xstart[0],nx) : xstart[0] xmax=n_elements(xmax) gt 0 ? xmax : cumulative ? xend : xend[nx-1] ;Clip the edges of the interval if they exceeed the input function's range xmin=xstart[0]>xmin<xend[nx-1] xmax=xstart[0]>xmax<xend[nx-1] nlocs=n_elements(xmin) minloc=value_locate(xstart,xmin) maxloc=(value_locate(xend,xmax)+1)<(nx-1) ;Area of the partial rectangles at the edges ret=(xend[minloc]-xmin)*iy[minloc]+(xmax-xend[maxloc])*iy[maxloc] ;Area where min and max fall in the same rectangle w=where((maxloc-minloc) eq 0L,nw) if (nw gt 0L) then ret[w]=(xmax[w]-xmin[w])*iy[maxloc[w]] ;Area of the full rectangles contained in the output intervals w=where((maxloc-minloc) gt 0L,nw) if (nw gt 0L) then begin areas=total((xend-xstart)*iy,/cumulative) ;The cumulative area of the rectangles ret[w]+=areas[maxloc[w]]-areas[minloc[w]] endif endif else begin ;Input function is interpreted literally ;Default interval xmin=n_elements(xmin) gt 0 ? xmin : cumulative ? replicate(ix[0],nx) : ix[0] xmax=n_elements(xmax) gt 0 ? xmax : cumulative ? ix : ix[nx-1] ;Clip the edges of the interval if they exceeed the input function's range xmin=ix[0]>xmin<ix[nx-1] xmax=ix[0]>xmax<ix[nx-1] minloc=value_locate(ix,xmin) maxloc=value_locate(ix,xmax) ymins=interpol(iy,ix,xmin,lsquadratic=lsquadratic,quadratic=quadratic,spline=spline) ymaxs=interpol(iy,ix,xmax,lsquadratic=lsquadratic,quadratic=quadratic,spline=spline) if (newton) then begin ;Use int_tabulated, that does a 5 point Newton-Cotes integration nlocs=n_elements(xmin) ret=dblarr(nlocs) for i=0L,nlocs-1 do begin xtmp=[xmin[i],ix[minloc[i]:maxloc[i]:minloc[i] le maxloc[i] ? 1L : -1L],xmax[i]] ytmp=[ymins[i],iy[minloc[i]:maxloc[i]:minloc[i] le maxloc[i] ? 1L : -1L],ymaxs[i]] ret[i]=int_tabulated(xtmp,ytmp,/double) endfor endif else begin ;Use trapezoidal integration ;Area of the partial trapezoids at the edges ret=(ix[(minloc+1)<(nx-1)]-xmin)*(ymins+iy[(minloc+1)<(nx-1)])*0.5d0 w=where(ix[minloc] eq xmin,nw) ; if (nw gt 0L) then ret[w]=0d0 ret+=(xmax-ix[maxloc])*(ymaxs+iy[maxloc])*0.5d0 ;Area where min and max fall in the same trapezoid w=where((maxloc-minloc) eq 0L,nw) if (nw gt 0L) then ret[w]=(xmax[w]-xmin[w])*(ymaxs[w]+ymins[w])*0.5d0 ;Area of the full trapezoids contained in the output intervals w=where((maxloc-minloc) gt 1L,nw) if (nw gt 0L) then begin areas=[0d0,total((ix[1:*]-ix[0:nx-1])*(iy[1:*]+iy[0:nx-1])*0.5d0,/cumulative)] ;The cumulative area of the trapezoids ret[w]+=areas[maxloc[w]]-areas[(minloc[w]+1)<(nx-1)] endif endelse endelse ;Create xbox and ybox arrays, with constant y values spanning the start and end of each x point if (~local && (arg_present(xbox) || arg_present(ybox))) then begin nix=n_elements(ix) xbox=dblarr(2L*nix) & ybox=dblarr(2L*nix) xbox[0:*:2]=xstart & xbox[1:*:2]=xend ybox[0:*:2]=iy & ybox[1:*:2]=iy endif return,ret end