splinc.f

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

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

Variables on this page (1)

json

Purpose

SPLINC computes the second derivatives of a natural cubic spline through the points $(x_i, y_i)$, $i = 1, \dots, N$, for interpolation with SPLINTD2. It is the tridiagonal algorithm of Numerical Recipes' spline, with the end conditions fixed to a natural spline.

How it fits

Who calls it: only HPOINT, which then calls SPLINTD2 with the result.
Outputs: Y2, the second derivatives at the points.
COMMON blocks: none.

The method

A natural spline has $y''_1 = y''_N = 0$. For the interior points the second derivatives solve

$$\frac{x_i - x_{i-1}}{6}\,y''_{i-1} + \frac{x_{i+1} - x_{i-1}}{3}\,y''_i + \frac{x_{i+1} - x_i}{6}\,y''_{i+1} = \frac{y_{i+1} - y_i}{x_{i+1} - x_i} - \frac{y_i - y_{i-1}}{x_i - x_{i-1}},$$

which the code solves by forward elimination (SIG, P, U) and back substitution.

Design

The arrays have the fixed size JSON (5000, from DIMENSIONS_MOD). splinj.f is the same routine with arrays sized N.

Definitions

subroutineSPLINC

lines 2–32 · GitHub

Natural cubic spline: the second derivatives Y2 at N points (at most 5000).

Uses modules DIMENSIONS_MOD

Called by HPOINT ×3

Change log

The changes since rp_1

Changed since rp_1. Declarations only: the array size JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C$$$$$$ 2 SUBROUTINE SPLINC(X,Y,Y2,N) 3C 10/15/26 MHP via Claude 4C code changed to take JSON from DIMENSIONS_MOD instead of declaring 5C PARAMETER(JSON=5000) here. 6 USE DIMENSIONS_MOD 7C 10/20/26 MHP via Claude 8C code changed to replace implicit typing with IMPLICIT NONE and 9C explicit declarations. 10 IMPLICIT NONE 11 REAL*8 P,QN,SIG,UN 12 INTEGER I,K,N 13 REAL*8 X(JSON),Y(JSON),Y2(JSON),U(JSON) 14 SAVE 15C NATURAL SPLINE 16 Y2(1) = 0.0D0 17 U(1) = 0.0D0 18 DO I = 2,N-1 19 SIG = (X(I)-X(I-1))/(X(I+1)-X(I-1)) 20 P = SIG*Y2(I-1)+2.0D0 21 Y2(I) = (SIG-1.0D0)/P 22 U(I) = (6.0D0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) 23 * * /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 24 END DO 25 QN = 0.0D0 26 UN = 0.0D0 27 Y2(N) = (UN-QN*U(N-1))/(QN*Y2(N-1)+1.0D0) 28 DO K = N-1,1,-1 29 Y2(K) = Y2(K)*Y2(K+1)+U(K) 30 END DO 31 RETURN 32 END

1"""Python version of splinc.f (YREC 7be4133): the second derivatives of a 2natural cubic spline, with arrays of the fixed size JSON.""" 3 4import numpy as np 5 6from pyrec import dimensions_mod 7 8JSON = dimensions_mod.json 9 10# The work array U (SAVE). 11u = np.zeros(JSON + 1) 12 13 14def splinc(x: np.ndarray, y: np.ndarray, y2: np.ndarray, n: int) -> None: 15 """The second derivatives y2 of the natural cubic spline through 16 (x, y) (y2 = 0 at both ends), for interpolation with SPLINTD2. 17 18 Translates: splinc.f#SPLINC (7be4133) 19 20 Arguments: x[1..n] (increasing), y[1..n]; y2[1..n], filled in place; 21 n (at most JSON). None is returned. 22 23 """ 24 25 # The natural spline's lower end (Fortran lines 16–17). 26 y2[1] = 0.0 27 u[1] = 0.0 28 29 # Forward elimination (Fortran lines 18–24). 30 for i in range(2, n): 31 sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]) 32 p = sig * y2[i - 1] + 2.0 33 y2[i] = (sig - 1.0) / p 34 u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) 35 - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) 36 / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p 37 38 # The natural spline's upper end (Fortran lines 25–27). 39 qn = 0.0 40 un = 0.0 41 y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0) 42 43 # Back substitution (Fortran lines 28–30). 44 for k in range(n - 1, 0, -1): 45 y2[k] = y2[k] * y2[k + 1] + u[k]