diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index d742bea7b..624fa9f93 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -2086,6 +2086,8 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { // no m_MassTransferTrackerHistory = m_Donor == m_Star1 ? MT_TRACKING::STABLE_1_TO_2_SURV : MT_TRACKING::STABLE_2_TO_1_SURV; // record what happened - for later printing double massGainAccretor = -massDiffDonor * m_FractionAccreted; // set accretor mass gain to mass loss * conservativeness + + double omegaDonor_pre_MT = m_Donor->Omega(); // used if full donor envelope is removed m_Donor->SetMassTransferDiffAndResolveWDShellChange(massDiffDonor); // set new mass of donor m_Accretor->SetMassTransferDiffAndResolveWDShellChange(massGainAccretor); // set new mass of accretor @@ -2095,7 +2097,12 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) { m_aMassTransferDiff = aFinal - aInitial; // set change in orbit (semi-major axis) STELLAR_TYPE stellarTypeDonor = m_Donor->StellarType(); // donor stellar type before resolving envelope loss - if (isEnvelopeRemoved) m_Donor->ResolveEnvelopeLossAndSwitch(); // if this was an envelope stripping episode, resolve envelope loss + + if (isEnvelopeRemoved) { // if this was an envelope stripping episode, resolve envelope loss + m_Donor->ResolveEnvelopeLossAndSwitch(); // resolve envelope loss for the donor and switch to new stellar type + m_Donor->SetOmega(omegaDonor_pre_MT); // keep the rotation frequency of the core equal to the pre-envelope-loss rotation frequency + } + if (m_Donor->StellarType() != stellarTypeDonor) { // stellar type change? (void)PrintDetailedOutput(m_Id, BSE_DETAILED_RECORD_TYPE::STELLAR_TYPE_CHANGE_DURING_MT); // yes - print (log) detailed output } @@ -2437,6 +2444,7 @@ void BaseBinaryStar::ResolveMassChanges() { if (utils::Compare(m_Star1->MassPrev(), m_Star1->Mass()) == 0) { // mass already updated? // no - resolve mass changes double massChange = m_Star1->MassLossDiff() + m_Star1->MassTransferDiff(); // mass change due to winds and mass transfer + if (utils::Compare(massChange, 0.0) != 0) { // winds/mass transfer changes mass? // yes - calculate new angular momentum; assume accretor is adding angular momentum from a circular orbit at the stellar radius double angularMomentumChange = (utils::Compare(massChange, 0.0) > 0) ? diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index ded560558..fa0273cef 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -2844,7 +2844,7 @@ void BaseStar::ResolveMassLoss(const bool p_UpdateMDt) { UpdateInitialMass(); // update effective initial mass (MS, HG & HeMS) UpdateAgeAfterMassLoss(); // update age (MS, HG & HeMS) ApplyMassTransferRejuvenationFactor(); // apply age rejuvenation factor - SetAngularMomentum(m_AngularMomentum + angularMomentumChange); + SetAngularMomentum(m_AngularMomentum + angularMomentumChange); } } diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 1602a0b31..fb48a2e0b 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -770,10 +770,12 @@ double GiantBranch::CalculateCoreMassAtBGB(const double p_Mass, const DBL_VECTOR #define gbParams(x) p_GBParams[static_cast(GBP::x)] // for convenience and readability - undefined at end of function #define massCutoffs(x) m_MassCutoffs[static_cast(MASS_CUTOFF::x)] // for convenience and readability - undefined at end of function + if (utils::Compare(p_Mass, massCutoffs(MHeF)) <= 0) return 0.0; // No McBGB for stars with mass below the helium flash threshold, see text above Eq. (44) of Hurley+ (2000) + double luminosity = GiantBranch::CalculateLuminosityAtPhaseBase_Static(massCutoffs(MHeF), m_AnCoefficients); double Mc_MHeF = BaseStar::CalculateCoreMassGivenLuminosity_Static(luminosity, p_GBParams); double c = (Mc_MHeF * Mc_MHeF * Mc_MHeF * Mc_MHeF) - (MC_L_C1 * PPOW(massCutoffs(MHeF), MC_L_C2)); // pow() is slow - use multiplication - + return std::min((0.95 * gbParams(McBAGB)), std::sqrt(std::sqrt(c + (MC_L_C1 * PPOW(p_Mass, MC_L_C2))))); // sqrt is much faster than PPOW() #undef massCutoffs diff --git a/src/HeHG.cpp b/src/HeHG.cpp index 34bfbd485..59ecb63f8 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -114,7 +114,7 @@ void HeHG::CalculateGBParams_Static(const double p_Mass0, DBL_VECTOR &p_GBParams) { #define gbParams(x) p_GBParams[static_cast(GBP::x)] // for convenience and readability - undefined at end of function - + GiantBranch::CalculateGBParams_Static(p_Mass, p_LogMetallicityXi, p_MassCutoffs, p_AnCoefficients, p_BnCoefficients, p_GBParams); // calculate common values (actually, all) // recalculate HeHG specific values @@ -125,6 +125,12 @@ void HeHG::CalculateGBParams_Static(const double p_Mass0, gbParams(Mx) = GiantBranch::CalculateCoreMass_Luminosity_Mx_Static(p_GBParams); // depends on B, D, p & q - recalculate if any of those are changed gbParams(Lx) = GiantBranch::CalculateCoreMass_Luminosity_Lx_Static(p_GBParams); // depends on B, D, p, q & Mx - recalculate if any of those are changed + // need to reset p and q to HeHG specific versions -- while they are set in GiantBranch::CalculateGBParams_Static(), that uses + // the GiantBranch versions of CalculateCoreMass_Luminosity_p_Static and CalculateCoreMass_Luminosity_q_Static + // should return to this and understand desired behavior - *ILYA* + gbParams(p) = CalculateCoreMass_Luminosity_p_Static(p_Mass, p_MassCutoffs); + gbParams(q) = CalculateCoreMass_Luminosity_q_Static(p_Mass, p_MassCutoffs); + gbParams(McBAGB) = p_Mass0; gbParams(McBGB) = GiantBranch::CalculateCoreMassAtBGB_Static(p_Mass, p_MassCutoffs, p_AnCoefficients, p_GBParams); diff --git a/src/WhiteDwarfs.cpp b/src/WhiteDwarfs.cpp index 3e49e3b07..494647678 100644 --- a/src/WhiteDwarfs.cpp +++ b/src/WhiteDwarfs.cpp @@ -155,6 +155,8 @@ double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) { // sanity check for mass - just return 0.0 if mass <= 0 if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0; + if (utils::Compare(p_Mass, MCH) >= 0) return NEUTRON_STAR_RADIUS; // only expected to come up if asking for the core or remnant radius of a giant star + double MCH_Mass_one_third = std::cbrt(MCH / p_Mass); double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third; diff --git a/src/changelog.h b/src/changelog.h index 7aedd9410..f6e40f14f 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1411,6 +1411,11 @@ // - Minor fixes, including to #1255, #1258 // 03.10.00 JR - Nov 29, 2024 - Enhancement: // - added functionality to allow stellar mergers (for BSE) to be logged to switchlog file (see documentation for details) -const std::string VERSION_STRING = "03.10.00"; +// 03.10.01 IM - Nov 30, 2024 - Defect repair: +// - corrected treatment of rotation to retain pre-mass-loss spin frequency, not angular momentum, on complete envelope removal during stable mass transfer +// - fixed issue with updating helium giants that manifested as supernovae with nan core mass (see #1245) +// - added check for exceeding Chandrasekhar mass when computing white dwarf radius (resolves issue #1264) +// - added check to only compute McBGB for stars with mass above MHeF, following text above Eq. 44 in Hurley+, 2000 (resolves issue #1256) +const std::string VERSION_STRING = "03.10.01"; # endif // __changelog_h__