getalex06.f

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

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

Variables on this page (10)

capalex06, cln, id, it, numd, numt, r0gr06, tgr06, xe, ze

Purpose

GETALEX06 returns the low-temperature (molecular) opacity from the low-temperature opacities of Ferguson et al. (2005), called ALEX06 in the code (READALEX06; "Alexander low T opacities," per its header, MHP 3/09) for a given $\rho$, $T$, $X$, $Z$, with its logarithmic derivatives.

How it fits

Who calls it: only GETOPAC, when $\log T \le$ TMOLMAX and LALEX06 is set.
Inputs: DL, TL ($\log\rho$, $\log T$), X, Z.
What it calls: ALEX06TAB, which builds a table for the current $X$, $Z$ from the full set, and INTERP (four-point Lagrange weights and derivative weights).
Outputs: O ($\kappa$, cm$^2$ g$^{-1}$), OL ($\log\kappa$), QOD $=\partial\log\kappa/\partial\log\rho|_T$, QOT $=\partial\log\kappa/\partial\log T|_\rho$.
Modules: GALOT06_MOD (the grids TGR06, R0GR06, XXG06, ZZG06) and ALOT06_MOD (the table CAPALEX06 at the composition XE, ZE, and the last indices IT, ID), shared with ALEX06TAB and READALEX06; $\ln 10$ from CONSTANTS_MOD.

The method

If $X$ or $Z$ differs from the stored XE, ZE by more than $10^{-8}$, rebuild the table with ALEX06TAB (lines 26–33 ↗).

Question 1 (lines 26–33 ↗). Table rebuilt per composition. Every change of $X$ or $Z$ beyond $10^{-8}$ rebuilds the whole table (85 × 19 values, from 16 tables). Where the points below TMOLMAX differ in composition (with settling, or with a gradient in the cool envelope), this happens at almost every call. Is that cost acceptable, or should the routine interpolate in $X$ and $Z$ directly? A profile of a run with diffusion would show the share of time spent in ALEX06TAB.
The table is in $\log T$ and $\log R$, with $\log R = \log\rho - 3(\log T - 6)$ (that is, $R = \rho/T_6^3$).
Find the 4×4 stencil from the last indices, get Lagrange weights in $\log T$ and $\log R$, and interpolate $\log\kappa$ and $\partial\log\kappa/\partial\log R$.

Question 2 (lines 37–59 ↗; lines 66–88 ↗). No lower-edge handling. Below the smallest $\log R$ ($-8$), or outside the $\log T$ range (2.70 to 4.50), the stencil stays at the edge and the Lagrange weights extrapolate the cubic silently; only above the largest $\log R$ is the extrapolation made linear (step 4). Should there be a check, or a linear extrapolation at the other edges too? (Above $\log T$ = 4.50 is reached only if TMOLMAX is set above it; its default is 4.1.)
Above the largest $\log R$ of the table, extrapolate $\log\kappa$ linearly from the last two columns (lines 90–96 ↗; lines 113–120 ↗); the derivative is not extrapolated (bug 1).
Convert the temperature derivative from fixed $R$ to fixed $\rho$: $$\left.\frac{\partial\log\kappa}{\partial\log T}\right|_\rho = \left.\frac{\partial\log\kappa}{\partial\log T}\right|_R - 3\left.\frac{\partial\log\kappa}{\partial\log R}\right|_T.$$

Suspected bugs, checked by running the code

Each bug is marked by what running the code showed: (changes results, large) or (small), the fix moving a run clearly beyond a null change; (stops the run); (wrong output), the fix changing only what is written; (harmless in the runs), within the noise; (never occurred), the code ran but the fault did not; (not reached); (cannot be reached), a guard prevents it; (withdrawn), the suspicion was wrong.

1. (harmless in the runs) The derivative is not extrapolated with the value.

