shape.f

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

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

Variables on this page (8)

acfpft, c4pi3l, cc13, cgl, cln, itfp1, itfp2, json

Purpose

Finds how much each shell of a rotating star is flattened. Its comment: "SHAPE FINDS THE DISTORTION OF A GIVEN SHELL FROM SPHERICAL SYMMETRY; IT CALCULATES THE CHARACTERISTIC RADIUS R0 AND ASSOCIATED PARAMETER ETA2. IN POLAR CO-ORDINATES THE RADIUS OF A SHELL IS R = R0(1 - AP2) … FOR MORE INFORMATION SEE ENDAL AND SOFIA."

How it fits

Who calls it: FPFT (the rotational correction factors $f_P$, $f_T$ of the structure equations) and SOLID (solid-body rotation).
Inputs: $\log\rho$, $\log r$, $\log m$ (HD, HR, HS), the angular velocities OMEGA, the range JSTART–JEND; from NMLPAR_MOD the tolerance ACFPFT and the iteration limits ITFP1, ITFP2; the constants from CONSTANTS_MOD. When JSTART > 1, the integration starts from the caller's ETA2(JSTART-1) and, in the second pass, R0(JSTART-1), so both arrays are inputs as well as outputs.
Outputs: ETA2 (Radau's $\eta_2$) and R0 (the characteristic radius) for each shell, which the callers use to compute the moments of inertia and the correction factors.

The method

With the mean density inside the shell $\bar\rho = m/(\tfrac43\pi r^3)$ and $r$ the radius of the sphere of equal volume (RPHI):

Radau's equation (as the comment states it), $$r\,\frac{d\eta_2}{dr} + 6\,\frac{\rho}{\bar\rho}\,(\eta_2 + 1) + \eta_2(\eta_2 - 1) = 6,$$ integrated outward with a fourth-order Runge–Kutta step per shell, using the averages of $\rho$ and $\bar\rho$ at the midpoint (lines 91–118 ↗). At the first shell, $\eta_2 = 6(1 - \rho/\bar\rho)$ (line 60 ↗; a suspected bug, suspected bug 1). $\eta_2$ is 0 at the center (the comment's boundary condition, line 93 ↗), and it sets how strongly a shell is flattened, $A \propto 1/(2+\eta_2)$ (step 2). The code does not define $\eta_2$ beyond Radau's equation; what it stands for is question 1.

Question 1 (lines 40–45 ↗; lines 92–93 ↗). What $\eta_2$ is. The code names $\eta_2$ only through Radau's equation. Is it the logarithmic derivative of the distortion amplitude, $\eta_2 = d\ln A/d\ln r_0$ (Radau's equation being Clairaut's equation for the shells' distortion rewritten in it, so that a shell's flattening depends on how the shells inside it are distorted)? FPFT's step 2 uses the same reading. If so, it should approach 3 at the surface of a centrally condensed star; does it in YREC's models? 2. The characteristic radius $r_0$ solves, by Newton's method (at most ITFP2 steps; the test compares the relative change of $r_0^3$ with ACFPFT itself, $10^{-36}$ by default, below double precision, so in practice the loop ends only when the update rounds to zero or after all ITFP2 steps: line 71 ↗, line 134 ↗, a suspected bug, suspected bug 2), $$r^3 = r_0^3\Bigl(1 + \tfrac35 A^2 - \tfrac{2}{35}A^3\Bigr),\qquad A = \frac{5\,\Omega^2 r_0^3}{3\,G m\,(2 + \eta_2)} .$$ 3. Iteration between the two (at most ITFP1 times, lines 122–162 ↗): the first integration assumes $r_0 = r$; if $r_0$ differs from that guess by more than ACFPFT$^{1/3}$, $\eta_2$ is integrated again in $r_0$ and $r_0$ recomputed.

Question 2 (lines 137–140 ↗; lines 159–160 ↗). Tolerances compared with different quantities. In the outer loop, the first test compares the relative change of $r_0$ with ACC = ACFPFT$^{1/3}$, and the second compares the absolute change of $\eta_2$ (ERR = ETA2T - ETA2(I)) with the same ACC. Is that intended?

