Equiprobable bins

For pairs of values (*x*, *P*(*x*)) with *P* a probability density, the
routine `eqpb` uses linear interpolation to compute `equi_dist`, an equidistributed probability density. There are two
basic applications of this routine: angular distributions (
= 1) data and distributions of energy carried off by
secondary particles (
= 4 data). But there is also
the
= 3 data for correlated distributions of angle
and energy. We begin with a description of the general method and
discuss later the special provisions made for the correlated
distributions for
= 3.

The method used is as follows. Suppose that we want
bins, each with probability
1/, with the bottom of the
lowest bin and the top of the highest bin fixed by the lowest and
highest input values of *x*. The equidistributed interior bin
boundaries are located at the *x*-values for which the accumulated
probabilities are
*k*/ for *k* = 1, 2, ...,
- 1. The geometry of the problem is that we want the set of *x*
points which give equal area under the graph of *P*. From a
mathematical point of view, what we want is the points (for *k* = 0, 1, ..., ` nbins`) such that is the first
library

Suppose that we seek the *k*-th equidistribution point, . We
start by adding the areas of the trapezoids given by the library (*x*, *P*(*x*)) data, and we eventually find the first time that the
accumulated sum exceeds the desired area
*k*/. This
brackets between two of the library *x*-values, and *P* is
linear on this interval. The integral of *P* is therefore quadratic
on the interval, and the determination of requires the
solution of a quadratic equation. Only one root of the quadratic
makes physical sense, and the algorithm takes advantage of this fact.

The extraneous root is easily understood by looking at the Fig. 1. In
that figure suppose that *B* is the equidistribution point we are
looking for, and suppose we know that *B* lies between the two library
*x*-values *A* and *C*. The library probability distribution is a
broken-line function *P*(*x*) which is equal to a linear function *L*(*x*)
on the interval
*A**x**C*. In the figure we have extrapolated
*L*(*x*) as a dashed line out to *D*. Note that