(lines 90–96 ↗; lines 108–111 ↗; lines 113–120 ↗; lines 126–129 ↗) Above the largest $\log R$ of the table ($\log R > 1$), $\log\kappa$ is extrapolated linearly from the last two columns, but QOPR, and so QOD $=\partial\log\kappa/\partial\log\rho|_T$, is the derivative of the cubic at the table edge, not the slope of the extrapolation. QOT is then converted to fixed $\rho$ with this QOD (line 129). What follows: There, the returned derivatives do not match the returned opacity; where they are used (the linearized equations) was not traced. Only with LALEX06 (off by default in the code; recommended by YREC's User Guide) and at $\log T \le$ TMOLMAX.

A possible fix:

getalex06.f

@@ lines 110–120 @@
      *           FXRD(3)*CAPALEX06(II,ID+2) +
      *           FXRD(4)*CAPALEX06(II,ID+3)
       END DO
+C 11/09/26 KC via Claude
+C code changed to give the extrapolated opacity its own derivative
+C in log R, the slope of the linear extrapolation.
       IF(LINEAR)THEN
          DO I = 1,4
             II = IT+I-1
             OP(I) = OP(I)+(RR-R)*(CAPALEX06(II,NUMD)-
      *      CAPALEX06(II,NUMD-1))/(R0GR06(NUMD)-R0GR06(NUMD-1))
+            QOPR(I) = (CAPALEX06(II,NUMD)-CAPALEX06(II,NUMD-1))/
+     *      (R0GR06(NUMD)-R0GR06(NUMD-1))
          END DO
          R = RR
       ENDIF

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): 34 runs, 26 of YREC's cases (from a 0.04 $M_\odot$ brown dwarf and the solar models to 9 $M_\odot$, rotation, and helium burning) and 8 made-up ones (each one of YREC's cases with a few settings changed), of which one stopped in its first second. The extrapolation ran in five of the 33 others, all low-mass stars: 12 million times in the 0.04 $M_\odot$ brown dwarf examples/run_from_dbl_to_zams/Test_m0040_feh0_noGS_ZAMS, fewer in the 0.1, 0.3, and 0.8 $M_\odot$ cases, at $\log T$ = 2.71 to 4.10 and $\log R$ = 1 to 15.6, that is, up to 14.6 beyond the table's last column. In the brown dwarf the derivative differed from the slope by more than 0.001 at 45 million of the 49 million temperature points (four a call). With the fix, the brown dwarf moves by up to $9.4\times10^{-5}$ dex in $\log L$ and $7.7\times10^{-6}$ dex in $\log T_c$, against $2.5\times10^{-5}$ and $1.5\times10^{-6}$ for the null change used as a noise floor (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$): four to five times the noise, so small, and not clearly real (748 models in both, end age 8,683 Gyr instead of 8,681). The track of testsuite/Test_m0p3feh+0p0GS98_scv_ZAMS is identical with the fix. The returned derivative disagrees with the returned opacity at 45 million of 49 million points, whatever the track does; its effect on the track is the same size as that of CONDOPACPINT's bug 4 on the same case, so likely the same amplified noise. Intent: none found.

Definitions

subroutineGETALEX06

lines 2–132 · GitHub

Ferguson et al. (2005) molecular opacity (ALEX06): four-point Lagrange interpolation in $\log T$ and $\log R$, with the derivative converted to fixed $\rho$ (method above).

Uses modules ALOT06_MOD, CONSTANTS_MOD, GALOT06_MOD

Called by GETOPAC

Flowchart of its calls, in order (3)

if (DIFZ.GT.1.0D-8 .OR. DIFX.GT.1.0D-8), lines 29–33

ALEX06TAB 32

INTERP 64→INTERP 100

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules, IMPLICIT NONE.

