interp.f

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

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

Purpose

INTERP returns the weights of 4-point Lagrange interpolation, and of its first derivative, at a point $x$.

How it fits

Called by EQSCVE and EQSCVG (the SCV equation of state), GETALEX06 and GETOPAL95 (opacity tables; YALO3D is removed since rp_1), and MICRODIFF_MTE and MODEL_TO_EQUAL (regridding for diffusion). No modules, no calls. The version without derivatives is INTRP2; the 3-point one is INTER3.

The header's "VandenBerg opacity tables" is from an older use. The comment calls it "THE INTERPOLATION ROUTINE FOR THE VANDENBERG OPACITY TABLES(CAPPA)" (lines 16–18 ↗); it now serves the SCV EOS (EQSCVE, EQSCVG), the opacity tables of GETOPAL95 and GETALEX06, and the diffusion grids (MICRODIFF_MTE, MODEL_TO_EQUAL), so the comment could describe a general four-point Lagrange interpolation.

The method

For nodes $a_1,\dots,a_4$,

$$B_i(x) = \prod_{j\ne i}\frac{x-a_j}{a_i-a_j},\qquad C_i(x) = \frac{dB_i}{dx},$$

returned in B and C; then $f(x)\approx\sum_i B_i f(a_i)$ and $f'(x)\approx\sum_i C_i f(a_i)$. Used in a table, the interpolant is piecewise: it is continuous where the four-point stencil shifts from one cell to the next, but its derivative, and so the derivative weights C the equation of state and opacities use, jump there.

Definitions

subroutineINTERP

lines 6–43 · GitHub

4-point Lagrange weights and derivative weights (above).

Called by EQSCVE ×2, EQSCVG ×2, GETALEX06 ×2, GETOPAL95 ×5, MICRODIFF_MTE ×2, MODEL_TO_EQUAL ×2

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE. One caller, YALO3D, is removed with the obsolete opacity options.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C INTERP 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE INTERP(A,B,C,X) 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 A,A21,A31,A32,A41,A42,A43,B,C,D1,D2,D3,D4,X,XA1,XA2,XA3,XA4 12 13 DIMENSION A(4),B(4),C(4) 14 SAVE 15 16C INTERP IS THE INTERPOLATION ROUTINE FOR THE VANDENBERG 17C OPACITY TABLES(CAPPA), AND IT USES A 4-POINT LAGRANGIAN 18C INTERPOLATION SCHEME. 19 A43 = A(4) - A(3) 20 A42 = A(4) - A(2) 21 A41 = A(4) - A(1) 22 A32 = A(3) - A(2) 23 A31 = A(3) - A(1) 24 A21 = A(2) - A(1) 25 D1 = -A21*A31*A41 26 D2 = A21*A32*A42 27 D3 = -A31*A32*A43 28 D4 = A41*A42*A43 29 XA1 = X - A(1) 30 XA2 = X - A(2) 31 XA3 = X - A(3) 32 XA4 = X - A(4) 33 B(1) = (XA2*XA3*XA4)/D1 34 B(2) = (XA1*XA3*XA4)/D2 35 B(3) = (XA1*XA2*XA4)/D3 36 B(4) = (XA1*XA2*XA3)/D4 37 C(1) = (XA3*XA4 + XA2*XA4 + XA2*XA3)/D1 38 C(2) = (XA3*XA4 + XA1*XA4 + XA1*XA3)/D2 39 C(3) = (XA2*XA4 + XA1*XA4 + XA1*XA2)/D3 40 C(4) = (XA2*XA3 + XA1*XA3 + XA1*XA2)/D4 41 42 RETURN 43 END

1"""Python version of interp.f (YREC 7be4133): 4-point Lagrange weights 2and derivative weights.""" 3 4import numpy as np 5 6 7def interp(a: np.ndarray, b: np.ndarray, c: np.ndarray, x: float) -> None: 8 """Set b to the 4-point Lagrange interpolation weights at x for the 9 nodes a[1:5], and c to the weights of the first derivative. 10 11 Translates: interp.f#INTERP (7be4133) 12 13 Arguments: a (the nodes), x. 14 Arrays changed in place: b, c (index 0 unused). 15 16 Returns nothing. 17 18 """ 19 20 # The node differences and the denominators (Fortran lines 19–28). 21 a43 = a[4] - a[3] 22 a42 = a[4] - a[2] 23 a41 = a[4] - a[1] 24 a32 = a[3] - a[2] 25 a31 = a[3] - a[1] 26 a21 = a[2] - a[1] 27 d1 = -a21 * a31 * a41 28 d2 = a21 * a32 * a42 29 d3 = -a31 * a32 * a43 30 d4 = a41 * a42 * a43 31 32 # The interpolation weights (Fortran lines 29–36). 33 xa1 = x - a[1] 34 xa2 = x - a[2] 35 xa3 = x - a[3] 36 xa4 = x - a[4] 37 b[1] = (xa2 * xa3 * xa4) / d1 38 b[2] = (xa1 * xa3 * xa4) / d2 39 b[3] = (xa1 * xa2 * xa4) / d3 40 b[4] = (xa1 * xa2 * xa3) / d4 41 42 # The derivative weights (Fortran lines 37–40). 43 c[1] = (xa3 * xa4 + xa2 * xa4 + xa2 * xa3) / d1 44 c[2] = (xa3 * xa4 + xa1 * xa4 + xa1 * xa3) / d2 45 c[3] = (xa2 * xa4 + xa1 * xa4 + xa1 * xa2) / d3 46 c[4] = (xa2 * xa3 + xa1 * xa3 + xa1 * xa2) / d4