cases.f

110 lines · 1 definitions · on GitHub (7be4133)

Python version: cases.py, 82 lines, on the right (switch to Python above the code).

Purpose

A helper for shape-preserving quadratic splines, from the package by D. F. McAllister and J. A. Roulier (coded by S. L. Dodd and M. Roulier, N.C. State University, per its header). Given two data points and the slopes there, it computes the knots of the quadratic spline on that interval.

How it fits

Inputs: the end points $(p_1, p_2)$ and $(q_1, q_2)$, the slopes $m_1$ and $m_2$ there, and NCASE (1 to 4), which fixes the number and place of the knots.
Who calls it: MEVAL, which evaluates the spline.
Outputs: the knots $(v_1,v_2)$, $(w_1,w_2)$, $(z_1,z_2)$, and for case 4 also $(y_1,y_2)$ and $(e_1,e_2)$.
No COMMON blocks.

The method

Each case builds the spline from straight segments of the given slopes through the end points. In case 1, for example, the two tangent lines meet at $$z_1 = \frac{p_2 - q_2 + m_2 q_1 - m_1 p_1}{m_2 - m_1},$$ and the knots are the midpoints $v = \tfrac12(p + z)$, $w = \tfrac12(z + q)$. Cases 2–4 place the middle knot at the interval's midpoint or at auxiliary points ($c_1 = p_1 + (q_2 - p_2)/m_1$ and $d_1 = q_1 + (p_2 - q_2)/m_2$), as the slopes require.

The divisions are safe for the cases CHOOSE picks, as long as $q_1 \ne p_1$: case 1 divides by $m_2 - m_1$ (line 55 ↗), and CHOOSE picks it only when the two slopes differ (opposite signs, or the same sign with $|s|$ strictly between $|m_1|$ and $|m_2|$, choose.f, lines 59–61 ↗, choose.f, lines 82–89 ↗, choose.f, lines 131–143 ↗); cases 3 and 4 divide by $m_1$ and $m_2$ (lines 74–75 ↗), which CHOOSE rules out as zero (choose.f, line 57 ↗), by $h_1 - p_1 = 2(q_2 - p_2)/m_1$, nonzero since $s \ne 0$ there, and in case 3 by $\bar m_2 - \bar m_1 = (m_2 - m_1)/2$ (line 84 ↗), nonzero since one slope exceeds $2|s|$ and the other does not exceed $(2-\epsilon)|s|$ (settled from the code, Nov 9, 2026).

Definitions

subroutineCASES

lines 5–110 · GitHub

Knots of a shape-preserving quadratic spline on one interval (McAllister and Roulier), for the four cases of their method.

Called by MEVAL ×4

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE.

