splinj.f

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

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

Purpose

SPLINJ computes the second derivatives of a natural cubic spline through $(x_i, y_i)$, $i = 1, \dots, N$ ($y''_1 = y''_N = 0$), for interpolation with SPLINT. The algorithm, with its equations, is that of splinc.f; only the arrays differ: here they are sized N (adjustable arrays), not 5000, and the routine has no SAVE. The work array U(N) (line 10 ↗) is a local with an adjustable size, an automatic array: not allowed in Fortran 77, standard in Fortran 90, which YREC is (its modules, IMPLICIT NONE and the F90 spellings already in use: docs/archive/F90_CONVERSION_PLAN.md, line 38) (settled from the code, Nov 9, 2026).

How it fits

Who calls it: BOOLE and CALCAD, which then call SPLINT.
COMMON blocks: none.

Definitions

subroutineSPLINJ

lines 2–29 · GitHub

Natural cubic spline's second derivatives, with arrays sized N; the same algorithm as SPLINC.

Called by BOOLE, CALCAD ×2

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE.

1C$$$$$$ 2 SUBROUTINE SPLINJ(X,Y,Y2,N) 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 REAL*8 P,QN,SIG,UN 8 INTEGER I,K,N 9C PARAMETER (JSON=5000) 10 REAL*8 X(N),Y(N),Y2(N),U(N) 11C SAVE 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 splinj.f (YREC 7be4133): the second derivatives of a 2natural cubic spline, with arrays sized N (the algorithm of SPLINC).""" 3 4import numpy as np 5 6 7def splinj(x: np.ndarray, y: np.ndarray, y2: np.ndarray, n: int) -> None: 8 """The second derivatives y2 of the natural cubic spline through 9 (x, y) (y2 = 0 at both ends), for interpolation with SPLINT. 10 11 Translates: splinj.f#SPLINJ (7be4133) 12 13 Arguments: x[1..n] (increasing), y[1..n]; y2[1..n], filled in place; 14 n. None is returned. 15 16 """ 17 18 # U is an automatic array of size N (Fortran line 10). 19 u = np.zeros(n + 1) 20 21 # The natural spline's lower end (Fortran lines 13–14). 22 y2[1] = 0.0 23 u[1] = 0.0 24 25 # Forward elimination (Fortran lines 15–21). 26 for i in range(2, n): 27 sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]) 28 p = sig * y2[i - 1] + 2.0 29 y2[i] = (sig - 1.0) / p 30 u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) 31 - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) 32 / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p 33 34 # The natural spline's upper end (Fortran lines 22–24). 35 qn = 0.0 36 un = 0.0 37 y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0) 38 39 # Back substitution (Fortran lines 25–27). 40 for k in range(n - 1, 0, -1): 41 y2[k] = y2[k] * y2[k + 1] + u[k]