1C
2C
3C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
4C SURFBC
5C$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
6 SUBROUTINE SURFBC(TRIT,TRIL,CFENV,PS,TS,RS,FTRI,ISTORE,ESTORE,
7C * LNEW,LRESET,LSBC,KSAHA,KENV,KATM,HSTOT,BL, ! KC 2025-05-31
8 * * LNEW,LRESET,KSAHA,KENV,KATM,HSTOT,BL,
9 * * TEFFL,X,Z,FPL,FTL,LFLP,PLIM,LC,M,CMIXL)
10C 10/15/26 MHP via Claude
11C code changed to take JSON from DIMENSIONS_MOD instead of declaring
12C PARAMETER(JSON=5000) here.
13 USE DIMENSIONS_MOD
14C 10/14/26 MHP via Claude
15C code changed to convert COMMON/COMP/ into COMP_MOD.
16 USE COMP_MOD
17C 10/14/26 MHP via Claude
18C code changed to convert COMMON/ALATM05/ into ALATM05_MOD.
19 USE ALATM05_MOD
20C 10/14/26 MHP via Claude
21C code changed to convert COMMON/NMLPAR/ into NMLPAR_MOD.
22 USE NMLPAR_MOD
23C 10/14/26 MHP via Claude
24C code changed to convert COMMON/CONSTANTS/ into CONSTANTS_MOD.
25 USE CONSTANTS_MOD
26C 10/14/26 MHP via Claude
27C code changed to convert COMMON/IOUNITS/ into IOUNITS_MOD.
28 USE IOUNITS_MOD
29C 10/20/26 MHP via Claude
30C code changed to replace implicit typing with IMPLICIT NONE and
31C explicit declarations.
32 IMPLICIT NONE
33 REAL*8 ATEFFL,B,BL,CFENV,CMIXL,ESTORE,FPL,FTL,FTRI,GL,HSTOT,PLIM,
34 * * PS,RL,RS,TEFFL,TEMP,TEMP1,TEMP2,TRIERR,TRIL,TRIT,TS,VTEFFL,
35 * * X,Z,DEN
36 INTEGER I,I1,I2,I3,ISTORE,IXX,J,KATM,KENV,KSAHA,M,NUMENV,NREFL
37 LOGICAL*4 LC,LFLP,LNEW,LPENV,LPULPT,LRENV,LRESET,LSBC0,LTRI
38
39C INPUTS LNEW = .T. START UP WITH 3 NEW ENVELOPES ABOUT(TEFFL,BL)
40C INPUTS LRESET = .T. REDO ALL 3 ENVELOPES AND RETRIANGULATE IF NEED
41C BOTH LNEW AND LRESET ARE RESET TO .FALSE.
42
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/
49 DIMENSION LTRI(3),TRIT(3),TRIL(3),CFENV(9),PS(3),TS(3),RS(3),
50 * * ESTORE(4),LC(JSON)
51C 10/11/26 MHP via Claude
52C code changed to remove LSCRIB from COMMON/CCOUT/: it was declared
53C positionally in all 17 files sharing this block but never read or
54C written in executable code anywhere, and no namelist file sets it.
55C 10/12/26 MHP via Claude
56C code changed to extend COMMON/NMLPAR/ to 234 members, adding the
57C namelist variables that the physics routines read from 27 COMMON
58C blocks which also held derived or working members. Those blocks keep
59C only their non-namelist members; five of them -- NEWMX, VMULT,
60C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone.
61C Derivation-only namelist variables from the same blocks went to
62C COMMON/NMLBYPASS/ instead, not here.
63C
64C Declared identically in all 121 files that hold any of these members.
65C Generated from one member list, so every copy is byte-identical; the
66C linker's 'size of symbol changed' warning is the check that they are.
67C 10/13/26 MHP via Claude: removed the unused COMMON/CCOUT/ declaration
68C (no member referenced here; residue of the obsolete-option removals).
69C 10/13/26 MHP via Claude: removed the unused COMMON/CENV/ declaration
70C (no member referenced here; residue of the obsolete-option removals).
71C 9/25/26 MHP via Claude
72C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/,
73C and COMMON/CONST3/ into a single COMMON/CONSTANTS/
74
75C 10/13/26 MHP via Claude: removed the unused COMMON/ATMOS/ declaration
76C (no member referenced here; residue of the obsolete-option removals).
77 DATA NUMENV/0/
78
79C G Somers END
80
81 SAVE
82C MHP 9/01
83 IF(KTTAU .EQ. 3) THEN
84 IF(TEFFL.GE. 3.95D0) THEN
85 write(*,*)
86 WRITE(*,5) TEFFL
87 5 5 FORMAT('LOG TEFF OF ',F7.3,' ABOVE 3.95 - SWITCH'
88 * * ,'TO GRAY ATMOSPHERE BOUNDARY CONDITION')
89 write(*,*)
90 KTTAU = 0
91 LNEW = .TRUE.
92C MHP 06/13 Remember that flag is switched
93 LTTAU = .TRUE.
94 ENDIF
95 ENDIF
96 IF(KTTAU. EQ. 4) THEN
97 IF(TEFFL.GE. AL_TEFFLmax) THEN
98 write(*,*)
99 WRITE(*,7) TEFFL,AL_TEFFLmax
100 7 7 FORMAT('LOG TEFF OF ',F7.3,' ABOVE Allard Table max ',F7.3
101 * * ,' - SWITCH TO GRAY ATMOSPHERE BOUNDARY CONDITION')
102 write(*,*)
103 KTTAU = 0
104 LNEW = .TRUE.
105C MHP 06/13 Remember that flag is switched
106 LTTAU = .TRUE.
107 ENDIF
108 ENDIF
109 IF(LTTAU)THEN
110 IF(KTTAU0.EQ.3.AND.TEFFL.LT.3.95D0)THEN
111 KTTAU = KTTAU0
112 LTTAU = .FALSE.
113C 10/30/26 MHP via Claude
114C code changed to start a new triangle on the switch back, as the switch
115C to gray does; the old vertices were integrated with the gray T(tau).
116 LNEW = .TRUE.
117 WRITE(*,9) TEFFL
118 9 9 FORMAT('LOG TEFF OF ',F7.3,' BELOW 3.95 - SWITCH'
119 * * ,' BACK TO KURUCZ ATMOSPHERE BOUNDARY CONDITION')
120 ELSE IF(KTTAU0.EQ.4.AND.TEFFL.LT.AL_TEFFLmax)THEN
121 KTTAU = KTTAU0
122 LTTAU = .FALSE.
123C 10/30/26 MHP via Claude: new triangle here too (see above).
124 LNEW = .TRUE.
125 WRITE(*,11) TEFFL,AL_TEFFLmax
126 11 11 FORMAT('LOG TEFF OF ',F7.3,' below Allard Table max ',F7.3
127 * * ,' - SWITCH BACK TO ALLARD ATMOSPHERE BOUNDARY CONDITION')
128 ENDIF
129 ENDIF
130 LFLP = .FALSE.
131C IF LNEW0, REDO ALL 3 ENVELOPES EVERY MODEL
132 IF(LNEW.OR.LRESET) THEN
133C REDO ALL THREE ENVELOPES
134 LTRI(1) = .FALSE.
135 LTRI(2) = .FALSE.
136 LTRI(3) = .FALSE.
137 ISTORE = 0
138 TRIERR = 0.0D0
139 LRESET = .FALSE.
140 ELSE
141 LTRI(1) = .TRUE.
142 LTRI(2) = .TRUE.
143 LTRI(3) = .TRUE.
144 TRIERR = 0.0D0
145 ENDIF
146 IF(LNEW) THEN
147C STARTING PROCEDURE
148 FTRI = +1.0D0
149 TRIT(3) = TEFFL
150 TRIT(1) = TRIT(3) - 0.5D0*TRIDT
151 TRIT(2) = TRIT(1) + TRIDT
152 TRIL(3) = BL + 0.5D0*TRIDL
153 TRIL(1) = TRIL(3) - TRIDL
154 TRIL(2) = TRIL(1)
155 LNEW = .FALSE.
156 ELSE
157C CHECK TRIANGULATION OF POINT (TEFFL,BL)
158C 10/30/26 MHP via Claude: limit the reflection walk. It never ended for
159C a degenerate triangle; the suite's decks need at most a few reflections.
160 NREFL = 0
161 10 10 CONTINUE
162 DO 20 I1 = 1,3
163 I2 = MOD(I1,3) + 1
164 I3 = MOD(I2,3) + 1
165 TEMP = FTRI*((TRIL(I2)-TRIL(I3))*(TEFFL-TRIT(I2)) +
166 * * (TRIT(I3)-TRIT(I2))*(BL -TRIL(I2)) )
167 IF(TEMP.LT.-TRIERR) THEN
168 FTRI = -FTRI
169 TRIT(I1) = TRIT(I2) + TRIT(I3) - TRIT(I1)
170 TRIL(I1) = TRIL(I2) + TRIL(I3) - TRIL(I1)
171 LTRI(I1) = .FALSE.
172 TRIERR = 0.0D0
173 NREFL = NREFL + 1
174 IF(NREFL.GT.1000)THEN
175 WRITE(*,60)TEFFL,BL
176 WRITE(ISHORT,60)TEFFL,BL
177 60 60 FORMAT(' SURFBC: ENVELOPE TRIANGLE WALK FAILED AFTER',
178 * * ' 1000 REFLECTIONS AT LOG TEFF',F9.5,' LOG L',
179 * * F9.5/' RUN STOPPED')
180 STOP 1
181 ENDIF
182 GOTO 10
183 ENDIF
184 20 20 CONTINUE
185 ENDIF
186C COMPUTE NEW ENVELOPES IF NECESSARY
187 LRENV = .FALSE.
188 DO 40 I = 1,3
189 IF(.NOT.LTRI(I)) THEN
190C NEW ENVELOPE NEEDED
191 LFLP = .TRUE.
192C 10/29/26 MHP via Claude
193C code changed to integrate the vertex at a local VTEFFL. Writing the
194C vertex Teff into TEFFL changed the model's own log Teff, which CRRECT
195C then passed on (TPGRAD's depth-dependent spot contrast) until it
196C recomputed TEFFL after the corrections.
197 VTEFFL = TRIT(I)
198 B = DEXP(CLN*TRIL(I))
199 RL = 0.5D0*(TRIL(I) + CLSUNL - 4.0D0*VTEFFL - C4PIL - CSIGL)
200 GL = CGL + STOTAL - RL - RL
201 NUMENV = NUMENV + 1
202C G Somers 11/14, LPENV FUNCTIONALITY NOW INCLUDED IN LSTENV.
203 LPENV = .FALSE.
204C LPENV = MOD(NUMENV,NPENV).EQ.0
205C IF(.NOT.LSTORE) LPENV = .FALSE.
206C IF(LPENV) THEN
207C WRITE(ISHORT,30) I,TRIT(I),TRIL(I)
208C 30 FORMAT(/,' NEW ENVELOPE NO.',I2,' LOG(TE) =',F9.5,
209C * ' LOG(L) =',F8.5)
210C ENDIF
211C G Somers END
212 IXX=I
213 LSBC0 = .TRUE.
214C DBG PULSE: DO NOT DO PULSE OUTPUT
215 LPULPT = .FALSE.
216C G Somers 10/14, FOR SPOTTED RUNS, FIND THE
217C PRESSURE AT THE AMBIENT TEMPERATURE ATEFFL
218 IF(LC(M).AND.SPOTF.NE.0.0.AND.SPOTX.NE.1.0)THEN
219 ATEFFL = VTEFFL - 0.25*LOG10(SPOTF * SPOTX**4.0 + 1.0 - SPOTF)
220 ELSE
221 ATEFFL = VTEFFL
222 ENDIF
223 CALL ENVINT(B,FPL,FTL,GL,HSTOT,IXX,LPENV,LSBC0,PLIM,RL,
224 * * ATEFFL,X,Z,ESTORE,ISTORE,KATM,KENV,KSAHA,PS,
225 * * RS,TS,LPULPT,CMIXL)
226C G Somers END
227 LRENV = .TRUE.
228 ENDIF
229 40 40 CONTINUE
230 IF(LRENV) THEN
231C RECOMPUTE COEFFICIENTS
232C 10/30/26 MHP via Claude
233C code changed to refuse a degenerate plane fit: two fitting points at
234C the same log T, or three collinear in (log P, log T), gave Inf or NaN
235C coefficients silently. The suite's decks stay above 6e-3 of the
236C differences divided, against the 1e-10 refused here.
237 IF(ABS(TS(2)-TS(1)).LE.1.0D-10*ABS(TS(3)-TS(1)))THEN
238 WRITE(*,61)(TS(I),PS(I),I=1,3)
239 WRITE(ISHORT,61)(TS(I),PS(I),I=1,3)
240 STOP 1
241 ENDIF
242 TEMP = 1.0D0/(TS(2)-TS(1))
243 TEMP1 = (TS(3)-TS(1))*TEMP
244 DEN = PS(3)-PS(1)-TEMP1*(PS(2)-PS(1))
245 IF(ABS(DEN).LE.1.0D-10*(ABS(PS(3)-PS(1))+
246 * * ABS(TEMP1*(PS(2)-PS(1)))))THEN
247 WRITE(*,61)(TS(I),PS(I),I=1,3)
248 WRITE(ISHORT,61)(TS(I),PS(I),I=1,3)
249 61 61 FORMAT(' SURFBC: DEGENERATE ENVELOPE TRIANGLE FIT; FITTING',
250 * * ' POINTS (LOG T, LOG P):'/3(2F12.7/)' RUN STOPPED')
251 STOP 1
252 ENDIF
253 TEMP2 = 1.0D0/DEN
254 CFENV(1) = ( RS(3)- RS(1)-TEMP1*( RS(2)- RS(1)))*TEMP2
255 CFENV(4) = (TRIL(3)-TRIL(1)-TEMP1*(TRIL(2)-TRIL(1)))*TEMP2
256 CFENV(7) = (TRIT(3)-TRIT(1)-TEMP1*(TRIT(2)-TRIT(1)))*TEMP2
257 CFENV(2) = ( RS(2)- RS(1)-CFENV(1)*(PS(2)-PS(1)))*TEMP
258 CFENV(5) = (TRIL(2)-TRIL(1)-CFENV(4)*(PS(2)-PS(1)))*TEMP
259 CFENV(8) = (TRIT(2)-TRIT(1)-CFENV(7)*(PS(2)-PS(1)))*TEMP
260 CFENV(3) = RS(1) - CFENV(1)*PS(1) - CFENV(2)*TS(1)
261 CFENV(6) = TRIL(1) - CFENV(4)*PS(1) - CFENV(5)*TS(1)
262 CFENV(9) = TRIT(1) - CFENV(7)*PS(1) - CFENV(8)*TS(1)
263 WRITE(ISHORT,50)(I,LTRI(I),TRIT(I),TRIL(I),PS(I),TS(I),RS(I),
264 * * (CFENV(I+I+I-3+J),J=1,3), I=1,3)
265 50 50 FORMAT(' ENVELOPE TRIANGLE LOG(TE) LOG(L) LOG(P) LOG(T)',
266 * * 3X,'LOG(R) COEFFICIENTS'/I14,L1,2X,5F9.5,8X,'LOG(R)',
267 * * 3F10.5/ I14,L1,2X,5F9.5,5X,'LOG(L/L0)',3F10.5/
268 * * I14,L1,2X,5F9.5,7X,'LOG(TE)',3F10.5)
269C 60 FORMAT(' COUNTS KATM',I5,' KENV',I5,' KSAHA',I5)
270 ELSE
271 WRITE(ISHORT,70)
272 70 70 FORMAT(' ENVELOPE TRIANGLE NOT CHANGED')
273 ENDIF
274 RETURN
275 END1"""Python version of surfbc.f (YREC 7be4133): the outer boundary condition
2from the envelope triangle."""
3
4import math
5
6import numpy as np
7
8from pyrec import alatm05_mod
9from pyrec import comp_mod
10from pyrec import constants_mod
11from pyrec import iounits_mod
12from pyrec import nmlpar_mod
13from pyrec.envint import envint
14
15# DATA NUMENV/0/ with SAVE: a count of the envelopes computed, kept between
16# calls (not the envelope's point count in envstruct_mod) (Fortran line 77).
17numenv = 0
18
19
20def surfbc(trit: np.ndarray, tril: np.ndarray, cfenv: np.ndarray, ps: np.ndarray, ts: np.ndarray, rs: np.ndarray,
21 ftri: float, istore: int, estore: np.ndarray, lnew: bool, lreset: bool, ksaha: int, kenv: int, katm: int,
22 hstot: float, bl: float, teffl: float, x: float, z: float, fpl: float, ftl: float, lflp: bool, plim: float,
23 lc: np.ndarray, m: int, cmixl: float) -> tuple[float, int, bool, bool, int, int, int, bool]:
24 """Keep three envelope integrations around (log Teff, log L), flip the
25 triangle when the model leaves it, and fit the plane cfenv that gives
26 log R, log L and log Teff from log P and log T at the fitting point;
27 switch between atmosphere tables and the gray atmosphere as Teff
28 changes.
29
30 Translates: surfbc.f#SURFBC (7be4133)
31
32 Arguments: the triangle trit, tril (log Teff, log L of its vertices,
33 index 1-3) and its envelopes' fitting-point values ps, ts, rs (log P,
34 log T, log R), cfenv (the plane's 9 coefficients), ftri (orientation),
35 istore, estore (passed to envint), lnew (start a new triangle), lreset
36 (recompute all three envelopes), the counters ksaha, kenv, katm, hstot
37 (log M), bl (log L/Lsun), teffl (log Teff), the surface x, z, the
38 rotation factors fpl, ftl, plim, lc (convective flags; lc[m] for the
39 spots), m, and cmixl (the mixing length, passed to envint). Arrays are
40 changed in place; lflp is output only. Stops (SystemExit) when the
41 triangle's reflection walk exceeds 1000 reflections or the plane fit is
42 degenerate.
43
44 Returns (ftri, istore, lnew, lreset, ksaha, kenv, katm, lflp): lflp is
45 True if an envelope was recomputed. teffl is no longer changed (7be4133:
46 each vertex is integrated at a local vteffl).
47
48 """
49
50 global numenv
51 c = constants_mod
52 nml = nmlpar_mod
53 ltri = [None, False, False, False]
54
55 # Above the atmosphere table's range, switch to the gray atmosphere and
56 # start a new triangle (MHP 9/01) (Fortran lines 82–108).
57 # Question (see the page): the limit 3.95 is hard-coded; KTTAU = 5 has
58 # no switch.
59 if nml.kttau == 3:
60 if teffl >= 3.95:
61 print()
62 print(f"LOG TEFF OF {teffl:7.3f} ABOVE 3.95 - SWITCH"
63 f"TO GRAY ATMOSPHERE BOUNDARY CONDITION") # format 5
64 print()
65 nml.kttau = 0
66 lnew = True
67 nml.lttau = True # remember that the flag is switched (MHP 06/13)
68 if nml.kttau == 4:
69 if teffl >= alatm05_mod.al_tefflmax:
70 print()
71 print(f"LOG TEFF OF {teffl:7.3f} ABOVE Allard Table max "
72 f"{alatm05_mod.al_tefflmax:7.3f}"
73 f" - SWITCH TO GRAY ATMOSPHERE BOUNDARY CONDITION")
74 print()
75 nml.kttau = 0
76 lnew = True
77 nml.lttau = True
78
79 # Switch back below the limit, with a new triangle, since the old
80 # vertices were integrated with the gray T(tau) (MHP 10/30/26)
81 # (Fortran lines 109–129).
82 if nml.lttau:
83 if nml.kttau0 == 3 and teffl < 3.95:
84 nml.kttau = nml.kttau0
85 nml.lttau = False
86 lnew = True
87 print(f"LOG TEFF OF {teffl:7.3f} BELOW 3.95 - SWITCH"
88 f" BACK TO KURUCZ ATMOSPHERE BOUNDARY CONDITION")
89 elif nml.kttau0 == 4 and teffl < alatm05_mod.al_tefflmax:
90 nml.kttau = nml.kttau0
91 nml.lttau = False
92 lnew = True
93 print(f"LOG TEFF OF {teffl:7.3f} below Allard Table max "
94 f"{alatm05_mod.al_tefflmax:7.3f}"
95 f" - SWITCH BACK TO ALLARD ATMOSPHERE BOUNDARY CONDITION")
96
97 # Which envelopes to redo (Fortran lines 130–145). Question (see the
98 # page): trierr is always zero.
99 lflp = False
100 if lnew or lreset:
101 # Redo all three envelopes.
102 ltri[1] = ltri[2] = ltri[3] = False
103 istore = 0
104 trierr = 0.0
105 lreset = False
106 else:
107 ltri[1] = ltri[2] = ltri[3] = True
108 trierr = 0.0
109
110 if lnew:
111 # A new triangle about (teffl, bl) (Fortran lines 146–155).
112 ftri = 1.0
113 trit[3] = teffl
114 trit[1] = trit[3] - 0.5 * nml.tridt
115 trit[2] = trit[1] + nml.tridt
116 tril[3] = bl + 0.5 * nml.tridl
117 tril[1] = tril[3] - nml.tridl
118 tril[2] = tril[1]
119 lnew = False
120 else:
121 # Check that (teffl, bl) is inside the triangle; flip a vertex
122 # across the opposite side until it is, at most 1000 times
123 # (MHP 10/30/26) (Fortran lines 156–185).
124 nrefl = 0
125 flipped = True
126 while flipped: # label 10
127 flipped = False
128 for i1 in range(1, 4): # loop 20
129 i2 = i1 % 3 + 1
130 i3 = i2 % 3 + 1
131 temp = ftri * ((tril[i2] - tril[i3]) * (teffl - trit[i2])
132 + (trit[i3] - trit[i2]) * (bl - tril[i2]))
133 if temp < -trierr:
134 ftri = -ftri
135 trit[i1] = trit[i2] + trit[i3] - trit[i1]
136 tril[i1] = tril[i2] + tril[i3] - tril[i1]
137 ltri[i1] = False
138 trierr = 0.0
139 nrefl = nrefl + 1
140 if nrefl > 1000:
141 text = (" SURFBC: ENVELOPE TRIANGLE WALK FAILED AFTER"
142 f" 1000 REFLECTIONS AT LOG TEFF{teffl:9.5f}"
143 f" LOG L{bl:9.5f}\n RUN STOPPED") # format 60
144 print(text)
145 iounits_mod.write(iounits_mod.ishort, text)
146 raise SystemExit(1)
147 flipped = True
148 break # GOTO 10
149
150 # Compute new envelopes where needed (Fortran lines 186–229).
151 lrenv = False
152 for i in range(1, 4): # loop 40
153 if not ltri[i]:
154 lflp = True
155 # The vertex's log Teff, kept local (MHP 10/29/26).
156 vteffl = trit[i]
157 b = math.exp(c.cln * tril[i])
158 rl = 0.5 * (tril[i] + c.clsunl - 4.0 * vteffl - c.c4pil
159 - c.csigl)
160 gl = c.cgl + comp_mod.stotal - rl - rl
161 numenv = numenv + 1
162 lpenv = False # part of LSTENV (G Somers 11/14)
163 ixx = i
164 lsbc0 = True
165 lpulpt = False # no pulsation output
166 # For spotted stars, the ambient temperature (G Somers 10/14).
167 if lc[m] and nml.spotf != 0.0 and nml.spotx != 1.0:
168 ateffl = vteffl - 0.25 * math.log10(
169 nml.spotf * nml.spotx**4.0 + 1.0 - nml.spotf)
170 else:
171 ateffl = vteffl
172 rl, istore, katm, kenv, ksaha = envint(
173 b, fpl, ftl, gl, hstot, ixx, lpenv, lsbc0, plim, rl, ateffl,
174 x, z, estore, istore, katm, kenv, ksaha, ps, rs, ts, lpulpt,
175 cmixl)
176 lrenv = True
177
178 if lrenv:
179 # The plane through the three envelopes, for log R, log L, log Teff;
180 # a degenerate fit stops the run (MHP 10/30/26) (Fortran lines
181 # 230–262).
182 if abs(ts[2] - ts[1]) <= 1.0e-10 * abs(ts[3] - ts[1]):
183 _degenerate_fit(ps, ts)
184 temp = 1.0 / (ts[2] - ts[1])
185 temp1 = (ts[3] - ts[1]) * temp
186 den = ps[3] - ps[1] - temp1 * (ps[2] - ps[1])
187 if abs(den) <= 1.0e-10 * (abs(ps[3] - ps[1])
188 + abs(temp1 * (ps[2] - ps[1]))):
189 _degenerate_fit(ps, ts)
190 temp2 = 1.0 / den
191 cfenv[1] = (rs[3] - rs[1] - temp1 * (rs[2] - rs[1])) * temp2
192 cfenv[4] = (tril[3] - tril[1] - temp1 * (tril[2] - tril[1])) * temp2
193 cfenv[7] = (trit[3] - trit[1] - temp1 * (trit[2] - trit[1])) * temp2
194 cfenv[2] = (rs[2] - rs[1] - cfenv[1] * (ps[2] - ps[1])) * temp
195 cfenv[5] = (tril[2] - tril[1] - cfenv[4] * (ps[2] - ps[1])) * temp
196 cfenv[8] = (trit[2] - trit[1] - cfenv[7] * (ps[2] - ps[1])) * temp
197 cfenv[3] = rs[1] - cfenv[1] * ps[1] - cfenv[2] * ts[1]
198 cfenv[6] = tril[1] - cfenv[4] * ps[1] - cfenv[5] * ts[1]
199 cfenv[9] = trit[1] - cfenv[7] * ps[1] - cfenv[8] * ts[1]
200
201 # Format 50 (Fortran lines 263–268).
202 labels = ["", " LOG(R)", " LOG(L/L0)", " LOG(TE)"]
203 text = (" ENVELOPE TRIANGLE LOG(TE) LOG(L) LOG(P) LOG(T)"
204 " LOG(R) COEFFICIENTS")
205 for i in range(1, 4):
206 row = "".join(f"{v:9.5f}"
207 for v in (trit[i], tril[i], ps[i], ts[i], rs[i]))
208 coef = "".join(f"{cfenv[3 * i - 3 + j]:10.5f}"
209 for j in range(1, 4))
210 text += (f"\n{i:14d}{'T' if ltri[i] else 'F'} {row}"
211 f"{labels[i]}{coef}")
212 iounits_mod.write(iounits_mod.ishort, text)
213 else:
214 # Fortran lines 270–273.
215 iounits_mod.write(iounits_mod.ishort,
216 " ENVELOPE TRIANGLE NOT CHANGED") # format 70
217
218 return ftri, istore, lnew, lreset, ksaha, kenv, katm, lflp
219
220
221def _degenerate_fit(ps, ts):
222 """Write format 61 to the screen and ISHORT, then STOP 1."""
223
224 # Fortran lines 237–252.
225 text = (" SURFBC: DEGENERATE ENVELOPE TRIANGLE FIT; FITTING POINTS"
226 " (LOG T, LOG P):\n"
227 + "".join(f"{ts[i]:12.7f}{ps[i]:12.7f}\n" for i in range(1, 4))
228 + " RUN STOPPED")
229 print(text)
230 iounits_mod.write(iounits_mod.ishort, text)
231 raise SystemExit(1)