getopac.f

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

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

Variables on this page (10)

ishort, l2z, lalex06, lcondopacp, lopal95, lpurez, optol, tmolmax, tmolmin, zenv

Purpose

GETOPAC is YREC's opacity routine: for a given $\rho$, $T$, $X$, $Z$ it combines the low-temperature (molecular) table and the interior table, ramps between them, handles metal-rich helium-burning cores, and adds conductive opacity. It returns $\kappa$ and its logarithmic derivatives.

How it fits

Who calls it: seven routines: COEFFT, ENVINT, PHYSIC, QATM, QENV, SCONVEC and STARIN.
Inputs: DL, TL ($\log\rho$, $\log T$), X, Z, and FXION (ionization fractions, passed on to the conductive opacity).

Question 1 (line 192 ↗; condopacpint.f, lines 32–40 ↗). FXION used only by conduction. Fine, but its meaning differs by equation of state (EQSAHA versus SCV); does CONDOPACPINT depend on which one supplied it?
What it calls: the low-temperature table GETALEX06 (the low-temperature tables of Ferguson et al. 2005, the code's ALEX06); the interior table GETOPAL95 (OPAL 1995); GTPURZ (pure-metal table); and CONDOPACPINT (Potekhin conductive opacity).
Outputs: O ($\kappa$), OL ($\log\kappa$), QOD $=\partial\log\kappa/\partial\log\rho$, QOT $=\partial\log\kappa/\partial\log T$.
Modules: NMLPAR_MOD (the namelist switches LALEX06, LOPAL95, LPUREZ, LcondOpacP, the ramp limits TMOLMIN, TMOLMAX, and the tolerance OPTOL, read by PARMIN); NEWOPAC_MOD (L2Z, set by SETUPOPAC from LZRAMP or LDIFZ); COMP_MOD (ZENV); IOUNITS_MOD (ISHORT, for the error messages).

The method

Molecular opacity for $\log T \le$ TMOLMAX, from the Ferguson et al. (2005) tables (LALEX06); if they are not switched on, the run stops ("NO MOLECULAR OPACITY TABLE CHOSEN").
Interior opacity, unless $\log T <$ TMOLMIN (where the molecular value is used alone):
for $Z > 0.1$ and $\log T > 7.0$ (helium-burning cores), linear interpolation in $Z$ between OPAL 1995 at $Z_\mathrm{tab} = 0.1$ and the pure-metal table at $Z = 1$: $\log\kappa = \log\kappa_{Z=1} + (Z-1)\,\dfrac{\log\kappa_\mathrm{tab} - \log\kappa_{Z=1}}{Z_\mathrm{tab} - 1}$, and the derivatives likewise; the run stops if no pure-$Z$ table is loaded (LPUREZ) or if OPAL 1995 is not switched on;
the run stops if $Z > 0.12$ below that temperature, or if $Z$ differs from the envelope's ZENV beyond OPTOL with a single-$Z$ table (neither L2Z nor OPAL 1995, which interpolates in $Z$ itself; bug 2 below);
otherwise OPAL 1995; if it is not switched on, the run stops ("NO OPACITY TABLE CHOSEN").
Ramp for TMOLMIN ≤ $\log T \le$ TMOLMAX: with $w = (\log T - \mathrm{TMOLMIN})/(\mathrm{TMOLMAX} - \mathrm{TMOLMIN})$, $\kappa = w\,\kappa_\mathrm{int} + (1-w)\,\kappa_\mathrm{mol}$, and the derivatives likewise. The ramp is linear in $\kappa$, while the derivatives (of $\log\kappa$) are averaged with the same weights, so they are not exactly those of the ramped $\kappa$: The weights should be $w\,\kappa_\mathrm{int}/\kappa$ and $(1-w)\,\kappa_\mathrm{mol}/\kappa$, and the term $(\partial w/\partial\log T)(\kappa_\mathrm{int} - \kappa_\mathrm{mol})/\kappa$ is dropped. This is known, and affects the Jacobian only (settled from YREC's notes, Nov 9, 2026: KNOWN_ISSUES.md, issue 40, "Opacity ramp derivatives"; lines 177–181 ↗).
Conduction. Potekhin's opacity if LcondOpacP is set and it succeeds; otherwise, for $\log T \ge 4.2$ and $\log\rho \ge 2\log T - 13$, the fit to Hubbard & Lampe (1969) as coded (the code names no source for the fit; see the question below): $$\log\kappa_c = \log(1 - 0.6X) - 14.6196 - (3.5853 + 0.1386\log\rho)\log\rho + (5.1324 - 0.3219\log T)\log T + 0.3901\log\rho\log T.$$ The two cuts leave out cool gas and gas of low density for its temperature ($\rho < 10^{-13}T^2$). The derivatives QODC, QOTC it receives from Potekhin's fit carry CONDOPACPINT's suspected bugs 1 and 2, and its extrapolation below $\log\rho = -6$ bugs 3 and 4.

Question 2 (lines 200–211 ↗). The conductive fit's source. The fallback fit to Hubbard & Lampe (1969) has no reference in the code. Whose fit is it, and what range of $\rho$, $T$ was it made for? Are the cuts ($\log T \ge 4.2$, $\rho \ge 10^{-13}T^2$) meant to leave out gas where the electrons are far from degenerate and conduction is negligible next to radiation?

Question 3 (lines 190–199 ↗). Conduction where no table applies. Below $\log T = 4.2$ or at low density the routine returns without conduction, which is fine; but when Potekhin is enabled and fails, Hubbard–Lampe is used without a message. Should the fallback be reported?
Combination: $1/\kappa = 1/\kappa_\mathrm{rad} + 1/\kappa_c$, and $$\frac{\partial\log\kappa}{\partial\log\rho} = \frac{\kappa_c\,\partial\log\kappa_\mathrm{rad} + \kappa_\mathrm{rad}\,\partial\log\kappa_c}{\kappa_\mathrm{rad} + \kappa_c},$$ and the same in $\log T$.

The tables' mixtures. Each table is computed for one heavy-element mixture, named in its file (FLIV95, FALEX06): in the test suite GS98.OP17 (OP, Grevesse & Sauval 1998, the YREC paper's base mixture), GN93hz, AAG2021zx0.0187.tab or Magg22_met_zx0.0226.tab, and alexmol06gs98.tab for the molecular table. The run's XENV, ZENV (and the equation of state's mixture) should match the tables', since the code only interpolates in $X$ and $Z$; the documentation (Control Options) recommends tying the CNO abundances "to preset mixtures used in opacity tables." The pure-metal table (FPUREZ) is input/opacity/lanl/PURECO.DBGLAOL, the table the User Guide names for helium burning ("Helium burning uses a different type of atomic opacity table"); so in a helium-burning core the step-2 interpolation runs from OPAL at $Z = 0.1$ to this table at $Z = 1$ (what mixture it holds is the question below).

