Quick Reference

2.3.10 Rank preserving (finite difference) range functions

Because Yorick arrays almost invariably represent function values, Yorick provides numercal equivalents to the common operations of differential and integral calculus. In order to handle functions of several variables in a straightforward manner, these operators are implemented as range functions. Unlike the statistical range functions, which return a scalar result, the finite difference range functions do not reduce the rank of the subscripted array. Instead, they preserve rank, in the same way that a simple index range start:stop preserves rank. The available finite difference functions are:

cum

returns the cumulative sums (sum of all values at smaller index, the first index of the result having a value of 0). The result dimension has a length one greater than the dimension in the subscripted array.

psum

returns the partial sum (sum of all values at equal or smaller index) up to each element. The result dimension has the same length as the dimension in the subscripted array. This is the same as cum, except that the leading 0 value in the result is omitted.

dif

returns the pairwise difference between successive elements, positive when the function values increase along the dimension, negative when they decrease. The result dimension has a length one less than the dimension of the subscripted array. The dif function is the inverse of the cum function.

zcen

returns the pairwise average between successive elements. The result dimension has a length one less than the dimension of the subscripted array. The name is short for "zone center".

pcen

returns the pairwise average between successive elements. However, the two endpoints are copied unchanged, so that the result dimension has a length one greater than the dimension of the subscripted array. The name is short for "point center".

uncp

returns the inverse of the pcen operation. (Gibberish will be the result of applying this function to an array dimension which was not produced by the pcen operation.) The result dimension has a length one less than the dimension of the subscripted array. The name is short for "uncenter point centered".

The derivative dy/dx of a function y(x), where y and x are represented by 1-D arrays y and x of equal length is:

deriv = y(dif)/x(dif);

Note that deriv has one fewer element than either y or x. The derivative is computed as if y(x) were the piecewise linear function passing through the given points (x,y); there is one fewer line segment (slope) than point.

The values x and y can be called "point-centered", while the values deriv can be called "zone-centered". The zcen and pcen functions provide a simple mechanism for moving back and forth between point-centered and zone-centered quantities. Usually, there will be several reasonable ways to point-center zone centered data and vice-versa. For example:

deriv_pc1 = deriv(pcen); deriv_pc2 = y(dif)(pcen)/x(dif)(pcen); deriv_pc3 = y(pcen)(dif)/x(pcen)(dif);

For a well-resolved function, the differences among these three arrays will be negligible. That is, the differences are second order in x(dif), which is often the order of the errors in the calculation that produced x and y in the first place. If x and y represent y(x) more accurately that this, then you must know a better model of the shape of y(x) than the simple piecewise linear model, and you should use that model to select deriv_pc1, deriv_pc2, and deriv_pc3, or some other expression.

An indefinite integral may be estimated using the trapezoidal rule:

indef_integ = (y(zcen)*x(dif))(psum);

Once again, indef_integ has one fewer point than x or y, because there is one fewer trapezoid than point. This time, however, indef_integ is not zone centered. Instead, indef_integ represents values at the upper (or lower) boundaries x(2:) (or x(:-1)). Often, you want to think of the integral of y(x) as a point centered array of definite integrals from x(1) up to each of the x(i). In this case (which actually arises more frequently), use the cum function instead of psum in order to produce a result def_integs with the same number of points as the x and y arrays:

def_integs = (y(zcen)*x(dif))(cum);

For single definite integrals, the matrix multiply syntax can be used in conjunction with the dif range function. For example, suppose you know the transmission fraction of the filter, ff, at several photon energies ef. That is, ff and ef are 1-D arrays of equal length, specifying a filter transmission function as the piecewise linear function connecting the given points. The final dimension of an array brightness represents an incident spectrum (any other dimensions represent different rays, say one ray per pixel of an imaging detector). The 1-D array gb represents the group boundary energies -- that is, the photon energies at the boundaries of the spectral groups represented in brightness. The following Yorick statements compute the detector response:

filter = integ(ff, ef, gb)(dif); response = brightness(..,+)*filter(+);

For this application, the correct interpolation routine is integ. The integrated transmission function is evaluated at the boundaries gb of the groups; the effective transmission fraction for each group is the difference between the integral at the upper bin boundary and the lower bin boundary. The dif range function computes these pairwise differences. Since the gb array naturally has one more element than the final dimension of brightness (there is one more group boundary energy than group), and since dif reduces the length of a dimension by one, the filter array has the same length as the final dimension of brightness.

Note that the cum function cannot be used to integrate ff(ef), because the points gb at which the integral must be evaluated are not the same as the points ef at which the integrand is known. Whenever the points at which the integral is required are the same (or a subset of) the points at which the integrand is known, you should perform the integral using the zcen, dif, and cum index functions, instead of the more general integ interpolator.