chkcal.f

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

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

Variables on this page (10)

calsolzx, cmixla, iowr, lcalsolzx, rescal, toll, tolr, tolz, xenv0a, zenv0a

Purpose

This file drives the solar calibration: after each solar model, it checks whether the model has the solar luminosity and radius (and, optionally, the solar surface $Z/X$) at the solar age; if not, it computes new guesses for the hydrogen abundance $X$ and the mixing length $\alpha$ (and $Z$) for the next model, by a Newton step.

How it fits

Inputs: the arguments BL ($\log L/L_\odot$ of the model), RL ($\log R/R_\odot$, computed by the caller from $L$ and $T_\mathrm{eff}$), N (the current kind card) and SOLZX (the model's surface $Z/X$); the tolerances TOLL, TOLR, TOLZ, the switch LCALSOLZX and the target CALSOLZX in NMLPAR_MOD, set in PARMIN and SETCAL.
Who calls it: the main program (main.f), when solar calibration is on (LCALS), at every third kind card (MOD(NK,3).EQ.0, main.f, lines 1147–1152 ↗): "solar models in 3 kind cards" (MHP 5/96).

The header predates the three-card cycle. It says "ONLY CALLED FOR EVEN NK, ASSUMES RESCALING ON ODD NK AND EVOLVING ON EVEN NK" (lines 22–23 ↗), but main.f calls it when NK is a multiple of 3 (main.f, line 1147 ↗), which matches the routine's use of N-2 and N+1..N+3 (lines 93–97 ↗; lines 153–163 ↗).
Outputs: the arguments LSOL and the derivatives DLDX, DRDX, DLDA, DRDA, which the main program keeps (they were in /CALSUN/, also used by WRTMONTE, removed since rp_1); the new $X$ (RESCAL(2,N+1), and XENV0A for cards $N+1$ to $N+3$), $\alpha$ (CMIXLA, cards $N+1$ to $N+3$) and $Z$ (RESCAL(3,N+1), ZENV0A) in RUNLOOP_MOD (RESCAL) and NMLPAR_MOD (the others); lines to unit IOWR.
Kept between calls: the previous model's BLP, RLP and the steps DX, DA, now local variables kept by SAVE.

The method

Converged if $|\log L/L_\odot| <$ TOLL and $|\log R/R_\odot| <$ TOLR and, with the $Z/X$ option, $|\log_{10}(Z/X) - \log_{10}(Z/X)_\odot| <$ TOLZ: set LSOL and return. The target is thus exactly $L_\odot$ and $R_\odot$. These are the code's CLSUN and CRSUN (in CONSTANTS_MOD, settable in the PHYSICS namelist), $3.8515\times10^{33}$ erg s⁻¹ and $6.9598\times10^{10}$ cm by default; BL and RL are measured in them (how these compare with other adopted solar values is the question below).

Question 1 (line 59 ↗; constants_mod.f, lines 44–49 ↗). The solar values targeted. Convergence means $|\log L/L_\odot| <$ TOLL with the default $L_\odot$ = CLSUN = $3.8515\times10^{33}$ erg s⁻¹ and $R_\odot$ = CRSUN = $6.9598\times10^{10}$ cm, while the IAU 2015 nominal values are $3.828\times10^{33}$ erg s⁻¹ and $6.957\times10^{10}$ cm. With TOLL = $10^{-5}$, the 0.6% difference in $L_\odot$ ($0.0027$ in $\log L$) is far larger than the tolerance, so the calibrated $X$ and $\alpha$ depend on which value is used. Are the defaults kept for continuity with published YREC solar models, and should the documentation of TOLL, CLSUN and CRSUN say that a calibration matches whatever $L_\odot$, $R_\odot$ the namelist sets?

Question 2 (lines 59–72 ↗). $Z/X$ not converged. If $L$ and $R$ are within tolerance but $Z/X$ is not, the routine falls through to the Newton step and also changes $X$ and $\alpha$ slightly. Is this intended? (A target CALSOLZX of 0 is bug 2.)
Newton step with fixed partial derivatives ("empirical," LLP 6/28/09): $\partial\log L/\partial X = -3.78$, $\partial\log R/\partial X = -0.89$, $\partial\log L/\partial\alpha = 0.0139$, $\partial\log R/\partial\alpha = -0.050$. It solves $$\log L + L_X\,\Delta X + L_\alpha\,\Delta\alpha = 0, \qquad \log R + R_X\,\Delta X + R_\alpha\,\Delta\alpha = 0,$$ that is $$\Delta\alpha = \frac{\log L\,R_X/L_X - \log R}{R_\alpha - L_\alpha R_X/L_X}, \qquad \Delta X = -\frac{\log L + L_\alpha\,\Delta\alpha}{L_X}.$$

Question 3 (lines 74–84 ↗; main.f, line 1163 ↗). Fixed derivatives. Why fixed "empirical" derivatives rather than secants from the previous cycles (BLP, RLP are already saved)? A user changing the physics (another opacity, say) would want to know whether convergence slows, or fails within the 15 cycles main.f allows, when the true derivatives differ from these.
The new values go to kind cards $N+1$ to $N+3$, starting from those of card $N-2$ (the start of the current cycle). With the $Z/X$ option, $Z$ is scaled: $Z_\mathrm{new} = Z_\mathrm{old}\,(Z/X)_\odot / (Z/X)_\mathrm{model}$.

Design

An older method, which measured the derivatives with extra runs at kind cards 3, 6 and 9, is still in the file but skipped by GOTO 1234 (lines 85–146 ↗; unused?; on that path DX and DA are never set, their assignments being commented out, line 91 ↗, line 110 ↗). The local ZLP is set (line 181 ↗) and never read (unused?).

From YREC's documentation

Control Options (Calibration): TOLL, TOLR are the tolerances of the solar calibration in $L$ and $R$; LCALSOLZX enforces the mixture's $Z/X$, CALSOLZX (0.02292 for GS98), and without it $Z$ is held constant in the calibration loop.

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. (wrong output) Checked by KC (Nov 10, 2026) The first cycle prints an "Old BL" and "Old RL" that were never set.

(lines 167–170 ↗; lines 179–180 ↗) BLP and RLP are saved locals with no DATA and are set only at the end of the Newton step, so on the first call (N = 3) the lines "New BL, Old BL, Delta BL" and "New RL, Old RL, Delta RL" print whatever they hold before any assignment (zero under most compilers, not promised by the standard). Only the output to IOWR is affected; the Newton step does not use them. Only with LCALS (off by default).

A possible fix:

chkcal.f

@@ lines 164–173 @@
          XENV0A(N+1) = RESCAL(2,N+1)
          XENV0A(N+2) = RESCAL(2,N+1)
          XENV0A(N+3) = RESCAL(2,N+1)
-         WRITE(IOWR,*) "New BL, Old BL, Delta BL: ",
-     *       BL, BLP, BL-BLP
-         WRITE(IOWR,*) "New RL, Old RL, Delta RL: ",
-     *       RL, RLP, RL-RLP
+C 11/09/26 KC via Claude
+C code changed to print the previous model's L and R only once there
+C is one.
+         IF(N.GT.3)THEN
+            WRITE(IOWR,*) "New BL, Old BL, Delta BL: ",
+     *          BL, BLP, BL-BLP
+            WRITE(IOWR,*) "New RL, Old RL, Delta RL: ",
+     *          RL, RLP, RL-RLP
+         ENDIF
          WRITE(IOWR,*) "New X, Old X, DX: ",
      *       RESCAL(2,N+1), RESCAL(2,N-2), DX
          WRITE(IOWR,*) "New A, Old A, DA: ",

Checked by running YREC (Nov 9, 2026; a line added at each suspected path counted how often it was taken, the outputs otherwise unchanged). YREC's solar calibration (examples/run_standard_solar_model/Test_solarcal_noGS_norot, the only case with LCALS on) printed these lines twice, the first time on the first call ($N$ = 3, with $\log L$ and $\log R$ off by about $1.4\times10^{-4}$ and $2.2\times10^{-5}$ dex), before BLP and RLP had been set: The suspected path is taken in every calibration that needs more than one cycle. No fix was run, since it changes only two printed lines and nothing computed reads them. Intent: none found; the comment at lines 42–45 ↗ explains why BLP and RLP are saved across calls, not what they hold on the first. Only the printed output is wrong.

2. (never occurred) A target $Z/X$ of 0 sets $Z = 0$.

(line 61 ↗; lines 156–157 ↗) The convergence test uses $Z/X$ only when LCALSOLZX is true and CALSOLZX > 0; the update of $Z$ tests only LCALSOLZX. With LCALSOLZX true and CALSOLZX = 0 (or negative), the test ignores $Z/X$ while the update sets RESCAL(3,N+1) = RESCAL(3,N-2)*CALSOLZX/SOLZX = 0, so the next three runs start with no metals. Not on the default path: LCALSOLZX is false by default and CALSOLZX is 0.02292 (nmlpar_mod.f, lines 53–54 ↗); a target of 0 is not a sensible input, so this is fragile code rather than a wrong result in practice.

A possible fix:

chkcal.f

@@ lines 153–159 @@
          CMIXLA(N+1) = CMIXLA(N-2)+DA
          CMIXLA(N+2) = CMIXLA(N+1)
          CMIXLA(N+3) = CMIXLA(N+1)
-         IF(LCALSOLZX)THEN
+C 11/09/26 KC via Claude
+C code changed to scale Z only under the test the convergence check uses.
+         IF(LCALSOLZX .AND. CALSOLZX.GT.0.0D0)THEN
             RESCAL(3,N+1) = RESCAL(3,N-2)*CALSOLZX/SOLZX
             DZ = RESCAL(3,N+1) - RESCAL(3,N-2)
             ZENV0A(N+1) = RESCAL(3,N+1)
@@ lines 172–178 @@
      *       RESCAL(2,N+1), RESCAL(2,N-2), DX
          WRITE(IOWR,*) "New A, Old A, DA: ",
      *       CMIXLA(N+1), CMIXLA(N-2), DA
-         IF(LCALSOLZX)THEN
+         IF(LCALSOLZX .AND. CALSOLZX.GT.0.0D0)THEN
             WRITE(IOWR,*) "New Z, Old Z, DZ: ",
      *       RESCAL(3,N+1), RESCAL(3,N-2), DZ
          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). The same calibration took the path twice, with CALSOLZX = 0.02292 both times (the model's $Z/X$ between 0.0229188 and 0.02292); the target was never 0 or below, and every one of YREC's namelists that sets it gives 0.02292, the GS98 value. The fix was not run, since no case reaches the condition. What would reach it: LCALS and LCALSOLZX on with CALSOLZX $\le 0$. Intent: none found.

Definitions

subroutineCHKCAL

lines 2–186 · GitHub

Checks a solar model against $L_\odot$, $R_\odot$ (and $Z/X$); if not calibrated, sets new $X$, $\alpha$ (and $Z$) for the next three kind cards by a Newton step with fixed derivatives.

Uses modules IOUNITS_MOD, NMLPAR_MOD, RUNLOOP_MOD

Called by (file level)

Change log

The changes since rp_1

Changed since rp_1. /CALSUN/ retired: the derivatives and LSOL are now arguments, and BLP, RLP, DX, DA saved locals; COMMON blocks to modules, IMPLICIT NONE.

Changed since a8a5950. Only the comment on DLDA changed: It read "+0.139," ten times the value used, and now reads "+0.0139"; a fit over the iterates of a recalibration of 11 solar decks gives $\partial\log L/\partial\alpha$ of order 0.01, so 0.0139 is the value meant (this settles former question 1).

1C $$$$$$ 2 SUBROUTINE CHKCAL(BL,RL,N,SOLZX,DLDX,DRDX,DLDA,DRDA,LSOL) 3C 10/14/26 MHP via Claude 4C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD. 5 USE NMLPAR_MOD 6C 10/14/26 MHP via Claude 7C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD. 8 USE IOUNITS_MOD 9C 10/14/26 MHP via Claude 10C code changed to convert COMMON/RUNLOOP/ into RUNLOOP_MOD. A blanket 11C USE imports exactly the names the COMMON declaration made global here, 12C so it cannot change meaning. See runloop_mod.f for the block history 13C and for why the module declares the types the implicit rules supplied. 14 USE RUNLOOP_MOD 15C 10/20/26 MHP via Claude 16C code changed to replace implicit typing with IMPLICIT NONE and 17C explicit declarations. 18 IMPLICIT NONE 19 REAL*8 BL,BLP,DA,DLDA,DLDX,DRDA,DRDX,DX,DZ,RL,RLP,SOLZX,ZL,ZLP 20 INTEGER N 21 LOGICAL*4 LSOL 22C ONLY CALLED FOR EVEN NK, ASSUMES RESCALING ON ODD NK AND EVOLVING 23C ON EVEN NK 24C 10/12/26 MHP via Claude 25C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 26C namelist variables that the physics routines read from 27 COMMON 27C blocks which also held derived or working members. Those blocks keep 28C only their non-namelist members; five of them -- NEWMX, VMULT, 29C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 30C Derivation-only namelist variables from the same blocks went to 31C COMMON/NMLBYPASS/ instead, not here. 32C 33C Declared identically in all 121 files that hold any of these members. 34C Generated from one member list, so every copy is byte-identical; the 35C linker's 'size of symbol changed' warning is the check that they are. 36C MHP 10/24 ADDED STOP CRITERIA FOR CENTRAL H,D,AND HE4 37C 10/13/26 MHP via Claude: removed the unused COMMON/SETT/ declaration 38C (no member referenced here; residue of the obsolete-option removals). 39C 10/13/26 MHP via Claude 40C code changed to replace COMMON/CALSUN/ with explicit arguments. This 41C routine wrote all nine members; only the four derivatives and LSOL are 42C read outside it, by main.f, so those are now arguments. BLP, RLP, DX 43C and DA never left this routine and are now locals, kept across calls by 44C the blanket SAVE below -- which they need, since BLP and RLP are stored 45C on one call and differenced against on the next. 46C MHP 6/13 ADD OPTION TO CALIBRATE SOLAR Z/X, SOLAR Z/X, SOLAR AGE 47C COMMON/CALS2/TOLL,TOLR,LCALS 48C MHP 10/02 MISSING LOGICAL UNIT IOWR ADDED 49C 9/25/26 MHP via Claude 50C code changed to consolidate all I/O logical unit numbers, previously 51C scattered across COMMON/LUOUT/, COMMON/LUNUM/, COMMON/ALEX06/, 52C COMMON/ACDPTH/, COMMON/NEWOPAC/, COMMON/LOPAL95/, COMMON/ATMOS2/, 53C COMMON/ALATM03/, COMMON/OPALEOS/, and COMMON/SCV2/, into a single 54C COMMON/IOUNITS/ 55 SAVE 56C DATA TOLL,TOLR/1.0D-5,1.0D-4/ 57C CHECK TO SEE IF THE MODEL IS A CALIBRATED SOLAR MODEL. 58C IF NOT, ESTIMATE CORRECTIONS TO X AND ALPHA. 59 IF(ABS(BL).LT.TOLL .AND. ABS(RL).LT.TOLR)THEN 60C ADD CHECK FOR Z/X 61 IF(LCALSOLZX .AND.CALSOLZX.GT.0.0D0)THEN 62 ZL = LOG10(SOLZX)-LOG10(CALSOLZX) 63 IF(ABS(ZL).LT.TOLZ)THEN 64 LSOL = .TRUE. 65 GOTO 9999 66 ENDIF 67 ELSE 68C CALIBRATED SOLAR MODEL. SET UP OUTPUT FLAGS AND EXIT 69 LSOL = .TRUE. 70 GOTO 9999 71 ENDIF 72 ENDIF 73 74c Added code to use pre-determined partial derivatives LLP 6/28/09 75 DLDX = -3.78D0 ! empirical result: -3.783 RMS error .070 76 DRDX = -0.89D0 ! empirical result: -0.890 RMS error .048 77C 10/29/26 MHP via Claude 78C DLDA's comment read +0.139, ten times the value used. 0.0139 is right: 79C fitting log L against X, alpha and ln Z over the iterates of the 80C 10/28/26 recalibration of 11 solar decks gives dlogL/dalpha of order 81C 0.01 (median 0.003) and dlogR/dalpha about -0.05, while 0.139 would 82C have moved log L 50 times above the noise. Only the comment changed. 83 DLDA = 0.0139D0 ! empirical result: +0.0139 RMS error .0022 84 DRDA = -0.050D0 ! empirical result: -0.0504 RMS error .0059 85 GOTO 1234 ! Bypass partial derivative code 86c mhp 5/96 added change to compute solar calibration for 3 kind cards 87 IF(N.EQ.3)THEN 88C SET UP RUN TO DETERMINE DERIVATIVE OF L AND R WITH RESPECT TO X. 89C 90C 3.7 is empirical average dL/dX, so trial DX is BL / (dL/dX) ! llp 6/18/09 91c DX = BL/3.7D0 92C DX = +.01 93 RESCAL(2,N+1) = RESCAL(2,1)+DX 94C STORE PREVIOUS L AND R. 95 XENV0A(N+1) = RESCAL(2,N+1) 96 XENV0A(N+2)=XENV0A(N+1) 97 XENV0A(N+3)=XENV0A(N+1) 98 BLP = BL 99 RLP = RL 100 GOTO 9999 101 ELSE IF(N.EQ.6)THEN 102C EVALUATE DERIVATIVE OF L AND R WITH RESPECT TO X. 103c DLDX = (BL - BLP)/DX 104c DRDX = (RL - RLP)/DX 105c WRITE(IOWR,*) "DX,DLDX,BL,BLP,DRDX,RL,RLP: ", 106c * DX,DLDX,BL,BLP,DRDX,RL,RLP 107C SET UP RUN TO DETERMINE DERIVATIVE OF L AND R WITH RESPECT TO ALPHA. 108 109C .042 is typical average dR/dA, so trial DA is RL / (dR/dA) ! llp 6/18/09 110c DA = RL/0.042D0 111 CMIXLA(N+1) = CMIXLA(1)+DA 112 CMIXLA(N+2) = CMIXLA(N+1) 113 CMIXLA(N+3) = CMIXLA(N+1) 114 RESCAL(2,N+1) = RESCAL(2,1) 115 XENV0A(N+1)=XENV0A(1) 116 XENV0A(N+2)=XENV0A(N+1) 117 XENV0A(N+3)=XENV0A(N+1) 118 BLP = BL 119 RLP = RL 120 GOTO 9999 121 ELSE IF(N.EQ.9)THEN 122C EVALUATE DERIVATIVE OF L AND R WITH RESPECT TO ALPHA. 123 DLDA = (BL - BLP)/DA 124 DRDA = (RL - RLP)/DA 125 WRITE(IOWR,*) "DA,DLDA,BL,BLP,DRDA,RL,RLP: ", 126 * * DA,DLDA,BL,BLP,DRDA,RL,RLP 127C USE DERIVATIVES OF L AND R WITH RESPECT TO X AND ALPHA TO 128C GET IMPROVED GUESSES FOR ALPHA AND X. 129 DA = ((BL*DRDX/DLDX-RL)/(DRDA-DLDA*DRDX/DLDX)) 130 DX = -(BL + DLDA*DA)/DLDX 131 CMIXLA(N+1) = CMIXLA(1)+DA 132 CMIXLA(N+2) = CMIXLA(N+1) 133 CMIXLA(N+3) = CMIXLA(N+1) 134 RESCAL(2,N+1) = RESCAL(2,1)+DX 135 XENV0A(N+1) = RESCAL(2,N+1) 136 XENV0A(N+2) = RESCAL(2,N+1) 137 XENV0A(N+3) = RESCAL(2,N+1) 138 BLP = BL 139 RLP = RL 140 WRITE(IOWR,*) "New X, Old X, Calc DX: ", 141 * * RESCAL(2,N+1), RESCAL(2,1), DX 142 WRITE(IOWR,*) "New A, Old A, Calc DA: ", 143 * * CMIXLA(N+1), CMIXLA(1), DA 144 GOTO 9999 145c ELSE 146 endif ! terrminate old partial derivative code 147 148 1234 1234 Continue 149C USE DERIVATIVES OF L AND R WITH RESPECT TO X AND ALPHA TO 150C GET IMPROVED GUESSES FOR ALPHA AND X. 151 DA = ((BL*DRDX/DLDX-RL)/(DRDA-DLDA*DRDX/DLDX)) 152 DX = -(BL + DLDA*DA)/DLDX 153 CMIXLA(N+1) = CMIXLA(N-2)+DA 154 CMIXLA(N+2) = CMIXLA(N+1) 155 CMIXLA(N+3) = CMIXLA(N+1) 156 IF(LCALSOLZX)THEN 157 RESCAL(3,N+1) = RESCAL(3,N-2)*CALSOLZX/SOLZX 158 DZ = RESCAL(3,N+1) - RESCAL(3,N-2) 159 ZENV0A(N+1) = RESCAL(3,N+1) 160 ZENV0A(N+2) = RESCAL(3,N+1) 161 ZENV0A(N+3) = RESCAL(3,N+1) 162 ENDIF 163 RESCAL(2,N+1) = RESCAL(2,N-2)+DX 164 XENV0A(N+1) = RESCAL(2,N+1) 165 XENV0A(N+2) = RESCAL(2,N+1) 166 XENV0A(N+3) = RESCAL(2,N+1) 167 WRITE(IOWR,*) "New BL, Old BL, Delta BL: ", 168 * * BL, BLP, BL-BLP 169 WRITE(IOWR,*) "New RL, Old RL, Delta RL: ", 170 * * RL, RLP, RL-RLP 171 WRITE(IOWR,*) "New X, Old X, DX: ", 172 * * RESCAL(2,N+1), RESCAL(2,N-2), DX 173 WRITE(IOWR,*) "New A, Old A, DA: ", 174 * * CMIXLA(N+1), CMIXLA(N-2), DA 175 IF(LCALSOLZX)THEN 176 WRITE(IOWR,*) "New Z, Old Z, DZ: ", 177 * * RESCAL(3,N+1), RESCAL(3,N-2), DZ 178 ENDIF 179 BLP = BL 180 RLP = RL 181 ZLP = ZL 182 GOTO 9999 183C ENDIF 184 9999 9999 CONTINUE 185 RETURN 186 END

1"""Python version of chkcal.f (YREC 7be4133): the solar calibration's 2check and Newton step.""" 3 4import math 5 6from pyrec import iounits_mod 7from pyrec import nmlpar_mod 8from pyrec import runloop_mod 9 10# CHKCAL has SAVE: BLP and RLP are stored on one call and differenced on 11# the next; the others keep their values too. 12blp = 0.0 13rlp = 0.0 14dx = 0.0 15da = 0.0 16dz = 0.0 17zl = 0.0 18zlp = 0.0 19 20 21def chkcal(bl: float, rl: float, n: int, solzx: float, dldx: float, drdx: float, dlda: float, drda: float, 22 lsol: bool) -> tuple[float, float, float, float, bool]: 23 """Check a solar model against L_sun, R_sun (and the solar Z/X); if it 24 is not calibrated, set new X, alpha (and Z) for kind cards n+1 to n+3 25 by a Newton step with fixed derivatives. 26 27 Translates: chkcal.f#CHKCAL (7be4133) 28 29 Arguments: bl, log L/L_sun; rl, log R/R_sun; n, the current kind card; 30 solzx, the model's surface Z/X; the derivatives dldx, drdx, dlda, drda 31 (set here); lsol, whether the calibration has converged. 32 33 Returns (dldx, drdx, dlda, drda, lsol). The new X, alpha and Z go to 34 runloop_mod.rescal and nmlpar_mod.xenv0a, cmixla, zenv0a. 35 36 """ 37 38 global blp, rlp, dx, da, dz, zl, zlp 39 nml = nmlpar_mod 40 rescal = runloop_mod.rescal 41 iowr = iounits_mod.iowr 42 43 # Converged or not (Fortran lines 59–72). 44 # Question (see the page): Z/X not converged. 45 if abs(bl) < nml.toll and abs(rl) < nml.tolr: 46 if nml.lcalsolzx and nml.calsolzx > 0.0: 47 zl = math.log10(solzx) - math.log10(nml.calsolzx) 48 if abs(zl) < nml.tolz: 49 lsol = True 50 return dldx, drdx, dlda, drda, lsol 51 else: 52 lsol = True 53 return dldx, drdx, dlda, drda, lsol 54 55 # Fixed, empirical derivatives (LLP 6/28/09) (Fortran lines 74–84). 56 dldx = -3.78 # empirical result: -3.783 RMS error .070 57 drdx = -0.89 # empirical result: -0.890 RMS error .048 58 dlda = 0.0139 # empirical result: +0.0139 RMS error .0022 59 drda = -0.050 # empirical result: -0.0504 RMS error .0059 60 61 # GOTO 1234 bypasses the older method, which measured the derivatives 62 # with extra runs at kind cards 3, 6 and 9; unused (Fortran lines 85–146). 63 if False: 64 if n == 3: 65 rescal[2, n + 1] = rescal[2, 1] + dx 66 nml.xenv0a[n + 1] = rescal[2, n + 1] 67 nml.xenv0a[n + 2] = nml.xenv0a[n + 1] 68 nml.xenv0a[n + 3] = nml.xenv0a[n + 1] 69 blp = bl 70 rlp = rl 71 return dldx, drdx, dlda, drda, lsol 72 elif n == 6: 73 nml.cmixla[n + 1] = nml.cmixla[1] + da 74 nml.cmixla[n + 2] = nml.cmixla[n + 1] 75 nml.cmixla[n + 3] = nml.cmixla[n + 1] 76 rescal[2, n + 1] = rescal[2, 1] 77 nml.xenv0a[n + 1] = nml.xenv0a[1] 78 nml.xenv0a[n + 2] = nml.xenv0a[n + 1] 79 nml.xenv0a[n + 3] = nml.xenv0a[n + 1] 80 blp = bl 81 rlp = rl 82 return dldx, drdx, dlda, drda, lsol 83 elif n == 9: 84 dlda = (bl - blp) / da 85 drda = (rl - rlp) / da 86 iounits_mod.write(iowr, f" DA,DLDA,BL,BLP,DRDA,RL,RLP: {da} " 87 f"{dlda} {bl} {blp} {drda} {rl} {rlp}") 88 da = ((bl * drdx / dldx - rl) / (drda - dlda * drdx / dldx)) 89 dx = -(bl + dlda * da) / dldx 90 nml.cmixla[n + 1] = nml.cmixla[1] + da 91 nml.cmixla[n + 2] = nml.cmixla[n + 1] 92 nml.cmixla[n + 3] = nml.cmixla[n + 1] 93 rescal[2, n + 1] = rescal[2, 1] + dx 94 nml.xenv0a[n + 1] = rescal[2, n + 1] 95 nml.xenv0a[n + 2] = rescal[2, n + 1] 96 nml.xenv0a[n + 3] = rescal[2, n + 1] 97 blp = bl 98 rlp = rl 99 iounits_mod.write(iowr, f" New X, Old X, Calc DX: " 100 f"{rescal[2, n + 1]} {rescal[2, 1]} {dx}") 101 iounits_mod.write(iowr, f" New A, Old A, Calc DA: " 102 f"{nml.cmixla[n + 1]} {nml.cmixla[1]} {da}") 103 return dldx, drdx, dlda, drda, lsol 104 105 # label 1234: the Newton step, from the values of card n-2 (Fortran lines 148–155). 106 da = ((bl * drdx / dldx - rl) / (drda - dlda * drdx / dldx)) 107 dx = -(bl + dlda * da) / dldx 108 nml.cmixla[n + 1] = nml.cmixla[n - 2] + da 109 nml.cmixla[n + 2] = nml.cmixla[n + 1] 110 nml.cmixla[n + 3] = nml.cmixla[n + 1] 111 112 # The new Z, with the Z/X option (Fortran lines 156–162). 113 # Suspected bug (see the page): tests only lcalsolzx, so a target of 0 sets Z = 0. 114 if nml.lcalsolzx: 115 rescal[3, n + 1] = rescal[3, n - 2] * nml.calsolzx / solzx 116 dz = rescal[3, n + 1] - rescal[3, n - 2] 117 nml.zenv0a[n + 1] = rescal[3, n + 1] 118 nml.zenv0a[n + 2] = rescal[3, n + 1] 119 nml.zenv0a[n + 3] = rescal[3, n + 1] 120 121 # The new X (Fortran lines 163–166). 122 rescal[2, n + 1] = rescal[2, n - 2] + dx 123 nml.xenv0a[n + 1] = rescal[2, n + 1] 124 nml.xenv0a[n + 2] = rescal[2, n + 1] 125 nml.xenv0a[n + 3] = rescal[2, n + 1] 126 127 # The report to IOWR (Fortran lines 167–178). 128 # Suspected bug (see the page): blp and rlp are not yet set on the first call. 129 iounits_mod.write(iowr, f" New BL, Old BL, Delta BL: {bl} {blp} " 130 f"{bl - blp}") 131 iounits_mod.write(iowr, f" New RL, Old RL, Delta RL: {rl} {rlp} " 132 f"{rl - rlp}") 133 iounits_mod.write(iowr, f" New X, Old X, DX: {rescal[2, n + 1]} " 134 f"{rescal[2, n - 2]} {dx}") 135 iounits_mod.write(iowr, f" New A, Old A, DA: {nml.cmixla[n + 1]} " 136 f"{nml.cmixla[n - 2]} {da}") 137 if nml.lcalsolzx: 138 iounits_mod.write(iowr, f" New Z, Old Z, DZ: {rescal[3, n + 1]} " 139 f"{rescal[3, n - 2]} {dz}") 140 141 # Kept for the next call (Fortran lines 179–181). 142 blp = bl 143 rlp = rl 144 zlp = zl 145 return dldx, drdx, dlda, drda, lsol