1C MHP 3/09 SUBROUTINE FOR CALCULATING ALEXANDER LOW T OPACITIES 2 SUBROUTINE GETALEX06(DL,TL,X,Z,O,OL,QOD,QOT) 3C 10/14/26 MHP via Claude 4C code changed to convert COMMON/GALOT06/ into GALOT06_MOD. 5 USE GALOT06_MOD 6C 10/14/26 MHP via Claude 7C code changed to convert COMMON/ALOT06/ into ALOT06_MOD. 8 USE ALOT06_MOD 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 11 USE CONSTANTS_MOD 12C 10/20/26 MHP via Claude 13C code changed to replace implicit typing with IMPLICIT NONE and 14C explicit declarations. 15 IMPLICIT NONE 16 REAL*8 DIFX,DIFZ,DL,O,OL,QOD,QOT,R,RR,TL,X,Z 17 INTEGER I,II,NUMD,NUMT 18 LOGICAL*4 LINEAR 19 PARAMETER (NUMT=85) 20 PARAMETER (NUMD=19) 21C 9/25/26 MHP via Claude 22C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 23C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 24 REAL*8 QR(4),FXT(4),FXTD(4),FXR(4),FXRD(4),OP(4),QOPR(4) 25 SAVE 26 DIFZ = ABS(Z-ZE) 27 DIFX = ABS(X-XE) 28C ENSURE THAT OPACITY TABLE HAS THE SAME X,Z VALUE AS THE POINT 29 IF(DIFZ.GT.1.0D-8 .OR. DIFX.GT.1.0D-8)THEN 30 ZE = Z 31 XE = X 32 CALL ALEX06TAB 33 ENDIF 34C COMPUTE R FOR GRID POINT 35 R = DL - 3.0D0*(TL-6.0D0) 36C Insure that index IT is within the required array bounds 37 IF(IT .LT. 1) IT=1 38 IF((IT+2) .GT. NUMT) IT=NUMT-2 39C FIND NEAREST GRID POINTS IN T. 40 IF(TL.LT.TGR06(IT+2))THEN 41 DO I = IT+1,2,-1 42 IF(TL.GT.TGR06(I))THEN 43 IT = I - 1 44 GOTO 10 45 ENDIF 46 END DO 47 IT = 1 48 10 10 CONTINUE 49 ELSE 50 DO I = IT+3,NUMT 51 IF(TL.LT.TGR06(I))THEN 52 IT = I - 2 53 IT = MIN(NUMT-3,IT) 54 GOTO 20 55 ENDIF 56 END DO 57 IT = NUMT - 3 58 20 20 CONTINUE 59 ENDIF 60C INTERPOLATION FACTORS IN LOG T 61 DO I = 1,4 62 QR(I) = TGR06(IT+I-1) 63 END DO 64 CALL INTERP(QR,FXT,FXTD,TL) 65C Insure that index ID is within the required array bounds 66 IF(ID .LT. 1) ID=1 67 IF((ID+2) .GT. NUMD) ID=NUMD-2 68C FIND NEAREST GRID POINTS IN R = RHO/T6**3 69 IF(R.LT.R0GR06(ID+2))THEN 70 DO I = ID+1,2,-1 71 IF(R.GT.R0GR06(I))THEN 72 ID = I - 1 73 GOTO 30 74 ENDIF 75 END DO 76 ID = 1 77 30 30 CONTINUE 78 ELSE 79 DO I = ID+3,NUMD 80 IF(R.LT.R0GR06(I))THEN 81 ID = I - 2 82 ID = MIN(NUMD-3,ID) 83 GOTO 40 84 ENDIF 85 END DO 86 ID = NUMD - 3 87 40 40 CONTINUE 88 ENDIF 89C INTERPOLATION FACTORS IN LOG R 90 IF(R.GT.R0GR06(NUMD))THEN 91 LINEAR = .TRUE. 92 RR = R 93 R = R0GR06(NUMD) 94 ELSE 95 LINEAR = .FALSE. 96 ENDIF 97 DO I = 1,4 98 QR(I) = R0GR06(ID+I-1) 99 END DO 100 CALL INTERP(QR,FXR,FXRD,R) 101C INTERPOLATE IN LOG R AT FIXED T 102 DO I = 1,4 103 II = IT + I - 1 104 OP(I) = FXR(1)*CAPALEX06(II,ID) + 105 * * FXR(2)*CAPALEX06(II,ID+1) + 106 * * FXR(3)*CAPALEX06(II,ID+2) + 107 * * FXR(4)*CAPALEX06(II,ID+3) 108 QOPR(I) = FXRD(1)*CAPALEX06(II,ID) + 109 * * FXRD(2)*CAPALEX06(II,ID+1) + 110 * * FXRD(3)*CAPALEX06(II,ID+2) + 111 * * FXRD(4)*CAPALEX06(II,ID+3) 112 END DO 113 IF(LINEAR)THEN 114 DO I = 1,4 115 II = IT+I-1 116 OP(I) = OP(I)+(RR-R)*(CAPALEX06(II,NUMD)- 117 * * CAPALEX06(II,NUMD-1))/(R0GR06(NUMD)-R0GR06(NUMD-1)) 118 END DO 119 R = RR 120 ENDIF 121C INTERPOLATE IN T 122 OL = FXT(1)*OP(1)+FXT(2)*OP(2)+FXT(3)*OP(3)+FXT(4)*OP(4) 123C D LN CAPPA/D LN T AT FIXED R 124 QOT = FXTD(1)*OP(1)+FXTD(2)*OP(2)+FXTD(3)*OP(3)+FXTD(4)*OP(4) 125C INTERPOLATE IN D LN CAPPA/ D LN R AT FIXED T 126 QOD = FXT(1)*QOPR(1) + FXT(2)*QOPR(2) + FXT(3)*QOPR(3) + 127 * * FXT(4)*QOPR(4) 128C CORRECT FROM DERIVATIVE AT FIXED R TO DERIVATIVE AT FIXED RHO 129 QOT = QOT - 3.0D0*QOD 130 O = EXP(CLN*OL) 131 RETURN 132 END