Question 3 (lines 99–100 ↗; lines 144–156 ↗). Mixed variables in the second pass. In the second Runge–Kutta pass, the step is in $r_0$, but the midpoint densities RHOA, RHOBA are the same averages as in the first pass, taken at the volume radii. Is this approximation intended?

Question 4 (lines 66–72 ↗; lines 124–162 ↗). Silent non-convergence. Neither the Newton loop (ITFP2 steps) nor the outer loop (ITFP1 passes) warns when it ends without meeting its tolerance; the last values are kept. Should it warn?

From YREC's documentation

Physics Options (Structural effects of rotation): ACFPFT (default 1.0D-36) is the fractional accuracy of the departures from spherical symmetry, given cubed (so $10^{-12}$ by default); ITFP1 (5) is the number of iterations between the equipotential shape and the effective radius, ITFP2 (20) that for the effective radius at a given shape parameter.

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 first shell's $\eta_2$ is seven times too large.

(line 60 ↗; line 93 ↗) The first shell takes $\eta_2 = 6(1-\rho/\bar\rho)$ (ETA2(1) = 6.D0*(1.0D0 - RHO/RHOBAR)). The series solution of Radau's equation as the comment states it (line 93 ↗) near the centre ($\rho = \rho_c(1-kr^2)$, so $\rho/\bar\rho \approx 1 - \tfrac25 kr^2$, and $\eta_2 = cr^2$) gives $7c = \tfrac{12}{5}k$, so $\eta_2 \approx \tfrac67(1-\rho/\bar\rho)$, a seventh of the code's value. What follows: The integration starts from a wrong $\eta_2$; the relative error decays outward as $r^{-7}$, so the effect on ETA2, R0, and so on FPFT's $f_P$, $f_T$ and MOMI's moments of inertia is probably small after a few shells. Only with rotation (LROT, off by default).

A possible fix:

shape.f

@@ lines 57–63 @@
          RHOBAR = DEXP(CLN*(HS(1) - C4PI3L - 3.0D0*HR(1)))
          RPHI = DEXP(CLN*HR(1))
          RHO = DEXP(CLN*HD(1))
