chkscal.f

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

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

Variables on this page (12)

c4pi, clsun, crsun, csig, endage, itrack, lendag, rescal, sr, xenv0a, xls, xlstol

Purpose

This file calibrates a track of another star (not the Sun) to a given radius and luminosity: it watches the track for the moment it crosses the target radius $R_*$, compares the luminosity there with the target $L_*$, and adjusts the hydrogen abundance $X$ (thus $Y$) for the next run until the track passes through $(R_*, L_*)$. Then it sets up a final run that stops at the age where the track reaches $R_*$.

How it fits

Inputs: the arguments BL ($\log L/L_\odot$), TEFFL ($\log T_\mathrm{eff}$), DAGE (age, Gyr), NK (kind card) and ALRI (the previous model's radius, kept by the main program and started by SETSCAL); in NMLPAR_MOD (set by PARMIN), the targets SR ($R_*/R_\odot$; computed by SETSCAL from STEFF when the option LTEFF gives the star's $T_\mathrm{eff}$ instead, so this routine needs only SR, settled from the code, Nov 9, 2026: setscal.f, lines 102–106 ↗), XLS ($L_*/L_\odot$) and the tolerance XLSTOL; the previous model's BLI, AGEI, local variables kept by SAVE.
Who calls it: the main program (main.f), for even kind cards, when the option LCALST is on and the star has not yet been found (LSTAR false), from the second model on (main.f, lines 1117–1123 ↗; the first model only sets TEFFI, which nothing reads: unused?).
Outputs: the arguments LPASSR (main.f then ends the run), LSTAR and ALRI; BLR, AGER, BLRP are saved locals; for the next kind cards, the new $X$ (RESCAL(2,·) in RUNLOOP_MOD, XENV0A in NMLPAR_MOD), or the stopping age (LENDAG in RUNLOOP_MOD, ENDAGE in NMLPAR_MOD). Lines to the screen and to ITRACK.

The method

The radius from $L$ and $T_\mathrm{eff}$, in solar units: $$\frac{R}{R_\odot} = \frac{\sqrt{L/(4\pi\sigma)}}{T_\mathrm{eff}^2\,R_\odot}.$$ If the track has not crossed $R_*$ since the previous model, store $L$, the age and $R$, and return.

The radius is linear, not a logarithm. The comments call ALR and ALRI "log(R/Rsun)" (lines 81–82 ↗), but the code computes $R/R_\odot$ itself (line 96 ↗) and compares it with SR, which is $R/R_\odot$ (line 79 ↗, line 98 ↗), so the comparison is consistent, the interpolation is linear in $R$, and the comments are stale.

Question 1 (lines 97–117 ↗). The first crossing only. The comment "Check if track has passed through Teff the right number of times" has no code behind it: The first crossing of $R_*$ (in either direction) is used. For tracks that cross $R_*$ more than once (pre-main sequence, then main sequence), is the first crossing the right one? 2. At the crossing, interpolate linearly in $R$ between the two models for $\log L$ and the age at $R_*$. 3. If $|L/L_\odot - L_*/L_\odot| \le$ XLSTOL, the star is found: the next two kind cards run with the same $X$ and stop at that age (in years, $10^9 \times$ Gyr). The test is in linear units, so XLSTOL is in $L_\odot$, while the secant step below works in $\log L$ (settled from the code, Nov 9, 2026: line 128 ↗, line 168 ↗). 4. Otherwise, at the first pass (NK = 2) lower $X$ by 0.01 ("Y=Y+0.01"); afterwards take a secant step in $X$ against $\log L$: $$X_\mathrm{new} = X + \frac{X - X_\mathrm{prev}}{\log L_R - \log L_{R,\mathrm{prev}}}\,\left(\log L_* - \log L_R\right).$$

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 false crossing at the start of a run that begins above $R_*$.

