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 x-value and
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 AxC. In the figure we have extrapolated L(x) as a dashed line out to D. Note that