Question 4 (lines 109–138 ↗). OPAL's C/O-rich tables. The YREC paper (§2.3.2) says OPAL's tables for carbon- and oxygen-rich mixtures (Iglesias & Rogers 1993) are used in helium-burning regions; the code instead interpolates in $Z$ between OPAL at $Z = 0.1$ and the table PURECO.DBGLAOL (step 2), which the tables' README describes as interpolating "between helium and a c/o mixture" (input/opacity/atomic/README, line 1): With what C/O ratio, and how does it interpolate between helium and C/O? Is the paper describing the older scheme, or is the C/O table meant to stand for it? A reader would ask how large the error in the core opacity of horizontal-branch and AGB models is when the C/O ratio of the core differs from the table's, and when $0.1 < Z < 1$ is reached by a straight line in $\log\kappa$.

From YREC's documentation

Control Options (Input files) and Physics Options (Opacity options): LALEX06 uses the molecular opacities of Ferguson et al. (2005) and LOPAL95 the atomic opacities in OP/OPAL format (both recommended); LPUREZ adds the pure Y/Z table (FPUREZ); LCONDOPACP uses the updated conductive opacities. Below $\log T$ = TMOLMIN (default 4.0) only the molecular opacities are used, above TMOLMAX (4.1) only the atomic ones, with a ramp between. The code builds a 2D table for the surface abundances, useful in outer convection zones and the envelope, and switches to the full 3D/4D interpolation when $X$ or $Z$ differ from it by more than OPTOL (1.0D-8). LOPAL95 has no default in NMLPAR_MOD (nor has ZOPAL951); the User Guide asks the user to set LOPAL95 = .TRUE. with the table's path; a namelist that omits it leaves it to the compiler (usually false), and the run then stops at a trap with status 3 (settled from the code and the documentation, Nov 9, 2026: docs/source/userguide.md, line 34; nmlpar_mod.f, line 243 ↗; lines 126–130 ↗, lines 162–166 ↗).

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. (not reached) A stale temperature in the pure-$Z$ trap's message.