-         ETA2(1) = 6.D0*(1.0D0 - RHO/RHOBAR)
+C 11/09/26 KC via Claude
+C code changed to start ETA2 from the series solution of Radau's
+C equation near the centre, ETA2 = 6/7*(1 - RHO/RHOBAR).
+         ETA2(1) = 6.D0/7.D0*(1.0D0 - RHO/RHOBAR)
 C ITERATE FOR R0 GIVEN RPHI AND ETA2
          GM = CG*DEXP(CLN*HS(1))
          RPHI3 = RPHI**3

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 first shell's $\eta_2$ was set in all five rotating test cases, from 26,314 to 90,623 times, with values from $-1.80$ to 0.68. With the fix, the 3 $M_\odot$ case (testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2) stays within the noise: $\log L$ moves by up to 0.0012 dex, against 0.0076 for a null change (nearly every opacity, all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$, times $1 + 10^{-10}$, on the same case), and central $Y$ by 0.018 against 0.014 (1,863 models instead of 1,853). In the 0.8 $M_\odot$ case (testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS) $\log L$ and $\log T_\mathrm{eff}$ stay within the noise, and central $X$ and $Y$ moved by up to 0.00015 at 0.395 Gyr, about four times a null change, as with TRAPZD's fix (TRAPZD, bug 1). Rechecked on Nov 10, 2026, with a build without counting lines in SHAPE and TRAPZD, and beside it a null change of the distortion itself (the first shell's $\eta_2$ times $1 + 10^{-10}$): Model by model, every change, the null one included, stays below $10^{-7}$ in central $X$ until about 0.046 Gyr, where central $X$ stops falling, and grows to some $10^{-5}$ by 0.16 Gyr, so this case amplifies whatever perturbs it there, and the size of the change depends on the build. In that build, SHAPE's fix moves central $X$ by up to $5.8\times10^{-5}$ at 0.222 Gyr, against $2.9\times10^{-5}$ for the null change at the same age ($6.8\times10^{-5}$ and $6.1\times10^{-5}$ model by model), while TRAPZD's fix again gives 0.00015 at 0.3955 Gyr: The shared number came not from one fix acting through the other but from the case's sensitivity, and SHAPE's fix is within the noise there. The 2 $M_\odot$ cases, which rotate faster, were not run with the fix. Intent: none found.

2. (harmless in the runs) The Newton test for $r_0$ compares with `ACFPFT`, not its cube root.

(line 71 ↗; line 122 ↗; line 134 ↗) The documentation describes ACFPFT (default $10^{-36}$) as the fractional accuracy of the departures from spherical symmetry "given cubed," and the outer test uses its cube root, ACC ($10^{-12}$). The Newton tests compare DABS(DELR03/R03) with ACFPFT itself, so they pass only when the update is exactly zero, or, in the calls from CHECKJ before its last iteration, which relax ACFPFT to $10^{-6}$ (checkj.f, lines 195–203 ↗), when it is below $10^{-6}$; otherwise all ITFP2 = 20 steps run. What follows: The setting is not applied as documented, and the time spent in SHAPE grows; the results are converged to machine precision, so they should change only at about $10^{-12}$ with the fix. With rotation only (LROT, off by default).

A possible fix:

shape.f

@@ lines 48–53 @@
       C2 = 2.0D0/3.5D1
       C3 = 3.0D0*C1
       C4 = 4.0D0*C2
+C 11/09/26 KC via Claude
+C code changed to compare the Newton update of R0**3 with ACC, the
+C cube root of ACFPFT, as the outer test does.
+      ACC = ACFPFT**CC13
       JBEGIN = JSTART
       IF(JSTART.EQ.1) THEN
 C  CALCULATE FIRST POINT ETA2 USING CENTRAL B.C.
@@ lines 68–74 @@
             DELR03 = (RPHI3-R03*(1.0D0 + C1*A**2 - C2*A**3))/
      *      (1.0D0 + C3*A**2 - C4*A**3)
             R03 = R03 + DELR03
-            IF(DABS(DELR03/R03).LE.ACFPFT)GOTO 20
+            IF(DABS(DELR03/R03).LE.ACC)GOTO 20
    10    CONTINUE
    20    R0(1) = R03**CC13
          IF(JEND.EQ.1) GOTO 9999
@@ lines 119–125 @@
          RPHI3 = RPHI**3
          R03 = RPHI3
          R0EST = RPHI
-         ACC = ACFPFT**CC13
 C ITERATE BETWEEN SOLUTION FOR ETA2 AND SOLUTION FOR R0 ITFP1 TIMES.
          DO 60 K = 1,ITFP1
             FACT = 5.0D0*CC13*OMEGA(I)**2/(GM*(2.0D0+ETA2(I)))
@@ lines 131–137 @@
                DELR03 = (RPHI3-R03*(1.0D0 + C1*A**2 - C2*A**3))/
      *         (1.0D0 + C3*A**2 - C4*A**3)
                R03 = R03 + DELR03
-               IF(DABS(DELR03/R03).LT.ACFPFT)GOTO 50
+               IF(DABS(DELR03/R03).LT.ACC)GOTO 50
    40       CONTINUE
    50       R0(I) = R03**CC13
             ERR = R0(I) - R0EST

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 Newton solve for $r_0$ ran about 29 to 95 million times in each rotating test case. It nearly always ended early (after 1 to 20 steps), with an update of zero, or of up to $10^{-6}$ in CHECKJ's relaxed calls; it ran all 20 steps, with updates of $1$ to $4\times10^{-16}$ that had long met the cube-root accuracy, in 47 (testsuite/Test_3Msolar_GS_rot_gray_OPALSCV_SF2_GS98_OP_CF10_RGp2diffrotI3cz2) to 15,909 (testsuite/Test_m2p0_feh+0p0_vzams150_IES1_TAMS) solves. With the fix, the 2 $M_\odot$ case stays within the noise (every track difference at most $1.6\times10^{-8}$, against $1.1\times10^{-7}$ for a null change, nearly every opacity (all with $3 \le \log T \le 9$ and $-6 \le \log\rho \le 9.75$) times $1 + 10^{-10}$ on the same case), and so does the 0.8 $M_\odot$ case (testsuite/Test_m0p8feh+0p0GS98_rotp8p0d0p1mag1e4std_TAMS), except central $\log T$ at its last model (0.00047 dex, against 0.0002). The setting is not applied as documented, but neither the results nor the time spent change measurably. In CHECKJ's relaxed calls the fix's ACC would be $10^{-2}$. If convergence to machine precision is meant, the fix is a comment instead. Intent: The outer test's comment asks for the cube root (lines 138–139 ↗); nothing is said of the inner one.

Definitions

subroutineSHAPE

lines 2–170 · GitHub

Uses modules CONSTANTS_MOD, DIMENSIONS_MOD, NMLPAR_MOD

Called by FPFT, SOLID

Change log

The changes since rp_1

Changed since rp_1. Declarations only: COMMON blocks to modules (/CONST1/, /CONST2/ to CONSTANTS_MOD, /ROT/ to NMLPAR_MOD), the array size JSON from DIMENSIONS_MOD, IMPLICIT NONE.

1C $$$$$$ 2 SUBROUTINE SHAPE(HD,HR,HS,JSTART,JEND,OMEGA,ETA2,R0) 3C 10/15/26 MHP via Claude 4C code changed to take JSON from DIMENSIONS_MOD instead of declaring 5C PARAMETER(JSON=5000) here. 6 USE DIMENSIONS_MOD 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/20/26 MHP via Claude 14C code changed to replace implicit typing with IMPLICIT NONE and 15C explicit declarations. 16 IMPLICIT NONE 17 REAL*8 A,ACC,C1,C2,C3,C4,CG,DELR03,DETA1,DETA2,DETA3,DETA4,DR,ERR, 18 * * ETA2,ETA2T,ETAA,FACT,GM,HD,HR,HS,OMEGA,R0,R03,R0A,R0EST, 19 * * RHO,RHOA,RHOBA,RHOBAR,RHOBP,RHOP,RPHI,RPHI3,RPHIA,RPHIP 20 INTEGER I,J,JBEGIN,JEND,JSTART,K 21 22C 9/25/26 MHP via Claude 23C code changed to merge COMMON/CONST/, COMMON/CONST1/, COMMON/CONST2/, 24C and COMMON/CONST3/ into a single COMMON/CONSTANTS/ 25C 10/12/26 MHP via Claude 26C code changed to extend COMMON/NMLPAR/ to 234 members, adding the 27C namelist variables that the physics routines read from 27 COMMON 28C blocks which also held derived or working members. Those blocks keep 29C only their non-namelist members; five of them -- NEWMX, VMULT, 30C PMMWIND, CCOUT1 and ENVGEN -- had nothing left and are gone. 31C Derivation-only namelist variables from the same blocks went to 32C COMMON/NMLBYPASS/ instead, not here. 33C 34C Declared identically in all 121 files that hold any of these members. 35C Generated from one member list, so every copy is byte-identical; the 36C linker's 'size of symbol changed' warning is the check that they are. 37 DIMENSION HD(JSON),HR(JSON),HS(JSON),OMEGA(JSON),ETA2(JSON), 38 * *R0(JSON) 39 SAVE 40C SHAPE FINDS THE DISTORTION OF A GIVEN SHELL FROM SPHERICAL SYMMETRY; 41C IT CALCULATES THE CHARACTERISTIC RADIUS R0 AND ASSOCIATED PARAMETER 42C ETA2. IN POLAR CO-ORDINATES THE RADIUS OF A SHELL IS 43C R = R0*(1 - A*P2), WHERE A IS A FUNCTION OF ROTATION RATE AND P2 IS 44C THE SECOND LEGENDRE POLYNOMIAL. FOR MORE INFORMATION SEE ENDAL AND 45C SOFIA. 46 CG = DEXP(CLN*CGL) 47 C1 = 0.6D0 48 C2 = 2.0D0/3.5D1 49 C3 = 3.0D0*C1 50 C4 = 4.0D0*C2 51 JBEGIN = JSTART 52 IF(JSTART.EQ.1) THEN 53C CALCULATE FIRST POINT ETA2 USING CENTRAL B.C. 54C NOTE THAT THE SUFFIX 'P' DENOTES A QUANTITY FROM THE PREVIOUS SHELL. 55C E.G. RHOP IS THE DENSITY OF THE LAST SHELL AND RHO IS THE DENSITY OF 56C THE CURRENT ONE. 57 RHOBAR = DEXP(CLN*(HS(1) - C4PI3L - 3.0D0*HR(1))) 58 RPHI = DEXP(CLN*HR(1)) 59 RHO = DEXP(CLN*HD(1)) 60 ETA2(1) = 6.D0*(1.0D0 - RHO/RHOBAR) 61C ITERATE FOR R0 GIVEN RPHI AND ETA2 62 GM = CG*DEXP(CLN*HS(1)) 63 RPHI3 = RPHI**3 64 R03 = RPHI3 65 FACT = 5.0D0*CC13*OMEGA(1)**2/(GM*(2.0D0+ETA2(1))) 66 DO 10 J = 1,ITFP2 67 A = FACT*R03 68 DELR03 = (RPHI3-R03*(1.0D0 + C1*A**2 - C2*A**3))/ 69 * * (1.0D0 + C3*A**2 - C4*A**3) 70 R03 = R03 + DELR03 71 IF(DABS(DELR03/R03).LE.ACFPFT)GOTO 20 72 10 10 CONTINUE 73 20 20 R0(1) = R03**CC13 74 IF(JEND.EQ.1) GOTO 9999 75 RHOBP = RHOBAR 76 RPHIP = RPHI 77 RHOP = RHO 78 JBEGIN = 2 79 ELSE 80 RHOBP = DEXP(CLN*(HS(JSTART-1) - C4PI3L - 3.0D0*HR(JSTART-1))) 81 RPHIP = DEXP(CLN*HR(JSTART-1)) 82 RHOP = DEXP(CLN*HD(JSTART-1)) 83 ENDIF 84C CALCULATE ETA2 AND R0 FOR REMAINING POINTS 85 DO 100 I = JBEGIN,JEND 86 RHOBAR = DEXP(CLN*(HS(I) - C4PI3L - 3.0D0*HR(I))) 87 RHO = DEXP(CLN*HD(I)) 88 RPHI = DEXP(CLN*HR(I)) 89 GM = CG*DEXP(CLN*HS(I)) 90 DR = RPHI - RPHIP 91C FIND ETA2 USING 4-POINT RUNGE-KUTTE TECHNIQUE 92C D(ETA2)/D(R0) IS COMPUTED USING RADAU'S EQUATION: 93C R*(DETA2/DR)+6*RHO*(ETA2+1)/RHOBAR+ETA2*(ETA2-1) = 6 ,AND ETA2(0) = 0 94C FOR A FIRST GUESS,R0 = RPHI IS ASSUMED. 95C FIRST EVALUATE DETA/DR AT THE BEGINNING OF THE INTERVAL 96 DETA1 = DR*(6.0D0 - 6.0D0*RHOP*(ETA2(I-1) + 1.0D0)/RHOBP 97 * * - ETA2(I-1)*(ETA2(I-1) - 1.0D0))/RPHIP 98C RHOA AND RHOBA ARE AVERAGES OF THE RHO,RHOBAR OF OLD AND NEW SHELLS 99 RHOA = 0.5D0*(RHO+RHOP) 100 RHOBA = 0.5D0*(RHOBAR+RHOBP) 101 ETAA = ETA2(I-1) + 0.5D0*DETA1 102 RPHIA = RPHIP + 0.5D0*DR 103C USING THE ESTIMATED ETA2 AT THE MIDPOINT,FIND DETA/DR AT THE 104C MIDPOINT. 105 DETA2 = DR*(6.0D0 - 6.0D0*RHOA*(ETAA + 1.0D0)/RHOBA 106 * * - ETAA*(ETAA - 1.0D0))/RPHIA 107 ETAA = ETA2(I-1) + 0.5D0*DETA2 108C USING THIS REFINED DERIVATIVE TO ESTIMATE ETA2 AT MIDPOINT, 109C FIND DETA/DR AT THE MIDPOINT AGAIN. 110 DETA3 = DR*(6.0D0 - 6.0D0*RHOA*(ETAA + 1.0D0)/RHOBA 111 * * - ETAA*(ETAA - 1.0D0))/RPHIA 112 ETAA = ETA2(I-1) + DETA3 113C USING DETA/DR AT THE MIDPOINT TO ESTIMATE ETA2 AT THE END OF 114C THE INTERVAL,GET DETA/DR AT THE END OF THE INTERVAL. 115 DETA4 = DR*(6.0D0 - 6.0D0*RHO*(ETAA + 1.0D0)/RHOBAR 116 * * - ETAA*(ETAA - 1.0D0))/RPHI 117C PERFORM 4TH ORDER RUNGE-KUTTE INTEGRATION USING ABOVE 4 DERIVS. 118 ETA2(I) = ETA2(I-1)+CC13*(0.5D0*DETA1+DETA2+DETA3+0.5D0*DETA4) 119 RPHI3 = RPHI**3 120 R03 = RPHI3 121 R0EST = RPHI 122 ACC = ACFPFT**CC13 123C ITERATE BETWEEN SOLUTION FOR ETA2 AND SOLUTION FOR R0 ITFP1 TIMES. 124 DO 60 K = 1,ITFP1 125 FACT = 5.0D0*CC13*OMEGA(I)**2/(GM*(2.0D0+ETA2(I))) 126C NOW ITERATE FOR R0 GIVEN RPHI AND ETA2, USING THE RELATION 127C RPHI**3 = R0**3(1.0 + 3/5A**2 - 2/35A**3) 128C WHERE A = OMEGA**2*R0**3*5/3GM(2+ETA2)) 129 DO 40 J = 1,ITFP2 130 A = FACT*R03 131 DELR03 = (RPHI3-R03*(1.0D0 + C1*A**2 - C2*A**3))/ 132 * * (1.0D0 + C3*A**2 - C4*A**3) 133 R03 = R03 + DELR03 134 IF(DABS(DELR03/R03).LT.ACFPFT)GOTO 50 135 40 40 CONTINUE 136 50 50 R0(I) = R03**CC13 137 ERR = R0(I) - R0EST 138C ETA2 IS A FUNCTION OF R0, AND R0=RPHI WAS USED TO CALCULATE ETA2 139C CORRECT ETA2 HERE IF DR/R0 > CUBE ROOT OF ACFPFT 140 IF(DABS(ERR)/R0(I).LE.ACC) GOTO 70 141C FIND ETA2 USING 4-POINT RUNGE-KUTTE TECHNIQUE AGAIN,BUT FINDING 142C ETA2 AT R0(I) RATHER THAN ASSUMING R0 = RPHI. 143 ETA2T = ETA2(I) 144 DR = R0(I) - R0(I-1) 145 DETA1 = DR*(6.0D0 - 6.0D0*RHOP*(ETA2(I-1) + 1.0D0)/RHOBP 146 * * - ETA2(I-1)*(ETA2(I-1) - 1.0D0))/R0(I-1) 147 ETAA = ETA2(I-1) + 0.5D0*DETA1 148 R0A = R0(I-1) + 0.5D0*DR 149 DETA2 = DR*(6.0D0 - 6.0D0*RHOA*(ETAA + 1.0D0)/RHOBA 150 * * - ETAA*(ETAA - 1.0D0))/R0A 151 ETAA = ETA2(I-1) + 0.5D0*DETA2 152 DETA3 = DR*(6.0D0 - 6.0D0*RHOA*(ETAA + 1.0D0)/RHOBA 153 * * - ETAA*(ETAA - 1.0D0))/R0A 154 ETAA = ETA2(I-1) + DETA3 155 DETA4 = DR*(6.0D0 - 6.0D0*RHO*(ETAA + 1.0D0)/RHOBAR 156 * * - ETAA*(ETAA - 1.0D0))/R0(I) 157 ETA2(I) = ETA2(I-1) + CC13*(0.5D0*DETA1+DETA2+DETA3 158 * * + 0.5D0*DETA4) 159 ERR = ETA2T - ETA2(I) 160 IF(DABS(ERR).LE.ACC) GOTO 70 161 R0EST = R0(I) 162 60 60 CONTINUE 163 70 70 CONTINUE 164 RHOP = RHO 165 RHOBP = RHOBAR 166 RPHIP = RPHI 167 100 100 CONTINUE 168 9999 9999 CONTINUE 169 RETURN 170 END