1 2C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 3C CASES 4C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 5 SUBROUTINE CASES(P1,P2,M1,M2,Q1,Q2,E1,E2,V1,V2,W1,W2,Z1,Z2, 6 * * Y1,Y2,NCASE) 7C 10/20/26 MHP via Claude 8C code changed to replace implicit typing with IMPLICIT NONE and 9C explicit declarations. 10 IMPLICIT NONE 11 INTEGER NCASE 12C 13C SHAPE PRESERVING QUADRATIC SPLINES 14C BY D.F.MCALLISTER & J.A.ROULIER 15C CODED BY S.L.DODD & M.ROULIER 16C N.C. STATE UNIVERSITY 17C 18 REAL*8 P1,P2,M1,M2,Q1,Q2,V1,V2,Z1,Z2,W1,W2,E1,E2, 19 * * MBAR1,MBAR2,MBAR3,C1,D1,H1,J1,Y1,Y2,K1,ZTWO 20C 21C CASES COMPUTES THE KNOTS AND OTHER PARAMETERS OF THE SPLINE ON THE 22C INTERVAL (P1,Q1). 23C 24C 25C ON INPUT-- 26C 27C (P1,P2) AND (Q1,Q2) ARE THE COORDINATES OF THE POINTS OF 28C INTERPOLATION. 29C 30C M1 IS THE SLOPE AT (P1,P2). 31C 32C M2 IS THE SLOPE AT (Q1,Q2) 33C 34C NCASE CONTROLS THE NUMBER AND LOCATION OF THE KNOTS. 35C 36C 37C ON OUTPUT-- 38C 39C (V1,V2),(W1,W2),(Z1,Z2), AND (E1,E2) ARE THE COORDINATES OF THE 40C KNOTS AND OTHER PARAMETERS OF THE SPLINE ON (P1,Q1). (E1,E2) 41C AND (Y1,Y2) ARE USED ONLY IF NCASE=4. 42C 43C AND 44C 45C CASES DOES NOT ALTER P1,P2,M1,M2,Q1,Q2. 46C 47C---------------------------------------------------------------------- 48C 49 save 50 51 IF ((NCASE .EQ. 3).OR.(NCASE .EQ. 4)) GO TO 20 52 IF (NCASE .EQ. 2) GO TO 10 53C 54C CALCULATE THE PARAMETERS FOR CASE 1. 55 Z1=(P2-Q2+M2*Q1-M1*P1)/(M2-M1) 56 ZTWO=P2+M1*(Z1-P1) 57 V1=(P1+Z1)/2.D0 58 V2=(P2+ZTWO)/2.D0 59 W1=(Z1+Q1)/2.D0 60 W2=(ZTWO+Q2)/2.D0 61 Z2=V2+((W2-V2)/(W1-V1))*(Z1-V1) 62 RETURN 63C 64C CALCULATE THE PARAMETERS FOR CASE 2. 65 10 10 Z1=(P1+Q1)/2.D0 66 V1=(P1+Z1)/2.D0 67 V2=P2+M1*(V1-P1) 68 W1=(Z1+Q1)/2.D0 69 W2=Q2+M2*(W1-Q1) 70 Z2=(V2+W2)/2.D0 71 RETURN 72C 73C CALCULATE THE PARAMETERS USED IN BOTH CASES 3 AND 4. 74 20 20 C1=P1+(Q2-P2)/M1 75 D1=Q1+(P2-Q2)/M2 76 H1=2.D0*C1-P1 77 J1=2.D0*D1-Q1 78 MBAR1=(Q2-P2)/(H1-P1) 79 MBAR2=(P2-Q2)/(J1-Q1) 80C 81 IF (NCASE .EQ. 4) GO TO 50 82C 83C CALCULATE THE PARAMETERS FOR CASE 3. 84 K1=(P2-Q2+Q1*MBAR2-P1*MBAR1)/(MBAR2-MBAR1) 85 IF (ABS(M1) .GT. ABS(M2)) GO TO 30 86 Z1=(K1+Q1)/2.D0 87 GO TO 40 88 30 30 Z1=(K1+P1)/2.D0 89 40 40 V1=(P1+Z1)/2.D0 90 V2=P2+M1*(V1-P1) 91 W1=(Q1+Z1)/2.D0 92 W2=Q2+M2*(W1-Q1) 93 Z2=V2+((W2-V2)/(W1-V1))*(Z1-V1) 94 RETURN 95C 96C CALCULATE THE PARAMETERS FOR CASE 4. 97 50 50 Y1=(P1+C1)/2.D0 98 V1=(P1+Y1)/2.D0 99 V2=M1*(V1-P1) + P2 100 Z1=(D1+Q1)/2.D0 101 W1=(Q1+Z1)/2.D0 102 W2=M2*(W1-Q1) + Q2 103 MBAR3=(W2-V2)/(W1-V1) 104 Y2=MBAR3*(Y1-V1) + V2 105 Z2=MBAR3*(Z1-V1) + V2 106 E1=(Y1+Z1)/2.D0 107 E2=MBAR3*(E1-V1) + V2 108 RETURN 109C 110 END