(lines 104–114 ↗) Since the helium-burning branch starts at $\log T > 7.0$ ("JCZ 211125," line 107), the message for a missing pure-$Z$ table still says "Z>0.10 T > 5 X 10^7 K," and the comment above the test still speaks of "50 million K." A run that stops here with $10^7 < T < 5\times10^7$ K reports a temperature it has not reached. The trap fires only when LPUREZ is off (its default, nmlpar_mod.f, line 244 ↗) and a core with $Z > 0.1$ passes $\log T = 7.0$.

A possible fix:

getopac.f

@@ lines 108–114 @@
 C     semiconvection+overshoot HB models, which can reach lower core temperatures
       IF((Z .GT. 0.1D0) .AND. (TL.GT.7.0D0)) THEN
          IF(.NOT.LPUREZ) THEN
-            WRITE(ISHORT, *)' ERROR: Z>0.10 T > 5 X 10^7 K',
+C 11/09/26 KC via Claude
+C code changed to give the branch's threshold of 10^7 K in the message.
+            WRITE(ISHORT, *)' ERROR: Z>0.10 T > 10^7 K',
      *        ' NEED PURE Z TABLE TO CONTINUE. Z,LOG T=',Z, TL
 C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38).
               STOP 1

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): not reached in any of 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 trap needs LPUREZ off, and every one of YREC's cases sets it on; what would reach it: LPUREZ off, with a helium-burning core whose $Z$ passes 0.1 above $\log T = 7.0$. The fix only changes a message, so it was not run. Intent: The comment at line 142 shows the other message was updated with the threshold ("JCZ 211125 changed to 10^7 K in message"), but not this one.

2. (not reached) The single-table test now only changes the stop.

(lines 140–141 ↗, lines 162–166 ↗) The test ABS(Z-ZENV) .GT. OPTOL .AND. .NOT.L2Z .AND. .NOT.LOPAL95 can be true only when OPAL 1995 is off, the one interior table left. Then a run whose $Z$ differs from ZENV stops with "OUTSIDE SINGLE TABLE USED" and status 1 (a run-time error), while the same run with $Z$ = ZENV stops at the final trap with "NO OPACITY TABLE CHOSEN" and status 3 (an input error), which is the true cause in both. L2Z is read nowhere else (set in setupopac.f, line 64 ↗). Taken only when LOPAL95 is off.

A possible fix:

getopac.f

@@ lines 137–148 @@
        SLOPE = (QOT-QOTZ)/(ZIT-1.0D0)
        QOT = QOTZ + (Z-1.0D0)*SLOPE
 C      ELSE
-      ELSE IF((Z.GT.0.12D0) .OR. ((ABS(Z-ZENV) .GT. OPTOL)
-     *      .AND..NOT.L2Z .AND. .NOT.LOPAL95))THEN
+C 11/09/26 KC via Claude
+C code changed to drop the single-table test, which can hold only when
+C no interior table is chosen; the final trap then stops the run.
+      ELSE IF(Z.GT.0.12D0)THEN
 C     JCZ 211125 changed to 10^7 K in message to reflect above change in logic.
          WRITE(ISHORT,*)' Z>0.12 T < 10^7 K',
-     *   ' OUTSIDE OPAL OPACITY TABLE RANGE OR Z',
-     *   ' OUTSIDE SINGLE TABLE USED.Z,ZENV,LOG T=',Z,ZENV,TL
+     *   ' OUTSIDE OPAL OPACITY TABLE RANGE.Z,ZENV,LOG T=',Z,ZENV,TL
 C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38).
          STOP 1
 C

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged): not reached in any of the 34 runs; every one of YREC's cases sets LOPAL95 on. What would reach it: LOPAL95 off, which stops the run at the first interior point either way. The fix only changes which stop fires, so it was not run. Intent: none found.