(lines 97–103 ↗; setscal.f, line 107 ↗; main.f, lines 1119–1123 ↗) ALRI starts at 0 (setscal.f, line 107 ↗) and is not updated at a crossing (lines 118–178 ↗ leave it as it was), and main.f's first call of each run is at its second model (the first only sets TEFFI). If that model's radius is above $R_*$ (a pre-main-sequence start with $R_*$ on the main sequence, say), the test ALR.GT.SR.AND.ALRI.LE.SR (line 98 ↗) holds at once: LPASSR ends the run at its second model, and $L$ and the age "at $R_*$" are interpolated between this model and BLI, AGEI, which are unset in the first run (saved locals with no DATA). Since ALRI stays 0, every later run stops the same way, and the secant step works on these values. When the track starts below $R_*$, the first comparison of each later run is made against the previous run's last stored model rather than the run's own, which is harmless as long as every run starts on the same side of $R_*$. Only with LCALST (off by default, nmlpar_mod.f, line 160 ↗) and a track that starts above $R_*$; a track that starts below $R_*$ and expands through it is not affected.

A possible fix (in main.f, then here):

main.f (another file)

@@ lines 1116–1123 @@
 C TEST IF MODEL IS NEAR DESIRED Teff AND L. IF NOT RESCALE AND TRY AGAIN.
          IF (LCALST .AND. .NOT. LSTAR) THEN
             IF(MOD(NK,2).EQ.0)THEN
+C 11/09/26 KC via Claude
+C code changed to mark the first model of each run, so that CHKSCAL
+C starts its search from the run's own models.
              IF(MODELN.EQ.1) THEN
                 TEFFI = 10.0D0**TEFFL
+                ALRI = -1.0D0
              ELSE
                 CALL CHKSCAL(BL, TEFFL, DAGE, NK,ALRI,LSTAR,LPASSR)
                 IF(LPASSR) GOTO 200

chkscal.f

@@ lines 93–99 @@
 C     If not store L and age and return.
       LPASSR=.FALSE.
       TEFF = 10.0D0**TEFFL
+C 11/09/26 KC via Claude
+C code changed to store the first model of a run without testing it.
       ALR = SQRT((10.0D0**BL)*CLSUN/(C4PI*CSIG))/(TEFF*TEFF*CRSUN)
+      IF(ALRI.LT.0.0D0) THEN
+         BLI = BL
+         AGEI = DAGE
+         ALRI = ALR
+         RETURN
+      ENDIF
       IF(ALR.GT.ALRI) THEN
          IF(.NOT.(ALR.GT.SR.AND.ALRI.LE.SR)) THEN
               BLI = BL

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged). No run reached CHKSCAL in the 34 cases run: LCALST is off by default (nmlpar_mod.f, line 160 ↗) and in all of YREC's cases, and no user document offers it (docs/archive/F90_CONVERSION_PLAN.md, line 247: "no case sets LCALST"). What would reach it: LCALST on, with a track whose second model lies above $R_*$, such as a pre-main-sequence start with $R_*$ on the main sequence. Intent: none found.

Definitions

subroutineCHKSCAL

lines 6–179 · GitHub

Calibrates $X$ (thus $Y$) so that a track passes through a target $(R_*, L_*)$: detects the crossing of $R_*$, interpolates $L$ and the age there, and sets $X$ by a secant step, or sets the final run's stopping age.

Uses modules CONSTANTS_MOD, IOUNITS_MOD, NMLPAR_MOD, RUNLOOP_MOD

Called by (file level)

Change log

The changes since rp_1

Changed since rp_1. /CALSTAR/ retired: ALRI, LSTAR, LPASSR are now arguments (owned by the main program), BLI, AGEI, BLR, AGER, BLRP saved locals, the targets in NMLPAR_MOD; COMMON blocks to modules, IMPLICIT NONE.

