splintd2.f

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

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

Variables on this page (2)

ishort, json

Purpose

SPLINTD2 evaluates a cubic spline at one point; it is splint.f's algorithm and equation, with arrays of the fixed size JSON (5000, from DIMENSIONS_MOD), for the second derivatives made by SPLINC. Its header comment still names it "SUBROUTINE SPLINTD" (line 3 ↗).

How it fits

Who calls it: only HPOINT, after SPLINC.
Modules: IOUNITS_MOD, for the unit ISHORT of its error message ("ERROR IN SPLINT ROUTINE: TWO ABSCISSAE COINCIDE," with the abscissa and $x$, also written to the terminal, then STOP 1).

Definitions

subroutineSPLINTD2

lines 11–59 · GitHub

Cubic-spline interpolation at one point, arrays sized 5000; the same as SPLINT.

Uses modules DIMENSIONS_MOD, IOUNITS_MOD

Called by HPOINT ×3

Compiler warnings (1)

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

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

Change log

The changes since rp_1

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

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). SPLINT changed the same way, so question 1 stands.

1C 2CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 3C SUBROUTINE SPLINTD 4CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 5C REAL*8 VERSION 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 SPLINTD2(XA, YA, N, Y2A, X, Y, KLO, KHI) 12C 10/15/26 MHP via Claude 13C code changed to take JSON from DIMENSIONS_MOD instead of declaring 14C PARAMETER(JSON=5000) here. 15 USE DIMENSIONS_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 18 USE IOUNITS_MOD 19C 10/20/26 MHP via Claude 20C code changed to replace implicit typing with IMPLICIT NONE and 21C explicit declarations. 22 IMPLICIT NONE 23 REAL*8 A,B,H,X,XA,Y,Y2A,YA 24 INTEGER K,KHI,KLO,N 25C 9/25/26 MHP via Claude 26C code changed to consolidate all I/O logical unit numbers, previously 27C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 28C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 29C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 30C COMMON/IOUNITS/ 31 DIMENSION XA(JSON), YA(JSON), Y2A(JSON) 32 SAVE 33 KLO = 1 34 KHI = N 35 1 1 IF (KHI-KLO .GT. 1) THEN 36 K = (KHI+KLO)/2 37 IF (XA(K) .GT. X) THEN 38 KHI = K 39 ELSE 40 KLO = K 41 END IF 42 GOTO 1 43 END IF 44 H = XA(KHI) - XA(KLO) 45 IF (H .EQ. 0D0) THEN 46C 10/30/26 MHP via Claude: report the values, and to the terminal too. 47 WRITE(*,900)XA(KLO),X 48 WRITE(ISHORT,900)XA(KLO),X 49 900 900 FORMAT(' ERROR IN SPLINT ROUTINE: TWO ABSCISSAE COINCIDE AT', 50 * * 1PE14.6,' X',E14.6/' RUN STOPPED') 51C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 52 STOP 1 53 ENDIF 54 A = (XA(KHI)-X)/H 55 B = (X - XA(KLO))/H 56 Y = A*YA(KLO)+B*YA(KHI)+ 57 1 1 ((A**3-A)*Y2A(KLO)+(B**3-B)*Y2A(KHI))*(H**2)/6D0 58 RETURN 59 END

1"""Python version of splintd2.f (YREC 7be4133): cubic-spline interpolation at 2one point, with arrays of the fixed size JSON (SPLINT's algorithm).""" 3 4from pyrec import iounits_mod 5import numpy as np 6 7 8def splintd2(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 SPLINC). Stops the run if two abscissas coincide; 12 outside the table it extrapolates the end interval. 13 14 Translates: splintd2.f#SPLINTD2 (7be4133) 15 16 Arguments: xa[1..n] (increasing), ya[1..n]; n (at most JSON); 17 y2a[1..n]; x; y, klo and khi, set only (klo and khi are reset on 18 entry). Returns (y, klo, khi), with klo, khi the bracketing indices. 19 20 """ 21 22 # Bisection for the bracketing interval (Fortran lines 33–43). 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 44–53). 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 54–57). 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