slopes.f

Computes the derivative at each node for the shape-preserving quadratic splines of McAllister and Roulier (coded by S. L. Dodd and M. Roulier, N.C. State University). Its comment: "the slopes provided will insure that an osculatory quadratic spline will have one additional knot between two adjacent points of interpolation. Convexity and monotonicity are preserved wherever these conditions are compatible with the data."

Called by OSPLIN before MEVAL evaluates the spline (with SEARCH). It contains a second routine, SAFEDIVIDE, added by KC. No COMMON blocks.

At each interior node $i$, with the slopes $m_1$, $m_2$ of the chords to its two neighbors:

• if $m_1 m_2 \le 0$ (a local extremum or a flat chord), the derivative is 0;
• otherwise the chord with the larger slope is extended until it reaches the level of the far neighbor, at $\bar x$, and the derivative is the slope of the line from the node to the far neighbor's level at $\hat x$, the midpoint between $\bar x$ and the far node. For $|m_1| > |m_2|$: $$\bar x = x_i + \frac{y_{i+1} - y_i}{m_1},\qquad \hat x = \tfrac12(\bar x + x_{i+1}),\qquad y'_i = \frac{y_{i+1} - y_i}{\hat x - x_i}.$$

At each end node, the derivative is the slope from the midpoint of the end interval, where it takes the value of the line through the neighbor node with that node's derivative, to the end node. It is set to 0 if that would reverse the sign of the end chord, or to $2m$, twice the end chord's slope $m$, if the two chords nearest that end have opposite signs.

Two divisions go through SAFEDIVIDE, which returns 0 instead of dividing by zero: the slope $m_2$ in the loop (line 96↗) and the derivative at the last node (line 105↗). A zero denominator means two equal abscissae, which would make the spline ill-defined anyway, so the 0 hides the problem, and a check in OSPLIN might be clearer; which caller can pass a table with two equal abscissae is CHOOSE's question 1.