Definitions

subroutineGETOPAC

lines 6–225 · GitHub

Ramps and combines the molecular, interior, pure-metal and conductive opacities (method above).

Uses modules COMP_MOD, IOUNITS_MOD, NEWOPAC_MOD, NMLPAR_MOD

Called by COEFFT, ENVINT, PHYSIC, QATM, QENV, SCONVEC ×3, STARIN

Flowchart of its calls, in order (5)

if (TL.LE.TMOLMAX), lines 79–95

if (LALEX06), lines 80–88

GETALEX06 81

if ((Z .GT. 0.1D0) .AND. (TL.GT.7.0D0)), lines 109–139

GTPURZ 116

if (LOPAL95), lines 117–125

GETOPAL95 120

else if (LOPAL95), lines 154–161

GETOPAL95 155

if (LcondOpacP), lines 190–192

CONDOPACPINT 192

Change log

The changes since rp_1

Changed since rp_1. The Alexander 1995, Kurucz 1990, OPAL 1992 and LAOL 1989 tables were removed, so OPAL 1995 is the only interior table and the Ferguson et al. (2005) tables (ALEX06) the only molecular ones; their two-$Z$ interpolations went with them, and so did the Kurucz failure path (the alternate return to label 100). New traps stop the run when no molecular table is chosen below TMOLMAX, and when no interior table is chosen in the helium-burning branch; this settles question 1 of rp_1 (ZIT undefined in that branch), which is removed. The helium-burning branch now starts at $\log T > 7.0$ instead of 7.7 ("JCZ 211125," for semiconvection and overshoot horizontal-branch models), and the message for $Z > 0.12$ says $10^7$ K. COMMON blocks became modules (NEWOPAC, NWLAOL and MISCOPAC to NMLPAR_MOD and NEWOPAC_MOD; OPTAB removed), IMPLICIT NONE.

Changed since a8a5950. The stops now exit with a nonzero status instead of a bare STOP (KNOWN_ISSUES.md, issue 38): status 3, the code for input errors, when no molecular or interior table is chosen (three traps), and status 1, the code for run-time errors, when a helium-burning core needs the missing pure-$Z$ table and when $Z > 0.12$ or $Z$ leaves a single-$Z$ table.

