Home
Manual
Packages
Global Index
Keywords
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 1D arrays y and x of equal length is:
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 "pointcentered", while the
values deriv can be called "zonecentered". The zcen and
pcen functions provide a simple mechanism for moving back and
forth between pointcentered and zonecentered quantities. Usually,
there will be several reasonable ways to pointcenter zone centered data
and viceversa. For example:
 deriv_pc1 = deriv(pcen);
deriv_pc2 = y(dif)(pcen)/x(dif)(pcen);
deriv_pc3 = y(pcen)(dif)/x(pcen)(dif);

For a wellresolved 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 1D 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 1D 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.
