### 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.