1"""Python version of cases.f (YREC 7be4133): the knots of a shape-preserving 2quadratic spline on one interval (McAllister and Roulier).""" 3 4 5def cases(p1: float, p2: float, m1: float, m2: float, q1: float, q2: float, e1: float, e2: float, v1: float, v2: float, 6 w1: float, w2: float, z1: float, z2: float, y1: float, y2: float, ncase: int) -> tuple[float, float, float, 7 float, float, float, float, float, float, float]: 8 """The knots of the quadratic spline on (p1, q1), for case ncase (1 to 9 4) of McAllister and Roulier's method. 10 11 Translates: cases.f#CASES (7be4133) 12 13 Arguments: (p1, p2) and (q1, q2), the points of interpolation; m1, m2, 14 the slopes there; ncase, which fixes the number and place of the 15 knots. The outputs e1 ... y2 are passed in (with any value). 16 17 Returns (e1, e2, v1, v2, w1, w2, z1, z2, y1, y2): the knots (v1, v2), 18 (w1, w2), (z1, z2), and for case 4 also (e1, e2) and (y1, y2); those a 19 case does not set are returned as passed in. 20 21 """ 22 23 # The divisions by m2 - m1, m1, and m2 are safe for the cases CHOOSE 24 # picks (see the page). 25 if ncase not in (3, 4): 26 if ncase == 2: 27 # Label 10: case 2 (Fortran lines 65–71). 28 z1 = (p1 + q1) / 2.0 29 v1 = (p1 + z1) / 2.0 30 v2 = p2 + m1 * (v1 - p1) 31 w1 = (z1 + q1) / 2.0 32 w2 = q2 + m2 * (w1 - q1) 33 z2 = (v2 + w2) / 2.0 34 return e1, e2, v1, v2, w1, w2, z1, z2, y1, y2 35 36 # Case 1, also any NCASE other than 2, 3, 4 (Fortran lines 51–62). 37 z1 = (p2 - q2 + m2 * q1 - m1 * p1) / (m2 - m1) 38 ztwo = p2 + m1 * (z1 - p1) 39 v1 = (p1 + z1) / 2.0 40 v2 = (p2 + ztwo) / 2.0 41 w1 = (z1 + q1) / 2.0 42 w2 = (ztwo + q2) / 2.0 43 z2 = v2 + ((w2 - v2) / (w1 - v1)) * (z1 - v1) 44 return e1, e2, v1, v2, w1, w2, z1, z2, y1, y2 45 46 # Label 20: the parameters used in both cases 3 and 4 (Fortran lines 47 # 74–81). 48 c1 = p1 + (q2 - p2) / m1 49 d1 = q1 + (p2 - q2) / m2 50 h1 = 2.0 * c1 - p1 51 j1 = 2.0 * d1 - q1 52 mbar1 = (q2 - p2) / (h1 - p1) 53 mbar2 = (p2 - q2) / (j1 - q1) 54 55 # Case 3 (Fortran lines 84–94). 56 if ncase == 3: 57 k1 = (p2 - q2 + q1 * mbar2 - p1 * mbar1) / (mbar2 - mbar1) 58 if abs(m1) > abs(m2): 59 z1 = (k1 + p1) / 2.0 # label 30 60 else: 61 z1 = (k1 + q1) / 2.0 62 # label 40 63 v1 = (p1 + z1) / 2.0 64 v2 = p2 + m1 * (v1 - p1) 65 w1 = (q1 + z1) / 2.0 66 w2 = q2 + m2 * (w1 - q1) 67 z2 = v2 + ((w2 - v2) / (w1 - v1)) * (z1 - v1) 68 return e1, e2, v1, v2, w1, w2, z1, z2, y1, y2 69 70 # Label 50: case 4 (Fortran lines 97–108). 71 y1 = (p1 + c1) / 2.0 72 v1 = (p1 + y1) / 2.0 73 v2 = m1 * (v1 - p1) + p2 74 z1 = (d1 + q1) / 2.0 75 w1 = (q1 + z1) / 2.0 76 w2 = m2 * (w1 - q1) + q2 77 mbar3 = (w2 - v2) / (w1 - v1) 78 y2 = mbar3 * (y1 - v1) + v2 79 z2 = mbar3 * (z1 - v1) + v2 80 e1 = (y1 + z1) / 2.0 81 e2 = mbar3 * (e1 - v1) + v2 82 return e1, e2, v1, v2, w1, w2, z1, z2, y1, y2