Defect repairs

[jeffriley]

@ilyamandel a couple more minor fixes related to Hurley coefficients and constants:

• Clamped B_GAMMA to [0.0, B_GAMMA] (per discussion just after eq 23 - confirmed in BSE Fortran source)
• Corrected calculation for Hurley Gamma constant C (C_GAMMA - see Hurley et al. 2000, just after eq 23, should use a(75) <= 1.0, not a(75) == 1.0 - confirmed in BSE Fortran source)

[ilyamandel]

Can I suggest that we reference the Fortran code in the comments along with Hurley Eq. (23) -- e.g., "// see Hurley et al. 2000, eq 23, but the code here follows the Fortran SSE implementation". Because the previous version implemented Eq. 23 as written in the paper, the new one does not (but I trust you that it implements the Fortran code).

[jeffriley]

Because the previous version implemented Eq. 23 as written in the paper, the new one does not

Do you mean the addition of the max()? That's in the discussion following eq 23 (as is the <= 1.0)

Edit: When I first made the change I added a max() to C_GAMMA (to ensure it is >= 0.0), but I wasn't sure the comment didn't refer only to B = gamma (when M = 1.0)... Checking MainSequence::CalculateGamma(const double p_Mass) now, we don't ensure gamma is >= 0.0 - I'm pretty sure Hurley does in the BSE code - I'll confirm that tomorrow.

[jeffriley]

Here's the state of affairs:

Current COMPAS code: Mainsequence.cpp::CalculateGamma(), BaseStar::CalculateAnCoefficients() (GammaConstants(B_GAMMA) and GammaConstants(C_GAMMA))
Legacy code: star.cpp, calculate_gamma(), calculate_B_gamma(), calculate_C_gamma()

BSE code zfuncs.f:

*
* A function to evaluate the radius gamma coefficent.
* (JH 24/11/97)
*
m1 = 1.d0
if(m.gt.(a(88)+0.1d0))then                                            <-- see discussion immediately prior to eq 23; endpoint is 0.1, not 1.0
    rgammf = 0.d0                                                     <-- gamma = 0.0
else
    b1 = MAX(0.d0,a(83) + a(84)*(m1-a(85))**a(86))                    <-- this is B (in paper) (B_GAMMA in COMPAS) (but is also used as C (in paper, C_GAMMA in COMPAS)).  MAX() missing in legacy code and current COMPAS code
    if(m.le.m1)then                                                   <-- eq 23 condition M <= 1.0
        rgammf = a(83) + a(84)*ABS(m-a(85))**a(86)                    <-- gamma - ABS() missing in legacy code and current COMPAS code
    elseif(m.le.a(88))then                                            <-- eq 23 condition 1.0 < M <= a[75]
        rgammf = b1 + (a(89) - b1)*((m - m1)/(a(88) - m1))**a(87)     <-- gamma
    else                                                              <-- eq 23 condition a[75] < M < a[75] + 0.1 (note 0.1, not 1.0)
        if(a(88).gt.m1) b1 = a(89)                                    <-- C = a[80] (C_GAMMA in COMPAS) - see discussion following eq 23; note reuse of b1 as C (was B)
        rgammf = b1 - 10.d0*b1*(m - a(88))                            <-- gamma
    endif
    rgammf = MAX(rgammf,0.d0)                                         <-- ensure gamma +ve.  MAX() missing in current COMPAS code; exists in legacy code
endif

return
end
    

BSE code a(83) = paper a[76]
BSE code a(84) = paper a[77]
BSE code a(85) = paper a[78]
BSE code a(86) = paper a[79]
BSE code a(87) = paper a[81]
BSE code a(88) = paper a[75]
BSE code a(89) = paper a[80]

So, I will push another change that adds ABS() appropriately, and clamps gamma to [0.0, gamma], in Mainsequence.cpp::CalculateGamma()

[ilyamandel]

GammaConstants(C_GAMMA) = (utils::Compare(a[75], 1.0) <= 0) ? GammaConstants(B_GAMMA) : a[80]; // Hurley et al. 2000, eq 23 and discussion immediately following - <= 1.0 confirmed in BSE Fortran code

But the Hurley paper actually says
"and C=a80 unless a75=1.0 when C=B"
Equal, not \leq.
That's why I was suggesting referencing the Fortran code directly in opposition to the paper -- but I am OK with the current comment wording if you prefer that.

[jeffriley]

"and C=a80 unless a75=1.0 when C=B"

Not the version I'm looking at...

hurleyp9.pdf