1C 2C 3CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 4C GETOPAC 5CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 6 SUBROUTINE GETOPAC(DL,TL,X,Z,O,OL,QOD,QOT,FXION) 7C 10/14/26 MHP via Claude 8C code changed to convert COMMON/NEWOPAC/ into NEWOPAC_MOD. 9 USE NEWOPAC_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/COMP/ into COMP_MOD. 12 USE COMP_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 15 USE NMLPAR_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 COND,DL,O,OC,OCL,OL,OLZ,OX,OZ,QOD,QODC,QODX,QODZ,QOT,QOTC, 24 * * QOTX,QOTZ,RMPWT,SLOPE,SO,SOL,SQOD,SQOT,TL,X,Z,ZIT 25 LOGICAL*4 LCONDO,LGOTATM 26 REAL*8 FXION(3) 27C MHP 8/25 Removed unused variables 28C CHARACTER*256 FLAOL, FPUREZ, FKUR2, FcondOpacP 29C OPACITY COMMON BLOCKS - modified 3/09 30C 9/23/26 MHP via Claude 31C code changed to remove the OPAL92 opacity tables 32C 9/23/26 MHP via Claude 33C code changed to remove the Kurucz 1990 opacity tables 34C 9/23/26 MHP via Claude 35C code changed to remove the LAOL89 opacity tables 36C 9/24/26 MHP via Claude 37C code changed to remove the Alexander 1995 opacity tables 38C 9/24/26 MHP via Claude 39C code changed to remove dead opacity-table Z-value inputs 40C 9/24/26 MHP via Claude 41C code changed to merge COMMON /MISCOPAC/ and COMMON/NWLAOL/ into 42C COMMON /NEWOPAC/ 43C 9/25/26 MHP via Claude 44C code changed to consolidate all I/O logical unit numbers, previously 45C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 46C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 47C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 48C COMMON/IOUNITS/ 49C 10/12/26 MHP via Claude 50C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 51C namelist variables that the physics routines read from 27 COMMON 52C blocks which also held derived or working members. Those blocks keep 53C only their non-namelist members; five of them -- NEWMX, VMULT, 54C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 55C Derivation-only namelist variables from the same blocks went to 56C COMMON/NMLBYPASS/ instead, not here. 57C 58C Declared identically in all 121 files that hold any of these members. 59C Generated from one member list, so every copy is byte-identical; the 60C linker's 'size of symbol changed' warning is the check that they are. 61C 10/12/26 MHP via Claude 62C code changed to remove COMMON/OPTAB/ entirely. Once OPTOL moved to the 63C consolidated namelist block, every remaining member was write-only: 64C ZSI appeared only in a DATA statement in parmin.f, and IDT and IDD were 65C assigned the constants 15 and 5 in coefft.f, envint.f, physic.f and 66C starin.f and read nowhere in the code. They are leftovers from the 67C removed opacity-table interpolation schemes. 68 SAVE 69C 70C THIS SUBROUTINE CALCULATES THE OPACITY FOR A GIVEN X AND Z. 71C IF LDIFZ=T OR LZRAMP=T THEN INTERPOLATE BETWEEN TWO Z TABLES. 72C IN A SMALL T RANGE THE ATMOSPHERE AND INTERIOR OPACITY ARE 73C RAMPED FROM ONE TO THE OTHER. 74C 75C 76C GET ATMOSPHERE OPACITY 77C 78 LGOTATM = .FALSE. 79 IF(TL.LE.TMOLMAX)THEN 80 IF (LALEX06) THEN 81 CALL GETALEX06(DL,TL,X,Z, SO, SOL,SQOD,SQOT) 82 LGOTATM = .TRUE. 83C 9/23/26 MHP via Claude 84C code changed to remove the Kurucz 1990 opacity tables 85C 9/23/26 MHP via Claude 86C added trap for no molecular (low-T) opacity table chosen 87C 9/24/26 MHP via Claude 88C code changed to remove the Alexander 1995 opacity tables 89 ELSE 90 WRITE(ISHORT,*)'NO MOLECULAR OPACITY TABLE CHOSEN', 91 * * ' RUN STOPPED. X Z TL=',X,Z,TL 92C 10/30/26 MHP via Claude: bare STOP -> STOP 3, the input-error code (issue 38). 93 STOP 3 94 END IF 95 ENDIF 96C 10/20/26 MHP via Claude 97C code changed to drop an unused "100 CONTINUE"; nothing branched to it. 98 99C GET INTERIOR OPACITY IF NEEDED 100 101 IF (TL .LT. TMOLMIN.AND.LGOTATM) GOTO 1000 102 103C HELIUM BURNING REGION (HB EVOLUTION) USE PURE Z TABLE 104C mhp 7/12 Altered logic of the opacities in the He burnng 105C regime. Switched to exclusive usage of OPAL below 50 million K 106C and switched the ramp to above Z = 0.1. 107c$$$ JCZ 211125 changing temperature limit to 7.0 to accommodate 108C semiconvection+overshoot HB models, which can reach lower core temperatures 109 IF((Z .GT. 0.1D0) .AND. (TL.GT.7.0D0)) THEN 110 IF(.NOT.LPUREZ) THEN 111 WRITE(ISHORT, *)' ERROR: Z>0.10 T > 5 X 10^7 K', 112 * * ' NEED PURE Z TABLE TO CONTINUE. Z,LOG T=',Z, TL 113C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 114 STOP 1 115 END IF 116 CALL GTPURZ(DL,TL,OZ,OLZ,QODZ,QOTZ) 117 IF (LOPAL95) THEN 118C MHP 7/12 INTERPOLATE TO MAXIMUM Z IN TABLE 119 ZIT = 0.1D0 120 CALL GETOPAL95(DL,TL,X,ZIT,O,OL,QOD,QOT) 121C 9/23/26 MHP via Claude 122C code changed to remove the OPAL92 opacity tables 123C 9/23/26 MHP via Claude 124C code changed to remove the LAOL89 opacity tables; added trap 125C for no opacity table chosen 126 ELSE 127 WRITE(ISHORT,*)'NO OPACITY TABLE CHOSEN', 128 * * ' RUN STOPPED. X Z TL=',X,Z,TL 129C 10/30/26 MHP via Claude: bare STOP -> STOP 3, the input-error code (issue 38). 130 STOP 3 131 END IF 132 SLOPE = (OL-OLZ)/(ZIT-1.0D0) 133 OL = OLZ + (Z-1.0D0)*SLOPE 134 O = 10.0D0**OL 135 SLOPE = (QOD-QODZ)/(ZIT-1.0D0) 136 QOD = QODZ + (Z-1.0D0)*SLOPE 137 SLOPE = (QOT-QOTZ)/(ZIT-1.0D0) 138 QOT = QOTZ + (Z-1.0D0)*SLOPE 139C ELSE 140 ELSE IF((Z.GT.0.12D0) .OR. ((ABS(Z-ZENV) .GT. OPTOL) 141 * * .AND..NOT.L2Z .AND. .NOT.LOPAL95))THEN 142C JCZ 211125 changed to 10^7 K in message to reflect above change in logic. 143 WRITE(ISHORT,*)' Z>0.12 T < 10^7 K', 144 * * ' OUTSIDE OPAL OPACITY TABLE RANGE OR Z', 145 * * ' OUTSIDE SINGLE TABLE USED.Z,ZENV,LOG T=',Z,ZENV,TL 146C 10/30/26 MHP via Claude: bare STOP -> STOP 1, the run-time error code (issue 38). 147 STOP 1 148C 149C NOT HELIUM BURNING REGION (HB EVOLUTION) OR L2Z=T AND 150C Z STILL NOT TOO LARGE IN CORE (<.15) SO CAN USE 151C SECOND Z TABLE RATHER THAN PURE Z TABLE 152C 153C IF (LOPAL95) THEN 154 ELSE IF (LOPAL95) THEN 155 CALL GETOPAL95(DL,TL,X,Z,O,OL,QOD,QOT) 156C 9/23/26 MHP via Claude 157C code changed to remove the OPAL92 opacity tables 158C 9/23/26 MHP via Claude 159C code changed to remove the LAOL89 opacity tables 160C MHP 7/12 INSERT FINAL TRAP - NO OPACITY COMPUTED 161C SHOULD NOT BE ABLE TO GET HERE. 162 ELSE 163 WRITE(ISHORT,*)'NO OPACITY TABLE CHOSEN', 164 * * ' RUN STOPPED. X Z TL=',X,Z,TL 165C 10/30/26 MHP via Claude: bare STOP -> STOP 3, the input-error code (issue 38). 166 STOP 3 167 END IF 168C END IF 169C 170 1000 1000 CONTINUE 171C DO A RAMP BETWEEN SURFACE AND INTERIOR OPACITY 172 173 IF (LGOTATM .AND. TL .LE. TMOLMAX) THEN 174 IF( TL.GE.TMOLMIN) THEN 175C RR = DL - 3.0D0*(TL-6.0D0) 176C WRITE(*,*)RR,TL,OL,SOL 177 RMPWT = (TL-TMOLMIN)/(TMOLMAX-TMOLMIN) 178 O = RMPWT*O + (1.0D0-RMPWT)*SO 179 OL = DLOG10(O) 180 QOD = RMPWT*QOD + (1.0D0-RMPWT)*SQOD 181 QOT = RMPWT*QOT + (1.0D0-RMPWT)*SQOT 182 ELSE 183 O = SO 184 OL = SOL 185 QOD = SQOD 186 QOT = SQOT 187 END IF 188 END IF 189C DO CONDUCTIVE OPACITY CORRECTION 190 IF (LcondOpacP) THEN 191C Get Potekhin conductive opacity 192 CALL condOpacPInt(DL,TL,X,Z,OC,OCL,QODC,QOTC,FXION,LCONDO) 193 ELSE 194 LCONDO = .FALSE. 195 END IF 196 197 IF (.NOT. LCONDO) THEN 198C If we get here we have no Potekhin opacity, so we 199C try for Hubbard Lampe 200 IF(TL.LT.4.2D0) THEN 201 RETURN 202 ELSE IF(DL.LT.(2.0D0*TL-13.0D0)) THEN 203 RETURN 204 ELSE 205C Do Hubbard Lampe conductive opacity calculation 206 COND = DLOG10(1.0D0-0.6D0*X) - 14.6196D0 - 207 * * (3.5853D0+0.1386D0*DL)*DL + 208 * * (5.1324D0-0.3219D0*TL)*TL + 0.3901D0*DL*TL 209 OC = 10.0D0**COND 210 QODC = 0.3901D0*TL - 0.2772D0*DL - 3.5853D0 211 QOTC = 0.3901D0*DL - 0.6438D0*TL + 5.1324D0 212 ENDIF 213 ENDIF 214 215 OX = O ! Save the radiative opacity stuff 216 QODX = QOD 217 QOTX = QOT 218C Add the opacities apppropriately 219 O = OX*OC/(OX + OC) ! e.g. 1/O = 1/OX + 1/OC 220 OL = DLOG10(O) 221 QOD = (QODX + QODC - (OX*QODX + OC*QODC)/(OX + OC)) 222 QOT = (QOTX + QOTC - (OX*QOTX + OC*QOTC)/(OX + OC)) 223 224 RETURN 225 END