1"""Python version of getalex06.f (YREC 7be4133): the Ferguson-Alexander 22006 low-temperature opacity.""" 3 4import math 5 6import numpy as np 7 8from pyrec import alot06_mod 9from pyrec import constants_mod 10from pyrec import galot06_mod 11from pyrec.alex06tab import alex06tab 12from pyrec.interp import interp 13 14numt = 85 15numd = 19 16 17 18def getalex06(dl: float, tl: float, x: float, z: float, o: float, ol: float, qod: float, qot: float) -> tuple[float, 19 float, float, float]: 20 """Interpolate the Ferguson-Alexander 2006 opacity in log T and 21 log R = log rho - 3 (log T - 6), with four-point Lagrange weights, 22 from the table for the current X, Z (rebuilt by ALEX06TAB when they 23 change); extrapolate linearly above the largest log R, and convert 24 the temperature derivative to fixed rho. 25 26 Translates: getalex06.f#GETALEX06 (7be4133) 27 28 Arguments: dl, tl (log rho, log T), x, z; o, ol, qod, qot are only 29 set. 30 31 Returns (o, ol, qod, qot): kappa, log kappa, dlog kappa/dlog rho at 32 fixed T, dlog kappa/dlog T at fixed rho. 33 34 """ 35 36 a = alot06_mod 37 g = galot06_mod 38 qr = np.zeros(4 + 1) 39 fxt = np.zeros(4 + 1) 40 fxtd = np.zeros(4 + 1) 41 fxr = np.zeros(4 + 1) 42 fxrd = np.zeros(4 + 1) 43 op = np.zeros(4 + 1) 44 qopr = np.zeros(4 + 1) 45 46 # Ensure that the table has the point's X, Z (Fortran lines 26–33). 47 # Question (see the page): the table is rebuilt at every change of 48 # composition. 49 difz = abs(z - a.ze) 50 difx = abs(x - a.xe) 51 if difz > 1.0e-8 or difx > 1.0e-8: 52 a.ze = z 53 a.xe = x 54 alex06tab() 55 56 # R for the grid point (Fortran line 35). 57 r = dl - 3.0 * (tl - 6.0) 58 59 # The 4-point stencil in log T, searched from the last one (Fortran 60 # lines 36–59). 61 # Question (see the page): no check below or above the table. 62 if a.it < 1: 63 a.it = 1 64 if a.it + 2 > numt: 65 a.it = numt - 2 66 if tl < g.tgr06[a.it + 2]: 67 for i in range(a.it + 1, 1, -1): 68 if tl > g.tgr06[i]: 69 a.it = i - 1 70 break # GOTO 10 71 else: 72 a.it = 1 73 else: 74 for i in range(a.it + 3, numt + 1): 75 if tl < g.tgr06[i]: 76 a.it = i - 2 77 a.it = min(numt - 3, a.it) 78 break # GOTO 20 79 else: 80 a.it = numt - 3 81 82 # Interpolation factors in log T (Fortran lines 60–64). 83 for i in range(1, 4 + 1): 84 qr[i] = g.tgr06[a.it + i - 1] 85 interp(qr, fxt, fxtd, tl) 86 87 # The same in log R (Fortran lines 65–88). 88 if a.id < 1: 89 a.id = 1 90 if a.id + 2 > numd: 91 a.id = numd - 2 92 if r < g.r0gr06[a.id + 2]: 93 for i in range(a.id + 1, 1, -1): 94 if r > g.r0gr06[i]: 95 a.id = i - 1 96 break # GOTO 30 97 else: 98 a.id = 1 99 else: 100 for i in range(a.id + 3, numd + 1): 101 if r < g.r0gr06[i]: 102 a.id = i - 2 103 a.id = min(numd - 3, a.id) 104 break # GOTO 40 105 else: 106 a.id = numd - 3 107 108 # Above the table's largest log R: interpolate at the edge, then 109 # extrapolate linearly (Fortran lines 89–100). 110 if r > g.r0gr06[numd]: 111 linear = True 112 rr = r 113 r = g.r0gr06[numd] 114 else: 115 linear = False 116 for i in range(1, 4 + 1): 117 qr[i] = g.r0gr06[a.id + i - 1] 118 interp(qr, fxr, fxrd, r) 119 120 # Interpolate in log R at fixed T (Fortran lines 101–112). 121 cap = a.capalex06 122 idd = a.id 123 for i in range(1, 4 + 1): 124 ii = a.it + i - 1 125 op[i] = (fxr[1] * cap[ii, idd] + fxr[2] * cap[ii, idd + 1] 126 + fxr[3] * cap[ii, idd + 2] + fxr[4] * cap[ii, idd + 3]) 127 qopr[i] = (fxrd[1] * cap[ii, idd] + fxrd[2] * cap[ii, idd + 1] 128 + fxrd[3] * cap[ii, idd + 2] + fxrd[4] * cap[ii, idd + 3]) 129 130 # The linear extrapolation in log R (Fortran lines 113–120). 131 if linear: 132 # Suspected bug (see the page): qopr, and so qod, is not 133 # extrapolated. 134 for i in range(1, 4 + 1): 135 ii = a.it + i - 1 136 op[i] = op[i] + (rr - r) * (cap[ii, numd] - cap[ii, numd - 1]) / ( 137 g.r0gr06[numd] - g.r0gr06[numd - 1]) 138 r = rr 139 140 # Interpolate in T (Fortran lines 121–122). 141 ol = fxt[1] * op[1] + fxt[2] * op[2] + fxt[3] * op[3] + fxt[4] * op[4] 142 143 # dlog kappa/dlog T at fixed R, then at fixed rho (Fortran lines 144 # 123–130). 145 qot = (fxtd[1] * op[1] + fxtd[2] * op[2] + fxtd[3] * op[3] 146 + fxtd[4] * op[4]) 147 qod = (fxt[1] * qopr[1] + fxt[2] * qopr[2] + fxt[3] * qopr[3] 148 + fxt[4] * qopr[4]) 149 qot = qot - 3.0 * qod 150 o = math.exp(constants_mod.cln * ol) 151 return o, ol, qod, qot