ksplint.f

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

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

Variables on this page (1)

ishort

Purpose

KSPLINT evaluates at $x$ the 4-point cubic spline whose second derivatives KSPLINE computed (Numerical Recipes' SPLINT, with NM = 4 fixed).

How it fits

Called by the same routines as KSPLINE: KCSURFP, SURFP, EQSCVE, TAUINTNEW. No calls; it uses IOUNITS_MOD for the unit ISHORT of its error message.

The method

It finds the interval $[x_\mathrm{lo}, x_\mathrm{hi}]$ by bisection, then with $h = x_\mathrm{hi}-x_\mathrm{lo}$, $A = (x_\mathrm{hi}-x)/h$, $B = (x-x_\mathrm{lo})/h$:

$$y = A\,y_\mathrm{lo} + B\,y_\mathrm{hi} + \frac{h^2}{6}\left[(A^3-A)\,y''_\mathrm{lo} + (B^3-B)\,y''_\mathrm{hi}\right].$$

Outside the 4 points it extrapolates with the end interval's cubic. If two nodes coincide it stops the run (STOP 1), with a message giving the node and $x$, written to the terminal and to ISHORT.

Definitions

subroutineKSPLINT

lines 2–43 · GitHub

Spline evaluation (above).

Uses modules IOUNITS_MOD

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

Compiler warnings (1)

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

line 29: Equality comparison for REAL(8) [-Wcompare-reals]

Change log

The changes since rp_1

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

Changed since a8a5950. The stop on two equal nodes now writes a message with the node and $x$ (it was "Ksplint failure," to the terminal only) to the terminal and to the .short file (the routine now uses IOUNITS_MOD), and exits with status 1, the code for run-time errors, instead of a bare STOP (KNOWN_ISSUES.md, issue 38).

1C$$$$$$ 2 SUBROUTINE KSPLINT(XA,YA,Y2A,X,Y) 3C 10/20/26 MHP via Claude 4C code changed to replace implicit typing with IMPLICIT NONE and 5C explicit declarations. 6C 10/30/26 MHP via Claude 7C code changed to USE IOUNITS_MOD for ISHORT, so the error stop below 8C also reports to the .short file (issue 38). 9 USE IOUNITS_MOD 10 IMPLICIT NONE 11 INTEGER K,KHI,KLO,NM 12 PARAMETER(NM=4) 13 REAL*8 XA(NM),YA(NM),Y2A(NM),X,Y,H,A,B 14 SAVE 15 KLO = 1 16 KHI = NM 17 1 1 CONTINUE 18 IF(KHI-KLO.GT.1)THEN 19 K = (KHI+KLO)/2 20 IF(XA(K).GT.X)THEN 21 KHI = K 22 ELSE 23 KLO = K 24 ENDIF 25 GOTO 1 26 ENDIF 27c WRITE(*,*)KHI,KLO,XA(KHI),XA(KLO),X 28 H = XA(KHI) - XA(KLO) 29 IF(H.EQ.0D0) THEN 30 WRITE(*,900)XA(KLO),X 31 WRITE(ISHORT,900)XA(KLO),X 32 900 900 FORMAT(' KSPLINT: TWO ABSCISSAE COINCIDE AT',1PE14.6, 33 * * ' X',E14.6/' RUN STOPPED') 34C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 35 STOP 1 36 ENDIF 37C IF(H.EQ.0D0) STOP 911 38 A = (XA(KHI)-X)/H 39 B = (X-XA(KLO))/H 40 Y = A*YA(KLO)+B*YA(KHI)+ 41 * * ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.0D0 42 RETURN 43 END

1"""Python version of ksplint.f (YREC development, 7be4133).""" 2 3from pyrec import iounits_mod 4import numpy as np 5 6NM = 4 7 8 9def ksplint(xa: np.ndarray, ya: np.ndarray, y2a: np.ndarray, x: float, y: float) -> float: 10 """Evaluate at ``x`` the 4-point cubic spline from KSPLINE. 11 12 *Numerical Recipes*' SPLINT with ``NM`` = 4: bisection for the 13 interval, then the cubic-spline formula. Outside the nodes it 14 extrapolates with the end interval's cubic. 15 16 Translates: ksplint.f#KSPLINT (7be4133) 17 18 Parameters 19 ---------- 20 xa, ya, y2a : ndarray, shape (NM + 1,) 21 Nodes, values and second derivatives, indices 1..4. 22 x : float 23 The point of evaluation. 24 y : float 25 Output only (passed in with any value). 26 27 Returns 28 ------- 29 y : float 30 The spline's value at ``x``. 31 32 """ 33 34 # Bisection for the bracketing interval (Fortran lines 15–26). 35 klo = 1 36 khi = NM 37 while khi - klo > 1: # label 1 38 k = (khi + klo) // 2 39 if xa[k] > x: 40 khi = k 41 else: 42 klo = k 43 44 # Stop if two nodes coincide (Fortran lines 28–36). 45 h = xa[khi] - xa[klo] 46 if h == 0.0: 47 text = (f" KSPLINT: TWO ABSCISSAE COINCIDE AT{xa[klo]:14.6E}" 48 f" X{x:14.6E}\n RUN STOPPED") # format 900 49 print(text) 50 iounits_mod.write(iounits_mod.ishort, text) 51 raise SystemExit(1) # STOP 1, the run-time error code 52 53 # The cubic spline's value (Fortran lines 38–41). 54 a = (xa[khi] - x) / h 55 b = (x - xa[klo]) / h 56 y = (a * ya[klo] + b * ya[khi] 57 + ((a ** 3 - a) * y2a[klo] + (b ** 3 - b) * y2a[khi]) * (h ** 2) / 6.0) 58 return y