1"""Python version of getopac.f (YREC development, 7be4133).""" 2 3import math 4import sys 5import numpy as np 6 7from pyrec import comp_mod 8from pyrec import newopac_mod 9from pyrec import nmlpar_mod 10from pyrec.condopacpint import condopacpint 11from pyrec.getalex06 import getalex06 12from pyrec.getopal95 import getopal95 13from pyrec.gtpurz import gtpurz 14 15 16def _stop(code, *items): 17 """WRITE(ISHORT,*) the message, then STOP with the exit code: 3 for an 18 input error (no table chosen), 1 for a run-time error. 19 20 The unit ISHORT (from iounits_mod) is the run's short output file; 21 printing to standard error stands in for it here. 22 23 """ 24 25 print(*items, file=sys.stderr) 26 sys.exit(code) 27 28 29def getopac(dl: float, tl: float, x: float, z: float, o: float, ol: float, qod: float, qot: float, 30 fxion: np.ndarray) -> tuple[float, float, float, float]: 31 """Opacity and its logarithmic derivatives for a given rho, T, X, Z. 32 33 Combines the molecular (Alexander 2006) and interior (OPAL 1995) 34 tables with a ramp between TMOLMIN and TMOLMAX, interpolates towards 35 the pure-metal table in metal-rich helium-burning cores, and adds 36 conductive opacity (Potekhin, or the Hubbard-Lampe fit). 37 38 Translates: getopac.f#GETOPAC (7be4133) 39 40 Parameters 41 ---------- 42 dl, tl : float 43 log10 density and log10 temperature. 44 x, z : float 45 Hydrogen and metal mass fractions. 46 o, ol, qod, qot : float 47 Output only. 48 fxion : ndarray, shape (4,) 49 Ionization fractions ``fxion[1..3]``, passed to the conductive 50 opacity. 51 52 Returns 53 ------- 54 o, ol, qod, qot : float 55 kappa, log10 kappa, dlog kappa/dlog rho, dlog kappa/dlog T. 56 57 """ 58 59 nml = nmlpar_mod 60 61 # Get atmosphere (molecular) opacity (Fortran lines 78–95). 62 lgotatm = False 63 if tl <= nml.tmolmax: 64 if nml.lalex06: 65 so, sol, sqod, sqot = getalex06(dl, tl, x, z, 0.0, 0.0, 0.0, 0.0) 66 lgotatm = True 67 else: 68 _stop(3, 'NO MOLECULAR OPACITY TABLE CHOSEN', 69 ' RUN STOPPED. X Z TL=', x, z, tl) 70 71 # Get interior opacity if needed (else go to label 1000) (Fortran line 101). 72 if not (tl < nml.tmolmin and lgotatm): 73 74 if z > 0.1 and tl > 7.0: 75 # Helium-burning region: interpolate towards the pure-Z table 76 # (Fortran lines 109–138). 77 if not nml.lpurez: 78 # Suspected bug (see the page): the message still says 79 # 5 x 10^7 K although the threshold is now log T > 7.0. 80 _stop(1, ' ERROR: Z>0.10 T > 5 X 10^7 K', 81 ' NEED PURE Z TABLE TO CONTINUE. Z,LOG T=', z, tl) 82 oz, olz, qodz, qotz = gtpurz(dl, tl, 0.0, 0.0, 0.0, 0.0) 83 if nml.lopal95: 84 # Interpolate to the maximum Z in the table. 85 zit = 0.1 86 o, ol, qod, qot = getopal95(dl, tl, x, zit, o, ol, qod, qot) 87 else: 88 _stop(3, 'NO OPACITY TABLE CHOSEN', 89 ' RUN STOPPED. X Z TL=', x, z, tl) 90 slope = (ol - olz) / (zit - 1.0) 91 ol = olz + (z - 1.0) * slope 92 o = 10.0 ** ol 93 slope = (qod - qodz) / (zit - 1.0) 94 qod = qodz + (z - 1.0) * slope 95 slope = (qot - qotz) / (zit - 1.0) 96 qot = qotz + (z - 1.0) * slope 97 98 # Z too large, or Z outside a single table (Fortran lines 140–147). 99 # Suspected bug (see the page): with OPAL 1995 off, the second 100 # test turns the "no table chosen" input error into this stop. 101 elif (z > 0.12 or (abs(z - comp_mod.zenv) > nml.optol 102 and not newopac_mod.l2z and not nml.lopal95)): 103 _stop(1, ' Z>0.12 T < 10^7 K', 104 ' OUTSIDE OPAL OPACITY TABLE RANGE OR Z', 105 ' OUTSIDE SINGLE TABLE USED.Z,ZENV,LOG T=', 106 z, comp_mod.zenv, tl) 107 108 # OPAL 1995 (Fortran lines 154–155). 109 elif nml.lopal95: 110 o, ol, qod, qot = getopal95(dl, tl, x, z, o, ol, qod, qot) 111 112 else: 113 # Final trap: should not be able to get here (Fortran lines 162–166). 114 _stop(3, 'NO OPACITY TABLE CHOSEN', 115 ' RUN STOPPED. X Z TL=', x, z, tl) 116 117 # label 1000: ramp between surface and interior opacity (Fortran lines 173–188). 118 if lgotatm and tl <= nml.tmolmax: 119 if tl >= nml.tmolmin: 120 rmpwt = (tl - nml.tmolmin) / (nml.tmolmax - nml.tmolmin) 121 o = rmpwt * o + (1.0 - rmpwt) * so 122 ol = math.log10(o) 123 qod = rmpwt * qod + (1.0 - rmpwt) * sqod 124 qot = rmpwt * qot + (1.0 - rmpwt) * sqot 125 else: 126 o = so 127 ol = sol 128 qod = sqod 129 qot = sqot 130 131 # Conductive opacity correction: Potekhin (Fortran lines 190–195). 132 if nml.lcondopacp: 133 # Potekhin conductive opacity. 134 oc, ocl, qodc, qotc, lcondo = condopacpint( 135 dl, tl, x, z, 0.0, 0.0, 0.0, 0.0, fxion, False) 136 else: 137 lcondo = False 138 139 if not lcondo: 140 # No Potekhin opacity: try Hubbard-Lampe (Fortran lines 197–213). 141 # Question (see the page): this fallback is silent. 142 if tl < 4.2: 143 return o, ol, qod, qot 144 elif dl < 2.0 * tl - 13.0: 145 return o, ol, qod, qot 146 else: 147 cond = (math.log10(1.0 - 0.6 * x) - 14.6196 148 - (3.5853 + 0.1386 * dl) * dl 149 + (5.1324 - 0.3219 * tl) * tl + 0.3901 * dl * tl) 150 oc = 10.0 ** cond 151 qodc = 0.3901 * tl - 0.2772 * dl - 3.5853 152 qotc = 0.3901 * dl - 0.6438 * tl + 5.1324 153 154 # Save the radiative opacity, then add the opacities: 1/O = 1/OX + 1/OC 155 # (Fortran lines 215–222). 156 ox = o 157 qodx = qod 158 qotx = qot 159 o = ox * oc / (ox + oc) 160 ol = math.log10(o) 161 qod = qodx + qodc - (ox * qodx + oc * qodc) / (ox + oc) 162 qot = qotx + qotc - (ox * qotx + oc * qotc) / (ox + oc) 163 164 return o, ol, qod, qot