Integrates a sequence of unequally or equally spaced points using Simpson's three-point formula.
function simpne ( x : numeric, y : numeric ) return_val : float or double
An array of one or more dimensions. The rightmost dimension is the dimension to be used for the integration. It must be strictly monotonically increasing. If y is multi-dimensional and x is one-dimensional then x must be the same size as the rightmost dimension of y. Otherwise, x and y must be the same size.y
An array of one or more dimensions. The rightmost dimension is the dimension to be integrated. This may require that the original variable be reordered. Given that Simpson's 3-point formula is used, there must be 3 or more non-missing values present.
The return value will have same dimensions as all but the rightmost dimension of y (or be a scalar if y is one-dimensional). The return type will be double if y is double, and float otherwise.
This function uses Simpson's 3-point formula for integration adapted for unequally sampled points. The number of points may be odd or even. The more points (finer resolution) the more accurate the answer. The x must be strictly monotonically increasing. Missing values are allowed. They are ignored.
The more points (finer resolution) the more accurate the answer.
Example 1 Evaluate exp(x) from x=1.8 to x=3.4.
x = (/1.8, 2.0, 2.33, 2.49, 3.122, 3.4/) ; even number pts (6) y = exp(x) simpne_6 = simpne(x,y) ; = 23.9757 ; even number pts (6) X = (/1.8, 2.0, 2.33, 2.49, 2.75, 3.122, 3.4/) Y = exp(X) simpne_7 = simpne(X,Y) ; = 23.911 ; odd number pts (7)The use of a function (here, exp) is for convenience. The fi (here, y) could be from a table or series of unequally spaced values.
Example 2: Assume x is dimensioned ntim and y is dimensioned nloc x ntim. Integrate over time at each location:
result = simpne(x,y) ; result(nloc)
Assume z is dimensioned lev and f is dimensioned ntim x klev x nlat x mlon with named dimensions "time", "lev", "lat", and "lon". Integrate over z:
result = simpne(z,f(time|:,lat|:,lon|:,lev|:)) ; result(ntim,nlat,mlon)
The values of x must be strictly monotonically increasing. Specifically, this means that every x(i+1).gt.x(i). Consider x[*] and y[*]:
if (isMonotonic(x).eq.1) then ; =1 means strictly monotonic result = simpne(x,y) else nx = dimsizes(x) imono = ind((x(1:)-x(0:nx-2)).gt.0) result = simpne(x(imono),y(imono)) end if