spline.f

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

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

Purpose

SPLINE evaluates a shape-preserving quadratic spline at one point, in one interval of the table. Its header credits the method to D. F. McAllister and J. A. Roulier, coded by S. L. Dodd and M. Roulier (N.C. State University). The knots and parameters of the spline in the interval are computed beforehand.

How it fits

Who calls it: only MEVAL, which controls the evaluation of the osculatory quadratic spline and computes the knots and parameters (through CHOOSE, by MEVAL's comments).
COMMON blocks: none; a function of its arguments only.

The method

In the interval $[x_0, x_1]$ (XTABS, XTABS1), with values $f_0, f_1$ (YTABS, YTABS1), the spline is made of quadratic pieces joined at one knot $z_1$ (Z1; cases 1–3) or two knots $\xi_1 < z_1$ (Y1, Z1; case 4, NCASE = 4). Each piece is a quadratic Bézier form; for example, for $x$ between $z_1$ and $x_1$:

$$s(x) = \frac{z_2\,(x_1 - x)^2 + 2 w_2\,(x - z_1)(x_1 - x) + f_1\,(x - z_1)^2}{(x_1 - z_1)^2},$$

where $z_2$ (Z2) is the value at the knot and $w_2$ (W2) the middle control value. For $x$ between $x_0$ and the knot, V2 plays the role of W2; between the two knots of case 4, E2 does, with Y2 the value at the first knot.

Design

Where a denominator would vanish (a knot at an end of the interval), linear interpolation between the interval's ends is used instead: "MODIFICATION DUE TO MARC PINSONNEAULT 6/87." Each fallback is reached only when XVALS lies outside [XTABS, XTABS1] (with the knot at XTABS1, only for XVALS > XTABS1, lines 66–73 ↗ and lines 139–146 ↗; with the knot at XTABS, only for XVALS < XTABS, lines 79–86 ↗ and lines 102–109 ↗), and the one between the two knots of case 4 (lines 126–133 ↗) never, since it needs Y1 < XVALS < Z1 with Z1 = Y1. MEVAL passes such points only when it extrapolates beyond the table's ends (meval.f, line 135 ↗, meval.f, line 259 ↗): The fallback is then a linear extrapolation of the end interval (settled from the code, Nov 9, 2026).

Definitions

functionSPLINE

lines 1–148 · GitHub

The value of a shape-preserving quadratic spline at XVALS, from one or two knots in the interval (equation above).

Called by MEVAL ×3

Compiler warnings (8)

gfortran -Wall -Wextra, compiling only; nothing was run.

line 60: Equality comparison for REAL(8) [-Wcompare-reals]
line 67: Inequality comparison for REAL(8) [-Wcompare-reals]
line 80: Inequality comparison for REAL(8) [-Wcompare-reals]
line 96: Equality comparison for REAL(8) [-Wcompare-reals]
line 103: Inequality comparison for REAL(8) [-Wcompare-reals]
line 120: Equality comparison for REAL(8) [-Wcompare-reals]
line 127: Inequality comparison for REAL(8) [-Wcompare-reals]
line 140: Inequality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE.

1 FUNCTION SPLINE(XVALS,Z1,Z2,XTABS,YTABS,XTABS1,YTABS1, 2 * * Y1,Y2,E2,W2,V2,NCASE) 3C 10/20/26 MHP via Claude 4C code changed to replace implicit typing with IMPLICIT NONE and 5C explicit declarations. 6 IMPLICIT NONE 7 INTEGER NCASE 8C 9C SHAPE PRESERVING QUDRATIC SPLINES 10C BY D.F.MCALLISTER & J.A.ROULIER 11C CODED BY S.L.DODD & M.ROULIER 12C N.C. STATE UNIVERSITY 13C 14 REAL*8 XVALS,Z1,Z2,XTABS,YTABS,XTABS1,YTABS1,V2, 15 * * W2,Y1,Y2,E2,SPLINE,FS 16C 17C SPLINE FINDS THE IMAGE OF A POINT IN XVAL. 18C 19C ON INPUT-- 20C 21C XVALS CONTAINS THE VALUE AT WHICH THE SPLINE IS EVALUATED. 22C 23C (XTABS,YTABS) ARE THE COORDINATES OF THE LEFT-HAND DATA POINT 24C USED IN THE EVALUATION OF XVALS. 25C 26C (XTABS1,YTABS1) ARE THE COORDINATES OF THE RIGHT-HAND DATA POINT 27C USED IN THE EVALUATION OF XVALS. 28C 29C Z1,Z2,Y1,Y2,E2,W2,V2 ARE THE PARAMETERS OF THE SPLINE. 30C 31C NCASE CONTROLS THE EVALUATION OF THE SPLINE BY INDICATING WHETHER 32C ONE OR TWO KNOTS WERE PLACED IN THE INTERVAL (XTABS,XTABS1). 33C 34C 35C ON OUTPUT-- 36C 37C SPLINE IS THE IMAGE OF XVALS. 38C 39C AND 40C 41C SPLINE DOES NOT ALTER ANY OF THE INPUT PARAMETERS. 42C 43C---------------------------------------------------------------------- 44C 45C *****MODIFICATION DUE TO MARC PINSONNEAULT 6/87***** 46C IF DIVISION BY ZERO WOULD BE CAUSED,LINEAR INTERPOLATION IS USED 47C INSTEAD OF THE SPLINE. 48C IF NCASE .EQ. 4, MORE THAN ONE KNOT WAS PLACED IN THE INTERVAL. 49 save 50 51 IF (NCASE .EQ. 4) GO TO 40 52C 53C CASES 1,2, OR 3. 54C 55C DETERMINE THE LOCATION OF XVALS RELATIVE TO THE KNOT. 56C KC 2025-05-30 fixed "Arithmetic IF statement" 57C IF (Z1 - XVALS) 10,20,30 58 IF (Z1 .LT. XVALS) THEN 59 GOTO 10 60 ELSE IF (Z1 .EQ. XVALS) THEN 61 GOTO 20 62 ELSE 63 GOTO 30 64 END IF 65C 66 10 10 CONTINUE 67 IF(XTABS1.NE.Z1)THEN 68 SPLINE=(Z2*(XTABS1-XVALS)**2+W2*2.D0*(XVALS-Z1)*(XTABS1-XVALS) 69 * * +YTABS1*(XVALS-Z1)**2)/(XTABS1-Z1)**2 70 ELSE 71 FS = (XVALS - XTABS)/(XTABS1 - XTABS) 72 SPLINE = YTABS + FS*(YTABS1 - YTABS) 73 ENDIF 74 RETURN 75C 76 20 20 SPLINE=Z2 77 RETURN 78C 79 30 30 CONTINUE 80 IF(Z1.NE.XTABS)THEN 81 SPLINE=(YTABS*(Z1-XVALS)**2+V2*2.D0*(XVALS-XTABS)*(Z1-XVALS)+ 82 * * Z2*(XVALS-XTABS)**2)/(Z1-XTABS)**2 83 ELSE 84 FS = (XVALS - XTABS)/(XTABS1 - XTABS) 85 SPLINE = YTABS + FS*(YTABS1 - YTABS) 86 ENDIF 87 RETURN 88C 89C CASE 4. 90C 91C DETERMINE THE LOCATION OF XVALS RELATIVE TO THE FIRST KNOT. 92C KC 2025-05-30 fixed "Arithmetic IF statement" 93C 40 IF (Y1 - XVALS) 70,60,50 94 40 40 IF (Y1 .LT. XVALS) THEN 95 GOTO 70 96 ELSE IF (Y1 .EQ. XVALS) THEN 97 GOTO 60 98 ELSE 99 GOTO 50 100 END IF 101C 102 50 50 CONTINUE 103 IF(Y1.NE.XTABS)THEN 104 SPLINE=(YTABS*(Y1-XVALS)**2+V2*2.D0*(XVALS-XTABS)*(Y1-XVALS)+ 105 * * Y2*(XVALS-XTABS)**2)/(Y1-XTABS)**2 106 ELSE 107 FS = (XVALS - XTABS)/(XTABS1 - XTABS) 108 SPLINE = YTABS + FS*(YTABS1 - YTABS) 109 ENDIF 110 RETURN 111C 112 60 60 SPLINE=Y2 113 RETURN 114C 115C DETERMINE THE LOCATION OF XVALS RELATIVE TO THE SECOND KNOT. 116C KC 2025-05-30 fixed "Arithmetic IF statement" 117C 70 IF (Z1 - XVALS) 100,90,80 118 70 70 IF (Z1 .LT. XVALS) THEN 119 GOTO 100 120 ELSE IF (Z1 .EQ. XVALS) THEN 121 GOTO 90 122 ELSE 123 GOTO 80 124 END IF 125C 126 80 80 CONTINUE 127 IF(Z1.NE.Y1)THEN 128 SPLINE=(Y2*(Z1-XVALS)**2+E2*2.D0*(XVALS-Y1)*(Z1-XVALS)+Z2*(XVALS 129 * * -Y1)**2)/(Z1-Y1)**2 130 ELSE 131 FS = (XVALS - XTABS)/(XTABS1 - XTABS) 132 SPLINE = YTABS + FS*(YTABS1 - YTABS) 133 ENDIF 134 RETURN 135C 136 90 90 SPLINE=Z2 137 RETURN 138C 139 100 100 CONTINUE 140 IF(XTABS1.NE.Z1)THEN 141 SPLINE=(Z2*(XTABS1-XVALS)**2+W2*2.D0*(XVALS-Z1)*(XTABS1-XVALS) 142 * * +YTABS1*(XVALS-Z1)**2)/(XTABS1-Z1)**2 143 ELSE 144 FS = (XVALS - XTABS)/(XTABS1 - XTABS) 145 SPLINE = YTABS + FS*(YTABS1 - YTABS) 146 ENDIF 147 RETURN 148 END

1"""Python version of spline.f (YREC 7be4133): the value of a 2shape-preserving quadratic spline (McAllister and Roulier) at one point of 3one interval.""" 4 5 6def _linear(xvals: float, xtabs: float, ytabs: float, xtabs1: float, ytabs1: float) -> float: 7 """The fallback where a denominator would vanish (MHP 6/87): linear 8 interpolation between the interval's ends. It runs only for xvals 9 outside [xtabs, xtabs1] (see the page). 10 11 """ 12 13 fs = (xvals - xtabs) / (xtabs1 - xtabs) 14 return ytabs + fs * (ytabs1 - ytabs) 15 16 17def spline(xvals: float, z1: float, z2: float, xtabs: float, ytabs: float, xtabs1: float, ytabs1: float, y1: float, 18 y2: float, e2: float, w2: float, v2: float, ncase: int) -> float: 19 """The value at xvals of the quadratic spline on [xtabs, xtabs1], made 20 of quadratic pieces joined at one knot z1 (ncase 1-3) or at two knots 21 y1 < z1 (ncase 4), with the parameters computed by MEVAL. 22 23 Translates: spline.f#SPLINE (7be4133) 24 25 Arguments: xvals; z1, z2 (the second knot and its value); xtabs, ytabs 26 and xtabs1, ytabs1 (the interval's ends); y1, y2 (the first knot and 27 its value, case 4); e2, w2, v2 (the middle control values); ncase. 28 Returns the spline's value. 29 30 """ 31 32 # Cases 1, 2 or 3: where xvals lies relative to the knot (Fortran 33 # lines 51–87). 34 if ncase != 4: 35 if z1 < xvals: # label 10 36 if xtabs1 != z1: 37 return (z2 * (xtabs1 - xvals) ** 2 38 + w2 * 2.0 * (xvals - z1) * (xtabs1 - xvals) 39 + ytabs1 * (xvals - z1) ** 2) / (xtabs1 - z1) ** 2 40 return _linear(xvals, xtabs, ytabs, xtabs1, ytabs1) 41 if z1 == xvals: # label 20 42 return z2 43 if z1 != xtabs: # label 30 44 return (ytabs * (z1 - xvals) ** 2 45 + v2 * 2.0 * (xvals - xtabs) * (z1 - xvals) 46 + z2 * (xvals - xtabs) ** 2) / (z1 - xtabs) ** 2 47 return _linear(xvals, xtabs, ytabs, xtabs1, ytabs1) 48 49 # Case 4 (label 40): where xvals lies relative to the first knot 50 # (Fortran lines 94–100); beyond it, label 70: where xvals lies 51 # relative to the second knot (Fortran lines 118–147). 52 if y1 < xvals: 53 if z1 < xvals: # label 100 54 if xtabs1 != z1: 55 return (z2 * (xtabs1 - xvals) ** 2 56 + w2 * 2.0 * (xvals - z1) * (xtabs1 - xvals) 57 + ytabs1 * (xvals - z1) ** 2) / (xtabs1 - z1) ** 2 58 return _linear(xvals, xtabs, ytabs, xtabs1, ytabs1) 59 if z1 == xvals: # label 90 60 return z2 61 if z1 != y1: # label 80 62 return (y2 * (z1 - xvals) ** 2 63 + e2 * 2.0 * (xvals - y1) * (z1 - xvals) 64 + z2 * (xvals - y1) ** 2) / (z1 - y1) ** 2 65 return _linear(xvals, xtabs, ytabs, xtabs1, ytabs1) 66 67 # Case 4, up to the first knot (Fortran lines 102–113). 68 if y1 == xvals: # label 60 69 return y2 70 if y1 != xtabs: # label 50 71 return (ytabs * (y1 - xvals) ** 2 72 + v2 * 2.0 * (xvals - xtabs) * (y1 - xvals) 73 + y2 * (xvals - xtabs) ** 2) / (y1 - xtabs) ** 2 74 return _linear(xvals, xtabs, ytabs, xtabs1, ytabs1)