splint.f

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

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

Variables on this page (1)

ishort

Purpose

SPLINT evaluates a cubic spline at one point $x$, from the table $(x_i, y_i)$ and the second derivatives $y''_i$ computed beforehand: by SPLINJ, or, for the opacity routine GTPURZ, by CSPLINE (computed there, or tabulated by ZSULAOL). The comments say it is taken from Numerical Recipes (Press et al., p. 89) and was changed from REAL*4, for the opacities, to REAL*8 (MHP, 10/02).

The header's "OUTPUT OF SPLINE ABOVE" (lines 7–10 ↗) is Numerical Recipes' text: Its spline (with end slopes) is CSPLINE here, while spline.f is an unrelated quadratic spline. The second derivatives come from SPLINJ in BOOLE and CALCAD (boole.f, lines 48–51 ↗; calcad.f, lines 347–348 ↗) and from CSPLINE in GTPURZ, directly (gtpurz.f, lines 115–116 ↗) and through the ZSDORL that ZSULAOL stores (gtpurz.f, lines 78–82 ↗; zsulaol.f, lines 40–44 ↗).

How it fits

Who calls it: BOOLE, CALCAD, and the low-temperature opacity routine GTPURZ (GTLAOL and GTLAOL2, and the table routine SULAOL, removed since rp_1, were the others).
Outputs: the value Y, and the bracketing indices KLO, KHI.
Modules: IOUNITS_MOD, for the unit ISHORT of its error message.

The method

A bisection finds $x_\mathrm{lo} \le x < x_\mathrm{hi}$ (KLO, KHI). With $h = x_\mathrm{hi} - x_\mathrm{lo}$, $A = (x_\mathrm{hi} - x)/h$ and $B = 1 - A$:

$$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].$$

If two abscissas coincide ($h = 0$) it writes an error, with the abscissa and $x$, to the terminal and to ISHORT, and stops the run (STOP 1). Outside the table it extrapolates the end interval's cubic, without a warning.

Design

KLO and KHI are arguments but are reset on entry, so a caller's previous bracket is not used as a starting guess.

Definitions

subroutineSPLINT

lines 11–56 · GitHub

Cubic-spline interpolation at one point (equation above), with arrays sized N; stops the run when two table abscissas coincide.

Uses modules IOUNITS_MOD

Called by BOOLE, CALCAD ×2, GTPURZ ×2

Compiler warnings (1)

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

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

Change log

The changes since rp_1

Changed since rp_1. Declarations only: /LUOUT/ to IOUNITS_MOD, IMPLICIT NONE; two callers, GTLAOL and GTLAOL2, were removed with their files.

Changed since a8a5950. The stop on two equal abscissas now gives the abscissa and $x$ in its message, writes it to the terminal as well as the .short file, and exits with status 1, the code for run-time errors, instead of a bare STOP (KNOWN_ISSUES.md, issue 38).

1C 2CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3C SUBROUTINE SPLINT 4CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 5C REAL*4 VERSION FOR OPACITIES CHANGED TO REAL*8 10/02 MHP 6C TAKEN FROM NUMERICAL RECIPES, PRESS, ET AL, P89. 7C GIVEN THE ARRAYS XA AND YA OF LENGTH N, WHICH TABULATE A FUNCTION 8C (WITH THE XA[I]'S IN ORDER), AND GIVEN THE ARRAY Y2A, WHICH IS THE 9C OUTPUT OF SPLINE ABOVE, AND GIVEN A VALUE OF X, THIS ROUTINE RETURNS 10C A CUBIC-SPLINE INTERPOLATED VALUE Y. 11 SUBROUTINE SPLINT(XA, YA, N, Y2A, X, Y, KLO, KHI) 12C 10/14/26 MHP via Claude 13C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 14 USE IOUNITS_MOD 15C 10/20/26 MHP via Claude 16C code changed to replace implicit typing with IMPLICIT NONE and 17C explicit declarations. 18 IMPLICIT NONE 19 REAL*8 A,B,H,X,XA,Y,Y2A,YA 20 INTEGER K,KHI,KLO,N 21C IMPLICIT REAL*4 (A-H,O-Z) 22C 9/25/26 MHP via Claude 23C code changed to consolidate all I/O logical unit numbers, previously 24C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 25C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 26C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 27C COMMON/IOUNITS/ 28 DIMENSION XA(N), YA(N), Y2A(N) 29 SAVE 30 KLO = 1 31 KHI = N 32 1 1 IF (KHI-KLO .GT. 1) THEN 33 K = (KHI+KLO)/2 34 IF (XA(K) .GT. X) THEN 35 KHI = K 36 ELSE 37 KLO = K 38 END IF 39 GOTO 1 40 END IF 41 H = XA(KHI) - XA(KLO) 42 IF (H .EQ. 0D0) THEN 43C 10/30/26 MHP via Claude: report the values, and to the terminal too. 44 WRITE(*,900)XA(KLO),X 45 WRITE(ISHORT,900)XA(KLO),X 46 900 900 FORMAT(' ERROR IN SPLINT ROUTINE: TWO ABSCISSAE COINCIDE AT', 47 * * 1PE14.6,' X',E14.6/' RUN STOPPED') 48C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 49 STOP 1 50 ENDIF 51 A = (XA(KHI)-X)/H 52 B = (X - XA(KLO))/H 53 Y = A*YA(KLO)+B*YA(KHI)+ 54 1 1 ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6.0D0 55 RETURN 56 END

1"""Python version of splint.f (YREC 7be4133): cubic-spline interpolation at 2one point (Numerical Recipes' SPLINT), with arrays sized N.""" 3 4from pyrec import iounits_mod 5import numpy as np 6 7 8def splint(xa: np.ndarray, ya: np.ndarray, n: int, y2a: np.ndarray, x: float, y: float, klo: int, 9 khi: int) -> tuple[float, int, int]: 10 """The cubic-spline value y at x, from the table (xa, ya) and its second 11 derivatives y2a (from SPLINJ or CSPLINE). Stops the run if two 12 abscissas coincide; outside the table it extrapolates the end interval. 13 14 Translates: splint.f#SPLINT (7be4133) 15 16 Arguments: xa[1..n] (increasing), ya[1..n]; n; y2a[1..n]; x; y, klo and 17 khi, set only (klo and khi are reset on entry). Returns (y, klo, khi), 18 with klo, khi the bracketing indices. 19 20 """ 21 22 # Bisection for the bracketing interval (Fortran lines 30–40). 23 klo = 1 24 khi = n 25 while khi - klo > 1: # label 1 26 k = (khi + klo) // 2 27 if xa[k] > x: 28 khi = k 29 else: 30 klo = k 31 32 # Stop if two abscissas coincide (Fortran lines 41–50). 33 h = xa[khi] - xa[klo] 34 if h == 0.0: 35 text = (f" ERROR IN SPLINT ROUTINE: TWO ABSCISSAE COINCIDE AT" 36 f"{xa[klo]:14.6E} X{x:14.6E}\n RUN STOPPED") # format 900 37 print(text) 38 iounits_mod.write(iounits_mod.ishort, text) 39 raise SystemExit(1) # STOP 1, the run-time error code 40 41 # The cubic spline's value (Fortran lines 51–54). 42 a = (xa[khi] - x) / h 43 b = (x - xa[klo]) / h 44 y = (a * ya[klo] + b * ya[khi] 45 + ((a ** 3 - a) * y2a[klo] + (b ** 3 - b) * y2a[khi]) * (h ** 2) / 6.0) 46 return y, klo, khi