cspline.f

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

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

Variables on this page (1)

json

Purpose

The second derivatives of a cubic spline through tabulated points, from Numerical Recipes (SPLINE), in double precision "for opacities."

How it fits

Inputs: X(N) (increasing), Y(N), and the end slopes YP1, YPN; a slope of $10^{30}$ or more asks for a natural spline (zero second derivative) at that end.
Who calls it: the opacity routines GTPURZ and ZSULAOL; at rp_1 also GTLAOL, GTLAOL2 and SULAOL (removed since rp_1).
Output: Y2(N), the second derivatives, for a spline evaluation routine.
Limits: the work array U has the fixed size JSON (5000, from DIMENSIONS_MOD; line 27 ↗), and N is not checked against it; the callers pass at most 104 points (the pure-Z table's densities and temperatures, ZLAOL_MOD).

The method

Continuity of the first derivative at interior points gives the tridiagonal system $$\frac{h_{i-1}}{6}y''_{i-1} + \frac{h_{i-1} + h_i}{3}y''_i + \frac{h_i}{6}y''_{i+1} = \frac{y_{i+1} - y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}},$$ with $h_i = x_{i+1} - x_i$, solved by forward elimination and back substitution.

Definitions

subroutineCSPLINE

lines 16–55 · GitHub

Numerical Recipes' cubic-spline setup: returns the second derivatives Y2, natural or with given end slopes.

Uses modules DIMENSIONS_MOD

Called by GTPURZ, ZSULAOL

Change log

The changes since rp_1

Changed since rp_1. Declarations only: JSON from DIMENSIONS_MOD, IMPLICIT NONE; three of its callers removed with their files.

1C 2CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3C SUBROUTINE CSPLINE 4CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 5C REAL*8 VERSION FOR OPACITIES 6C TAKEN FROM NUMERICAL RECIPES, PRESS ET AL, P88 7C GIVEN ARRAYS X AND Y OF LENGTH N CONTAINING A TABULATED FUNCITION, 8C I.E. Y[I] = F(X[I]) WITH X[1] < X[2] < ... X[N], AND GIVEN VALUES 9C YP1 AND YPN FOR THE FIRST DERIVATIVE OF THE INTERPOLATING FUNCTION 10C AT POINTS 1 AND N, RESPECTIVELY, THIS ROUTINE RETURNS AND ARRAY Y2 11C OF LENGTH N WHICH CONTAINS THE SECOND DERIVATIVES OF THE INTERPOLATING 12C FUNCTION AT THE TABULATED POINTS X[I]. IF YP1 AND/OR YPN ARE EQUAL 13C TO 1.0E30 OR LARGER, THE ROUTINE IS SIGNALLED TO SET THE CORRESPONDING 14C BOUNDARY CONDITION FOR A NATURAL SPLINE, WITH ZERO SECOND DERIVATIVE 15C ON THAT BOUNDARY. 16 SUBROUTINE CSPLINE(X, Y, N, YP1, YPN, Y2) 17C 10/15/26 MHP via Claude 18C code changed to take JSON from DIMENSIONS_MOD instead of declaring 19C PARAMETER(JSON=5000) here. 20 USE DIMENSIONS_MOD 21C 10/20/26 MHP via Claude 22C code changed to replace implicit typing with IMPLICIT NONE and 23C explicit declarations. 24 IMPLICIT NONE 25 REAL*8 P,QN,SIG,U,UN,X,Y,Y2,YP1,YPN 26 INTEGER I,K,N 27 DIMENSION X(N), Y(N), Y2(N), U(JSON) 28 SAVE 29 IF (YP1 .GT. 0.99D30) THEN 30 Y2(1) = 0.0D0 31 U(1) = 0.0D0 32 ELSE 33 Y2(1) = -0.5D0 34 U(1) = (3.D0/(X(2)-X(1)))*((Y(2)-Y(1))/(X(2)-X(1))-YP1) 35 ENDIF 36 DO I=2, N-1 37 SIG = (X(I)-X(I-1))/(X(I+1)-X(I-1)) 38 P = SIG*Y2(I-1)+2.0D0 39 Y2(I) = (SIG-1.0D0)/P 40 U(I) = (6.0D0*((Y(I+1)-Y(I))/(X(I+1)-X(I))-(Y(I)-Y(I-1)) 41 1 1 /(X(I)-X(I-1)))/(X(I+1)-X(I-1))-SIG*U(I-1))/P 42 END DO 43 IF (YPN .GT. 0.99D30) THEN 44 QN = 0.0D0 45 UN = 0.0D0 46 ELSE 47 QN = 0.5D0 48 UN = (3.D0/(X(N)-X(N-1)))*(YPN-(Y(N)-Y(N-1))/(X(N)-X(N-1))) 49 END IF 50 Y2(N) = (UN-QN*U(N-1))/(QN*Y2(N-1)+1.0D0) 51 DO K=N-1, 1, -1 52 Y2(K) = Y2(K)*Y2(K+1)+U(K) 53 ENDDO 54 RETURN 55 END

1"""Python version of cspline.f (YREC 7be4133): the second derivatives of a 2cubic spline (Numerical Recipes' SPLINE), for the opacity tables.""" 3 4import numpy as np 5 6from pyrec import dimensions_mod 7 8JSON = dimensions_mod.json 9 10# The work array U (SAVE), of fixed size JSON; N is not checked against it. 11u = np.zeros(JSON + 1) 12 13 14def cspline(x: np.ndarray, y: np.ndarray, n: int, yp1: float, ypn: float, y2: np.ndarray) -> None: 15 """The second derivatives y2 of the cubic spline through (x, y), with 16 end slopes yp1, ypn; a slope of 0.99e30 or more gives a natural 17 spline (zero second derivative) at that end. 18 19 Translates: cspline.f#CSPLINE (7be4133) 20 21 Arguments: x[1..n] (increasing), y[1..n]; n; yp1, ypn; y2[1..n], 22 filled in place. None is returned. 23 24 """ 25 26 # The lower boundary condition (Fortran lines 29–35). 27 if yp1 > 0.99e30: 28 y2[1] = 0.0 29 u[1] = 0.0 30 else: 31 y2[1] = -0.5 32 u[1] = (3.0 / (x[2] - x[1])) * ((y[2] - y[1]) / (x[2] - x[1]) - yp1) 33 34 # Forward elimination (Fortran lines 36–42). 35 for i in range(2, n): 36 sig = (x[i] - x[i - 1]) / (x[i + 1] - x[i - 1]) 37 p = sig * y2[i - 1] + 2.0 38 y2[i] = (sig - 1.0) / p 39 u[i] = (6.0 * ((y[i + 1] - y[i]) / (x[i + 1] - x[i]) 40 - (y[i] - y[i - 1]) / (x[i] - x[i - 1])) 41 / (x[i + 1] - x[i - 1]) - sig * u[i - 1]) / p 42 43 # The upper boundary condition (Fortran lines 43–50). 44 if ypn > 0.99e30: 45 qn = 0.0 46 un = 0.0 47 else: 48 qn = 0.5 49 un = (3.0 / (x[n] - x[n - 1])) * (ypn - (y[n] - y[n - 1]) 50 / (x[n] - x[n - 1])) 51 y2[n] = (un - qn * u[n - 1]) / (qn * y2[n - 1] + 1.0) 52 53 # Back substitution (Fortran lines 51–53). 54 for k in range(n - 1, 0, -1): 55 y2[k] = y2[k] * y2[k + 1] + u[k]