1C 2C 3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 4C CHKSCAL 5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ 6 SUBROUTINE CHKSCAL(BL, TEFFL, DAGE, NK,ALRI,LSTAR,LPASSR) 7C 10/14/26 MHP via Claude 8C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 9 USE NMLPAR_MOD 10C 10/14/26 MHP via Claude 11C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD. 12 USE CONSTANTS_MOD 13C 10/14/26 MHP via Claude 14C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 15 USE IOUNITS_MOD 16C 10/14/26 MHP via Claude 17C code changed to convert COMMON/RUNLOOP/ into RUNLOOP_MOD. A blanket 18C USE imports exactly the names the COMMON declaration made global here, 19C so it cannot change meaning. See runloop_mod.f for the block history 20C and for why the module declares the types the implicit rules supplied. 21 USE RUNLOOP_MOD 22C 10/20/26 MHP via Claude 23C code changed to replace implicit typing with IMPLICIT NONE and 24C explicit declarations. 25 IMPLICIT NONE 26 REAL*8 AGEI,AGER,ALR,ALRI,BL,BLI,BLR,BLRP,DADR,DAGE,DLDR,DXDL, 27 * * TEFF,TEFFL,X,XP 28 INTEGER NK 29 LOGICAL*4 LPASSR,LSTAR 30C ONLY CALLED FOR EVEN NK, ASSUMES RESCALING ON ODD NK AND EVOLVING 31C ON EVEN NK 32C 10/13/26 MHP via Claude 33C code changed to replace COMMON/CALSTAR/ with explicit arguments. 34C setscal.f initializes ALRI, LSTAR and LPASSR, chkscal.f reads and 35C updates ALRI and sets the two flags, and main.f reads the flags; main.f 36C calls both routines, so it now owns all three. ALRI is read as well as 37C written inside chkscal.f and is carried across calls there, which 38C main.f's blanket SAVE preserves exactly as the block did. 39C 9/25/26 MHP via Claude 40C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 41C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 42C 10/12/26 MHP via Claude 43C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 44C namelist variables that the physics routines read from 27 COMMON 45C blocks which also held derived or working members. Those blocks keep 46C only their non-namelist members; five of them -- NEWMX, VMULT, 47C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 48C Derivation-only namelist variables from the same blocks went to 49C COMMON/NMLBYPASS/ instead, not here. 50C 51C Declared identically in all 121 files that hold any of these members. 52C Generated from one member list, so every copy is byte-identical; the 53C linker's 'size of symbol changed' warning is the check that they are. 54C MHP 10/24 ADDED STOP CRITERIA FOR CENTRAL H,D,AND HE4 55C 10/12/26 MHP via Claude 56C code changed to demote BLI, AGER, BLR, BLRP and AGEI from 57C COMMON/CALSTAR/ to locals of this routine. All five are written and 58C read only here -- BLI/AGEI at the interpolation setup, BLR/AGER/BLRP in 59C the radius-matching iteration below -- while main.f, parmin.f and 60C setscal.f merely carried them positionally. The blanket SAVE below 61C keeps their values across calls exactly as COMMON storage did, which 62C BLRP relies on, being set on one call and read on the next. Only 63C ALRI, LSTAR and LPASSR are genuinely shared; setscal.f initializes 64C those three. 65C 9/25/26 MHP via Claude 66C code changed to consolidate all I/O logical unit numbers, previously 67C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 68C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 69C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 70C COMMON/IOUNITS/ 71 SAVE 72C LSTAR T - have got a star at Teff and L 73C LPASSR T - on run have just passed Teff 74C XLS Luminosity (L/Lsun) wanted by adjusting Y 75C XLSTOL tolerance wanted for luminosity 76C LTEFF T - specify Teff for star 77C F - specify R/Rsun for star 78C STEFF Effective temperature of star (K) or... 79C SR Radius of star (R/Rsun) 80C TEFF Teff of current model 81C ALR log(R/Rsun) of current model 82C ALRI log(R/Rsun) of previous model 83C DAGE age of current model (Gyr) 84C AGEI age of previous model (Gyr) 85C AGER age of model at R* 86C BL luminosity of current model log(L/Lsun) 87C BLI luminosity of previous model 88C BLR luminosity of model at R 89C BLRP luminosity of model at R* of previous run 90C XP X of previous run = RESCAL(2, NK-1) 91C 92C Check if star has passed R*. 93C If not store L and age and return. 94 LPASSR=.FALSE. 95 TEFF = 10.0D0**TEFFL 96 ALR = SQRT((10.0D0**BL)*CLSUN/(C4PI*CSIG))/(TEFF*TEFF*CRSUN) 97 IF(ALR.GT.ALRI) THEN 98 IF(.NOT.(ALR.GT.SR.AND.ALRI.LE.SR)) THEN 99 BLI = BL 100 AGEI = DAGE 101 ALRI = ALR 102 RETURN 103 END IF 104 ELSE 105 IF (.NOT.(ALR.LT.SR.AND. ALRI.GE.SR)) THEN 106 BLI = BL 107 AGEI = DAGE 108 ALRI = ALR 109 RETURN 110 END IF 111 ENDIF 112C 113C Check if track has passed through Teff the right number of 114C times. If not store L and age and return. 115C ZZZ 116 WRITE(*,*) ' Just passed R*' 117 WRITE(ITRACK,*) '#Just passed R*' 118 LPASSR = .TRUE. 119C 120C Have previous L,Age and current L,Age (one before R* and 121C one after R*). Interpolate to get L,Age at R* 122 DLDR = (BL-BLI)/(ALR-ALRI) 123 BLR = BL + DLDR*(SR-ALR) 124 DADR = (DAGE-AGEI)/(ALR-ALRI) 125 AGER = DAGE + DADR*(SR-ALR) 126 WRITE(*,*) ' X, LogL/Lsun at R* =', RESCAL(2,NK-1), BLR 127 WRITE(ITRACK,*) '#X, LogL/Lsun at R* =', RESCAL(2,NK-1), BLR 128 IF (ABS(10.0D0**BLR-XLS) .LE. XLSTOL) THEN 129C Get here then have track that passes through specified 130C L and R. Use age at R for final 131C run to stop at that age. Do one more run 132C stopping at interpolated age. 133 LSTAR=.TRUE. 134 LENDAG(NK+1) = .TRUE. 135 ENDAGE(NK+1) = AGER*1.0D9 136 LENDAG(NK+2) = .TRUE. 137 ENDAGE(NK+2) = AGER*1.0D9 138 RESCAL(2,NK+1) = RESCAL(2,NK-1) 139 XENV0A(NK+1) = RESCAL(2,NK+1) 140 XENV0A(NK+2) = XENV0A(NK+1) 141C 142 WRITE(*,*) ' Have hit R* & L*, prepare final run to age:', 143 1 1 AGER 144 WRITE(ITRACK,*)'#Have hit R* & L*, prepare final run to age:', 145 1 1 AGER 146 RETURN 147 ELSE 148 IF (NK .EQ. 2) THEN 149C First time through. Save L and X at R*. 150C Add 0.01 to Y. Start next run. 151 BLRP = BLR 152 XP = RESCAL(2,NK-1) 153 X = XP - 0.01D0 154 RESCAL(2,NK+1) = X 155 XENV0A(NK+1) = X 156 XENV0A(NK+2) = X 157C 158 WRITE(*,*) ' NK=2, Y=Y+0.01, Setup next run, X=', X 159 WRITE(ITRACK,*) '#NK=2, Y=Y+0.01, Setup next run, X=', X 160 RETURN 161 ELSE 162C If NK=4,6,8,... (second and more times through) then 163C Use current and previous values of L at R and X to calculate 164C dX/dlogL. Save L. Start next run. 165 X = RESCAL(2, NK-1) 166 XP = RESCAL(2, NK-3) 167 DXDL = (X-XP)/(BLR-BLRP) 168 X = DXDL*(LOG10(XLS)-BLR)+X 169C 170 WRITE(*,*) ' Setup next run, NK, X =', NK, X 171 WRITE(ITRACK, *) ' Setup next run, NK, X =', NK, X 172 BLRP = BLR 173 RESCAL(2,NK+1) = X 174 XENV0A(NK+1) = X 175 XENV0A(NK+2) = X 176 END IF 177 ENDIF 178 RETURN 179 END