1"""Python version of shape.f (YREC 7be4133): the distortion of each shell 2of a rotating star (Radau's eta_2 and the characteristic radius r_0).""" 3 4import math 5import numpy as np 6 7from pyrec import constants_mod 8from pyrec import nmlpar_mod 9 10 11def _deta(dr, rho, rhobar, eta, r): 12 """One Runge-Kutta increment of Radau's equation over dr, at radius r.""" 13 14 return dr * (6.0 - 6.0 * rho * (eta + 1.0) / rhobar 15 - eta * (eta - 1.0)) / r 16 17 18def shape(hd: np.ndarray, hr: np.ndarray, hs: np.ndarray, jstart: int, jend: int, omega: np.ndarray, eta2: np.ndarray, 19 r0: np.ndarray) -> None: 20 """For shells jstart..jend: integrate Radau's equation outward for 21 eta_2, solve for the characteristic radius r_0 (the shell's radius 22 being r_0 (1 - A P_2)) by Newton's method, and iterate the two. 23 24 Translates: shape.f#SHAPE (7be4133) 25 26 Arguments: hd, hr, hs (log rho, log r, log m, cgs); jstart, jend; 27 omega (angular velocity per shell); eta2 and r0, changed in place 28 (when jstart > 1, eta2[jstart-1] and r0[jstart-1] are read as the 29 starting values). None is returned. 30 31 """ 32 33 c = constants_mod 34 nml = nmlpar_mod 35 36 # The constants of the r_0 relation (Fortran lines 46–51). 37 cg = math.exp(c.cln * c.cgl) 38 c1 = 0.6 39 c2 = 2.0 / 3.5e1 40 c3 = 3.0 * c1 41 c4 = 4.0 * c2 42 jbegin = jstart 43 44 if jstart == 1: 45 # The first shell: eta_2 from its leading term near the centre 46 # (Fortran lines 52–60). 47 # Suspected bug (see the page): the factor 6 should be 6/7. 48 rhobar = math.exp(c.cln * (hs[1] - c.c4pi3l - 3.0 * hr[1])) 49 rphi = math.exp(c.cln * hr[1]) 50 rho = math.exp(c.cln * hd[1]) 51 eta2[1] = 6.0 * (1.0 - rho / rhobar) 52 53 # Newton's method for r_0 at the first shell (Fortran lines 61–73). 54 gm = cg * math.exp(c.cln * hs[1]) 55 rphi3 = rphi ** 3 56 r03 = rphi3 57 fact = 5.0 * c.cc13 * omega[1] ** 2 / (gm * (2.0 + eta2[1])) 58 # Question (see the page): silent non-convergence. 59 for j in range(1, nml.itfp2 + 1): 60 a = fact * r03 61 delr03 = ((rphi3 - r03 * (1.0 + c1 * a ** 2 - c2 * a ** 3)) 62 / (1.0 + c3 * a ** 2 - c4 * a ** 3)) 63 r03 = r03 + delr03 64 # Suspected bug (see the page): compared with ACFPFT, not its 65 # cube root. 66 if abs(delr03 / r03) <= nml.acfpft: 67 break 68 r0[1] = r03 ** c.cc13 # label 20 69 70 # The values at the previous shell (Fortran lines 74–83). 71 if jend == 1: 72 return 73 rhobp = rhobar 74 rphip = rphi 75 rhop = rho 76 jbegin = 2 77 else: 78 rhobp = math.exp(c.cln * (hs[jstart - 1] - c.c4pi3l 79 - 3.0 * hr[jstart - 1])) 80 rphip = math.exp(c.cln * hr[jstart - 1]) 81 rhop = math.exp(c.cln * hd[jstart - 1]) 82 83 for i in range(jbegin, jend + 1): 84 85 # The shell's values (Fortran lines 86–89). 86 rhobar = math.exp(c.cln * (hs[i] - c.c4pi3l - 3.0 * hr[i])) 87 rho = math.exp(c.cln * hd[i]) 88 rphi = math.exp(c.cln * hr[i]) 89 gm = cg * math.exp(c.cln * hs[i]) 90 91 # Runge-Kutta step in the volume radius (Fortran lines 90–118). 92 dr = rphi - rphip 93 deta1 = _deta(dr, rhop, rhobp, eta2[i - 1], rphip) 94 rhoa = 0.5 * (rho + rhop) 95 rhoba = 0.5 * (rhobar + rhobp) 96 etaa = eta2[i - 1] + 0.5 * deta1 97 rphia = rphip + 0.5 * dr 98 deta2 = _deta(dr, rhoa, rhoba, etaa, rphia) 99 etaa = eta2[i - 1] + 0.5 * deta2 100 deta3 = _deta(dr, rhoa, rhoba, etaa, rphia) 101 etaa = eta2[i - 1] + deta3 102 deta4 = _deta(dr, rho, rhobar, etaa, rphi) 103 eta2[i] = eta2[i - 1] + c.cc13 * (0.5 * deta1 + deta2 + deta3 104 + 0.5 * deta4) 105 106 # Iterate between eta_2 and r_0 (Fortran lines 119–163). 107 rphi3 = rphi ** 3 108 r03 = rphi3 109 r0est = rphi 110 acc = nml.acfpft ** c.cc13 111 for k in range(1, nml.itfp1 + 1): 112 fact = 5.0 * c.cc13 * omega[i] ** 2 / (gm * (2.0 + eta2[i])) 113 for j in range(1, nml.itfp2 + 1): 114 a = fact * r03 115 delr03 = ((rphi3 - r03 * (1.0 + c1 * a ** 2 - c2 * a ** 3)) 116 / (1.0 + c3 * a ** 2 - c4 * a ** 3)) 117 r03 = r03 + delr03 118 # Suspected bug (see the page): compared with ACFPFT, not 119 # its cube root. 120 if abs(delr03 / r03) < nml.acfpft: 121 break 122 r0[i] = r03 ** c.cc13 # label 50 123 err = r0[i] - r0est 124 # Question (see the page): the two tolerance tests use one ACC. 125 if abs(err) / r0[i] <= acc: 126 break 127 128 # Runge-Kutta step again, in r_0 (Fortran lines 141–162). 129 # Question (see the page): RHOA, RHOBA are the first pass's. 130 eta2t = eta2[i] 131 dr = r0[i] - r0[i - 1] 132 deta1 = _deta(dr, rhop, rhobp, eta2[i - 1], r0[i - 1]) 133 etaa = eta2[i - 1] + 0.5 * deta1 134 r0a = r0[i - 1] + 0.5 * dr 135 deta2 = _deta(dr, rhoa, rhoba, etaa, r0a) 136 etaa = eta2[i - 1] + 0.5 * deta2 137 deta3 = _deta(dr, rhoa, rhoba, etaa, r0a) 138 etaa = eta2[i - 1] + deta3 139 deta4 = _deta(dr, rho, rhobar, etaa, r0[i]) 140 eta2[i] = eta2[i - 1] + c.cc13 * (0.5 * deta1 + deta2 + deta3 141 + 0.5 * deta4) 142 err = eta2t - eta2[i] 143 if abs(err) <= acc: 144 break 145 r0est = r0[i] 146 147 # label 70: the values for the next shell (Fortran lines 163–166). 148 rhop = rho 149 rhobp = rhobar 150 rphip = rphi

shape.f

Purpose

How it fits

The method

From YREC's documentation

Suspected bugs, checked by running the code

1. (harmless in the runs) The first shell's $\eta_2$ is seven times too large.

2. (harmless in the runs) The Newton test for $r_0$ compares with ACFPFT, not its cube root.

Definitions

subroutineSHAPE

Change log

2. (harmless in the runs) The Newton test for $r_0$ compares with `ACFPFT`, not its cube root.