polint.f

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

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

Variables on this page (1)

ishort

Purpose

Polynomial interpolation through N points by Neville's algorithm (the routine POLINT of Numerical Recipes, by its name and form): it returns the value Y at X and an error estimate DY, the last correction added.

How it fits

Who calls it: ALSURFP (the Allard model atmospheres' surface pressures) and FPFT (the rotational distortion factors).
Modules: IOUNITS_MOD, for the unit ISHORT of its error message.
Limits: at most 20 points (the arrays are fixed at 20, "for consistency with fpft," MHP 10/02).

The method

Starting from the nearest tabulated point, the tableau of differences $C_{i}$, $D_{i}$ is built column by column, $$D_i \leftarrow \frac{(x_{i+j}-x)\,(C_{i+1}-D_i)}{x_i - x_{i+j}},\qquad C_i \leftarrow \frac{(x_i-x)\,(C_{i+1}-D_i)}{x_i - x_{i+j}},$$ and the path through it that stays closest to the centre adds one correction per column to $y$.

Two equal abscissas ($|x_i - x_{i+j}| < 10^{-20}$) stop the run (STOP 1) with a message to the terminal and to ISHORT giving the two abscissas, $x$ and N.

Definitions

subroutinePOLINT

lines 6–61 · GitHub

Uses modules IOUNITS_MOD

Called by ALSURFP ×6, FPFT

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE.

Changed since a8a5950. The stop on two equal abscissas now writes a message with the values to the terminal and to the .short file (the routine now uses IOUNITS_MOD), and exits with status 1, the code for run-time errors, instead of a bare STOP (KNOWN_ISSUES.md, issue 38). This settles former question 1 (a stop without a message), which is removed.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C POLINT 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE POLINT(XA,YA,N,X,Y,DY) 7C 10/20/26 MHP via Claude 8C code changed to replace implicit typing with IMPLICIT NONE and 9C explicit declarations. 10C 10/30/26 MHP via Claude 11C code changed to USE IOUNITS_MOD for ISHORT, so the error stop below 12C also reports to the .short file (issue 38). 13 USE IOUNITS_MOD 14 IMPLICIT NONE 15 REAL*8 C,D,DEN,DIF,DIFT,DY,HO,HP,W,X,XA,Y,YA 16 INTEGER I,J,N,NS 17C MHP 10/02 dimensions changed for consistency with fpft 18 DIMENSION XA(20),YA(20),C(20),D(20) 19 SAVE 20 21 NS = 1 22 DIF = DABS(X-XA(1)) 23 DO 10 I = 1,N 24 DIFT = DABS(X-XA(I)) 25 IF(DIFT.LT.DIF) THEN 26 NS = I 27 DIF = DIFT 28 ENDIF 29 C(I) = YA(I) 30 D(I) = YA(I) 31 10 10 CONTINUE 32 Y = YA(NS) 33 NS = NS - 1 34 DO 30 J = 1,N-1 35 DO 20 I = 1,N-J 36 HO = XA(I)-X 37 HP = XA(I+J)-X 38 W = C(I+1) - D(I) 39 DEN = HO - HP 40 IF(dabs(DEN).LT.1.0D-20)THEN 41 WRITE(*,900)XA(I),XA(I+J),X,N 42 WRITE(ISHORT,900)XA(I),XA(I+J),X,N 43 900 900 FORMAT(' POLINT: TWO ABSCISSAE COINCIDE',1P3E14.6, 44 * * ' N',I3/' RUN STOPPED') 45C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 46 STOP 1 47 ENDIF 48 DEN = W/DEN 49 D(I) = HP*DEN 50 C(I) = HO*DEN 51 20 20 CONTINUE 52 IF(2*NS.LT.N-J) THEN 53 DY = C(NS+1) 54 ELSE 55 DY = D(NS) 56 NS = NS-1 57 ENDIF 58 Y = Y + DY 59 30 30 CONTINUE 60 RETURN 61 END

1"""Python version of polint.f (YREC 7be4133): polynomial interpolation by 2Neville's algorithm (Numerical Recipes' POLINT).""" 3 4import numpy as np 5 6from pyrec import iounits_mod 7 8# Local work arrays (SAVE), fixed at 20 points. 9c = np.zeros(21) 10d = np.zeros(21) 11 12 13def polint(xa: np.ndarray, ya: np.ndarray, n: int, x: float, y: float, dy: float) -> tuple[float, float]: 14 """Interpolate the polynomial through the n points (xa, ya) at x, 15 starting from the nearest tabulated point. 16 17 Translates: polint.f#POLINT (7be4133) 18 19 Arguments: xa[1..n], ya[1..n] (n <= 20); x; y, dy (set). Returns 20 (y, dy): the value at x and the last correction added, an error 21 estimate. 22 23 """ 24 25 # The nearest tabulated point, and the tableau's first column 26 # (Fortran lines 21–33). 27 ns = 1 28 dif = abs(x - xa[1]) 29 for i in range(1, n + 1): 30 dift = abs(x - xa[i]) 31 if dift < dif: 32 ns = i 33 dif = dift 34 c[i] = ya[i] 35 d[i] = ya[i] 36 y = ya[ns] 37 ns = ns - 1 38 39 # Neville's tableau, one column per pass, adding one correction to y 40 # each (Fortran lines 34–59). 41 for j in range(1, n): 42 for i in range(1, n - j + 1): 43 ho = xa[i] - x 44 hp = xa[i + j] - x 45 w = c[i + 1] - d[i] 46 den = ho - hp 47 if abs(den) < 1.0e-20: 48 text = (f" POLINT: TWO ABSCISSAE COINCIDE{xa[i]:14.6E}" 49 f"{xa[i + j]:14.6E}{x:14.6E} N{n:3d}\n" 50 " RUN STOPPED") # format 900 51 print(text) 52 iounits_mod.write(iounits_mod.ishort, text) 53 raise SystemExit(1) # STOP 1, the run-time error code 54 den = w / den 55 d[i] = hp * den 56 c[i] = ho * den 57 if 2 * ns < n - j: 58 dy = c[ns + 1] 59 else: 60 dy = d[ns] 61 ns = ns - 1 62 y = y + dy 63 return y, dy