intrp2.f

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

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

Purpose

INTRP2 returns the weights of 4-point Lagrange interpolation at a point $x$, without the derivative weights: the same as INTERP minus its C output.

How it fits

Called by the opacity-table routines ALEX06TAB, GETOPAL95, OP95XTAB and OP95ZTAB (ALX8TH, ALXZTAB and YALO3D also called it; removed since rp_1), by the diffusion regridding EQUAL_TO_MODEL and MICRODIFF_ETM, and by PHYSIC. No COMMON blocks, no calls.

The header has a typo and a copied comment. The banner says "INTEP2" (line 4 ↗), and the comment, copied from INTERP, still speaks of the VandenBerg opacity tables only (lines 15–17 ↗).

The method

$B_i(x) = \prod_{j\ne i}(x-a_j)/(a_i-a_j)$ for $i = 1..4$.

Definitions

subroutineINTRP2

lines 6–37 · GitHub

4-point Lagrange weights only.

Called by ALEX06TAB ×2, EQUAL_TO_MODEL, GETOPAL95 ×3, MICRODIFF_ETM, OP95XTAB, OP95ZTAB ×3, PHYSIC

Change log

The changes since rp_1

Changed since rp_1. Declarations only: IMPLICIT NONE with explicit declarations. Three of its callers (ALX8TH, ALXZTAB, YALO3D) are removed since rp_1.

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

1"""Python version of intrp2.f (YREC development, 7be4133).""" 2 3import numpy as np 4 5 6def intrp2(a: np.ndarray, b: np.ndarray, x: float) -> None: 7 """Weights of 4-point Lagrange interpolation at ``x``, no derivatives. 8 9 The same as INTERP without its derivative weights. Used by the opacity 10 tables, the diffusion regridding and PHYSIC. 11 12 Translates: intrp2.f#INTRP2 (7be4133) 13 14 Parameters 15 ---------- 16 a : ndarray, shape (5,) 17 The four nodes, ``a[1..4]``. 18 b : ndarray, shape (5,) 19 Set in place: the weights ``b[1..4]``, with 20 ``f(x) ~ sum_i b[i] f(a[i])``. 21 x : float 22 The point of interpolation. 23 24 """ 25 26 # SAVE in the Fortran is harmless: every local is set before use. 27 # The node differences and the denominators (Fortran lines 18–27). 28 a43 = a[4] - a[3] 29 a42 = a[4] - a[2] 30 a41 = a[4] - a[1] 31 a32 = a[3] - a[2] 32 a31 = a[3] - a[1] 33 a21 = a[2] - a[1] 34 d1 = -a21 * a31 * a41 35 d2 = a21 * a32 * a42 36 d3 = -a31 * a32 * a43 37 d4 = a41 * a42 * a43 38 39 # The interpolation weights (Fortran lines 28–35). 40 xa1 = x - a[1] 41 xa2 = x - a[2] 42 xa3 = x - a[3] 43 xa4 = x - a[4] 44 b[1] = (xa2 * xa3 * xa4) / d1 45 b[2] = (xa1 * xa3 * xa4) / d2 46 b[3] = (xa1 * xa2 * xa4) / d3 47 b[4] = (xa1 * xa2 * xa3) / d4