1"""Python version of chkscal.f (YREC 7be4133): calibrating a star's track 2to a target radius and luminosity.""" 3 4import math 5 6from pyrec import constants_mod 7from pyrec import iounits_mod 8from pyrec import nmlpar_mod 9from pyrec import runloop_mod 10 11# CHKSCAL has SAVE (formerly part of COMMON/CALSTAR/): the previous model's 12# L and age, and L and age at R*; BLRP is set in one run and read in the 13# next. 14bli = 0.0 15agei = 0.0 16blr = 0.0 17ager = 0.0 18blrp = 0.0 19 20 21def chkscal(bl: float, teffl: float, dage: float, nk: int, alri: float, lstar: bool, 22 lpassr: bool) -> tuple[float, bool, bool]: 23 """Calibrate X (thus Y) so that a track passes through a target 24 (R*, L*): detect the crossing of R*, interpolate L and the age there, 25 and set X for the next run by a secant step, or set the final run's 26 stopping age. 27 28 Translates: chkscal.f#CHKSCAL (7be4133) 29 30 Arguments: bl, log L/L_sun; teffl, log Teff; dage, the age (Gyr); nk, 31 the kind card; alri, the previous model's R/R_sun (despite the comments, 32 not a log); lstar, the star found; lpassr, R* just passed. 33 34 Returns (alri, lstar, lpassr). The new X goes to runloop_mod.rescal and 35 nmlpar_mod.xenv0a; the stopping age to runloop_mod.lendag and 36 nmlpar_mod.endage. 37 38 """ 39 40 global bli, agei, blr, ager, blrp 41 c = constants_mod 42 nml = nmlpar_mod 43 rescal = runloop_mod.rescal 44 itrack = iounits_mod.itrack 45 46 # Has the star passed R*? If not, store L and age and return (Fortran lines 94–111). 47 # Question (see the page): linear or log radius. 48 # Suspected bug (see the page): alri starts at 0 and is not reset at a 49 # crossing, so a track that starts above R* crosses at once. 50 lpassr = False 51 teff = 10.0**teffl 52 alr = math.sqrt((10.0**bl) * c.clsun / (c.c4pi * c.csig)) / ( 53 teff * teff * c.crsun) 54 if alr > alri: 55 if not (alr > nml.sr and alri <= nml.sr): 56 bli = bl 57 agei = dage 58 alri = alr 59 return alri, lstar, lpassr 60 else: 61 if not (alr < nml.sr and alri >= nml.sr): 62 bli = bl 63 agei = dage 64 alri = alr 65 return alri, lstar, lpassr 66 67 # Passed (Fortran lines 112–118). 68 # Question (see the page): the first crossing only. 69 print(" Just passed R*") 70 iounits_mod.write(itrack, " #Just passed R*") 71 lpassr = True 72 73 # Interpolate L and the age at R*, linearly in R (Fortran lines 120–127). 74 dldr = (bl - bli) / (alr - alri) 75 blr = bl + dldr * (nml.sr - alr) 76 dadr = (dage - agei) / (alr - alri) 77 ager = dage + dadr * (nml.sr - alr) 78 print(f" X, LogL/Lsun at R* = {rescal[2, nk - 1]} {blr}") 79 iounits_mod.write(itrack, f" #X, LogL/Lsun at R* = {rescal[2, nk - 1]}" 80 f" {blr}") 81 82 # Found: one more run (Fortran lines 128–146). 83 if abs(10.0**blr - nml.xls) <= nml.xlstol: 84 # Found: one more run, with the same X, stopping at that age. 85 lstar = True 86 runloop_mod.lendag[nk + 1] = True 87 nml.endage[nk + 1] = ager * 1.0e9 88 runloop_mod.lendag[nk + 2] = True 89 nml.endage[nk + 2] = ager * 1.0e9 90 rescal[2, nk + 1] = rescal[2, nk - 1] 91 nml.xenv0a[nk + 1] = rescal[2, nk + 1] 92 nml.xenv0a[nk + 2] = nml.xenv0a[nk + 1] 93 print(f" Have hit R* & L*, prepare final run to age: {ager}") 94 iounits_mod.write(itrack, " #Have hit R* & L*, prepare final run " 95 f"to age: {ager}") 96 return alri, lstar, lpassr 97 98 # First time through: save L at R*, lower X by 0.01 (Fortran lines 148–160). 99 if nk == 2: 100 blrp = blr 101 xp = rescal[2, nk - 1] 102 x = xp - 0.01 103 rescal[2, nk + 1] = x 104 nml.xenv0a[nk + 1] = x 105 nml.xenv0a[nk + 2] = x 106 print(f" NK=2, Y=Y+0.01, Setup next run, X= {x}") 107 iounits_mod.write(itrack, f" #NK=2, Y=Y+0.01, Setup next run, X= {x}") 108 return alri, lstar, lpassr 109 110 # NK = 4, 6, 8, ...: a secant step in X against log L (Fortran lines 161–176). 111 x = rescal[2, nk - 1] 112 xp = rescal[2, nk - 3] 113 dxdl = (x - xp) / (blr - blrp) 114 x = dxdl * (math.log10(nml.xls) - blr) + x 115 print(f" Setup next run, NK, X = {nk} {x}") 116 iounits_mod.write(itrack, f" Setup next run, NK, X = {nk} {x}") 117 blrp = blr 118 rescal[2, nk + 1] = x 119 nml.xenv0a[nk + 1] = x 120 nml.xenv0a[nk + 2] = x 121 return alri, lstar, lpassr