kspline.f

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

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

Purpose

KSPLINE computes the second derivatives of a natural cubic spline through exactly 4 points (the size NM = 4 is fixed), as in Numerical Recipes' SPLINE. KSPLINT evaluates the spline.

How it fits

Called by the surface-pressure routines KCSURFP and SURFP, the SCV equation of state EQSCVE, and TAUINTNEW. No COMMON blocks, no calls.

The method

The tridiagonal system for the second derivatives $y''_i$ of a natural spline ($y''_1 = y''_4 = 0$), solved by forward elimination and back-substitution.

Definitions

subroutineKSPLINE

lines 2–29 · GitHub

Natural cubic spline on 4 points: returns $y''$ in Y2.

Called by EQSCVE ×8, KCSURFP ×4, SURFP ×4, TAUINTNEW ×4

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE with explicit declarations.

1C$$$$$$ 2 SUBROUTINE KSPLINE(X,Y,Y2) 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 I,K,N,NM 8 PARAMETER (NM=4) 9 REAL*8 X(NM),Y(NM),Y2(NM),U(NM),SIG,QN,UN,P 10 SAVE 11 N = NM 12C NATURAL SPLINE 13 Y2(1) = 0.0D0 14 U(1) = 0.0D0 15 DO I = 2,N-1 16 SIG = (X(I)-X(I-1))/(X(I+1)-X(I-1)) 17 P = SIG*Y2(I-1)+2.0D0 18 Y2(I) = (SIG-1.0D0)/P 19 U(I) = (6.0D0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) 20 * * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 21 END DO 22 QN = 0.0D0 23 UN = 0.0D0 24 Y2(N) = (UN-QN*U(N-1))/(QN*Y2(N-1)+1.0D0) 25 DO K = N-1,1,-1 26 Y2(K) = Y2(K)*Y2(K+1)+U(K) 27 END DO 28 RETURN 29 END

1"""Python version of kspline.f (YREC development, 7be4133).""" 2 3import numpy as np 4 5NM = 4 6 7 8def kspline(x: np.ndarray, y: np.ndarray, y2: np.ndarray) -> None: 9 """Second derivatives of a natural cubic spline through 4 points. 10 11 *Numerical Recipes*' SPLINE with the size fixed at ``NM`` = 4 and 12 natural end conditions (``y2[1] = y2[4] = 0``). KSPLINT evaluates it. 13 14 Translates: kspline.f#KSPLINE (7be4133) 15 16 Parameters 17 ---------- 18 x, y : ndarray, shape (NM + 1,) 19 The nodes and values, ``x[1..4]``, ``y[1..4]``. 20 y2 : ndarray, shape (NM + 1,) 21 Set in place: the second derivatives at the nodes. 22 23 """ 24 25 # The work array U is SAVEd in the Fortran but fully reset each call. 26 u = np.zeros(NM + 1) 27 n = NM 28 29 # The natural spline's lower end (Fortran lines 11–14). 30 y2[1] = 0.0 31 u[1] = 0.0 32 33 # Forward elimination (Fortran lines 15–21). 34 for i in range(2, n): 35 sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]) 36 p = sig * y2[i - 1] + 2.0 37 y2[i] = (sig - 1.0) / p 38 u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) 39 - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) 40 / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p 41 42 # The natural spline's upper end (Fortran lines 22–24). 43 qn = 0.0 44 un = 0.0 45 y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0) 46 47 # Back substitution (Fortran lines 25–27). 48 for k in range(n - 1, 0, -1): 49 y2[k] = y2[k] * y2[k + 1] + u[k]