diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index a635a9b4b..bda5dec0b 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -578,8 +578,8 @@ void BaseStar::CalculateAnCoefficients(DBL_VECTOR &p_AnCoefficients, RConstants(C_BETA_R) = (a[69] * 16384.0) / (a[70] + PPOW(16.0, a[71])); // Hurley et al. 2000, eq 22a RConstants(B_DELTA_R) = (a[38] + a[39] * 8.0 * M_SQRT2) / (a[40] * 8.0 + PPOW(2.0, a[41])) - 1.0; // Hurley et al. 2000, eq 17 - GammaConstants(B_GAMMA) = a[76] + (a[77] * PPOW((1.0 - a[78]), a[79])); // Hurley et al. 2000, eq 23 - GammaConstants(C_GAMMA) = (utils::Compare(a[75], 1.0) == 0) ? GammaConstants(B_GAMMA) : a[80]; // Hurley et al. 2000, eq 23 + GammaConstants(B_GAMMA) = max(0.0, a[76] + (a[77] * PPOW((1.0 - a[78]), a[79]))); // Hurley et al. 2000, eq 23 and discussion immediately following - max() confirmed in BSE Fortran code + 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 #undef GammaConstants #undef RConstants diff --git a/src/MainSequence.cpp b/src/MainSequence.cpp index b06640564..f8546ff92 100644 --- a/src/MainSequence.cpp +++ b/src/MainSequence.cpp @@ -228,14 +228,14 @@ double MainSequence::CalculateGamma(const double p_Mass) const { #define B_GAMMA m_GammaConstants[static_cast(GAMMA_CONSTANTS::B_GAMMA)] // for convenience and readability - undefined at end of function #define C_GAMMA m_GammaConstants[static_cast(GAMMA_CONSTANTS::C_GAMMA)] // for convenience and readability - undefined at end of function - double gamma; +double gamma = 0.0; // default return value - if (utils::Compare(p_Mass, 1.0) <= 0) gamma = a[76] + (a[77] * PPOW(p_Mass - a[78], a[79])); - else if (utils::Compare(p_Mass, a[75]) <= 0) gamma = B_GAMMA + (a[80] - B_GAMMA) * PPOW((p_Mass - 1.0) / (a[75] - 1.0), a[81]); - else if (utils::Compare(p_Mass, (a[75] + 0.1)) <= 0) gamma = C_GAMMA - (10.0 * (p_Mass - a[75]) * C_GAMMA); // included = case, missing from Hurley+ 2000 - else gamma = 0.0; // this really is zero + if (utils::Compare(p_Mass, 1.0) <= 0) gamma = a[76] + (a[77] * PPOW(std::abs(p_Mass - a[78]), a[79])); // BSE Fortran code has abs() +else if (utils::Compare(p_Mass, a[75]) <= 0) gamma = B_GAMMA + (a[80] - B_GAMMA) * PPOW((p_Mass - 1.0) / (a[75] - 1.0), a[81]); +else if (utils::Compare(p_Mass, (a[75] + 0.1)) <= 0) gamma = C_GAMMA - (10.0 * (p_Mass - a[75]) * C_GAMMA); // see discussion just prior to eq 23 - the end point is wrong in the arxiv version of Hurley et al. 2000 (should be 0.1, not 1.0) - confirmed in BSE Fortran code +else gamma = 0.0; // see discussion just prior to eq 23 - confirmed in BSE Fortran code - return gamma; +return std::max(0.0, gamma); // see discussion following eq 23 - confirmed in BSE Fortran code #undef C_GAMMA #undef B_GAMMA diff --git a/src/changelog.h b/src/changelog.h index 175768c45..026aa60d2 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1659,6 +1659,11 @@ // - Corrected calculations for Hurley A(n) and B(n) coefficients (see Hurley et al. 2000, appendix) // - Changed utils::GetGSLVersion() to avoid compiler warning "warning: ignoring attributes on template argument ‘int (*)(FILE*)’ [-Wignored-attributes]" // - Reverted Makefile line "SOURCES := $(wildcard *.cpp)" to listing actual source files +// 03.25.02 JR - August 20, 2025 - Defect repairs: +// - 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) +// - Added abs() to gamma calculation in Mainsequence.cpp::CalculateGamma() (per BSE Fortran source) +// - Clamped gamma to [0.0, gamma] in Mainsequence.cpp::CalculateGamma() (per discussion just after eq 23 - confirmed in BSE Fortran source) // // Version string format is MM.mm.rr, where // @@ -1669,7 +1674,7 @@ // if MM is incremented, set mm and rr to 00, even if defect repairs and minor enhancements were also made // if mm is incremented, set rr to 00, even if defect repairs were also made -const std::string VERSION_STRING = "03.25.01"; +const std::string VERSION_STRING = "03.25.02"; # endif // __changelog_h__