Recalculate timescales when updating stellar age after mass loss

src/BaseBinaryStar.cpp

@@ -1995,6 +1995,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
     double zetaEquilibrium = m_Donor->CalculateZetaEquilibrium();
 
     m_ZetaLobe = CalculateZetaRocheLobe(jLoss, betaNuclear);                                                                    // try nuclear timescale mass transfer first
+
     if (m_Donor->IsOneOf(ALL_MAIN_SEQUENCE) && utils::Compare(zetaEquilibrium, m_ZetaLobe) > 0) {
         m_MassLossRateInRLOF    = donorMassLossRateNuclear;
         m_FractionAccreted      = betaNuclear;
@@ -2097,6 +2098,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
         m_Donor->SetRLOFOntoNS();                                                                                               // donor donated mass to a neutron star
         m_Accretor->SetRecycledNS();                                                                                            // accretor is (was) a recycled NS
 	}
+
 }
 
 

src/BaseStar.h

@@ -325,7 +325,7 @@ class BaseStar {
     virtual double          ResolveCommonEnvelopeAccretion(const double p_FinalMass,
                                                            const double p_CompanionMass     = 0.0,
                                                            const double p_CompanionRadius   = 0.0,
-                                                           const double p_CompanionEnvelope = 0.0)              { return p_FinalMass - Mass(); }                                    // LvS: todo: more consistent super eddington accretion during CE should also affect e.g. MS stars
+                                                           const double p_CompanionEnvelope = 0.0)              { return p_FinalMass - Mass(); }                                    // Overwritten in NS.h; for now, no accretion on stars other than compact objects during CE
 
     virtual STELLAR_TYPE    ResolveEnvelopeLoss(bool p_Force = false)                                           { return m_StellarType; }
 

src/CH.cpp

@@ -277,10 +277,11 @@ void CH::UpdateAgeAfterMassLoss() {
     // Calculate the ratio of the lifetime of a CH star to a normal MS star
     double lifetimeRatio = CalculateLifetimeRatio(m_Mass);
 
-    tBGBprime *= lifetimeRatio;
     tMSprime  *= lifetimeRatio;
 
     m_Age *= tMSprime / tMS;
+
+    CalculateTimescales(m_Mass, m_Timescales);
 }
 
 

src/HG.cpp

@@ -1017,6 +1017,8 @@ void HG::UpdateAgeAfterMassLoss() {
     double tMSprime  = MainSequence::CalculateLifetimeOnPhase(m_Mass0, tBGBprime);
 
     m_Age = tMSprime + (((tBGBprime - tMSprime) / (tBGB - tMS)) * (m_Age - tMS));
+
+    CalculateTimescales(m_Mass0, m_Timescales);
 }
 
 

src/HeGB.h

@@ -52,14 +52,14 @@ class HeGB: virtual public BaseStar, public HeHG {
     void Initialise(const STELLAR_TYPE p_PreviousStellarType) {
         CalculateTimescales();                                                                                                                                          // Initialise timescales
         if (p_PreviousStellarType != STELLAR_TYPE::NAKED_HELIUM_STAR_HERTZSPRUNG_GAP)                                                                                   // If not evolving from HeHG...
-            m_Age = CalculateAgeOnPhase_Static(m_Mass, m_COCoreMass, m_Timescales[static_cast<int>(TIMESCALE::tHeMS)], m_GBParams);                                     // ... Set age appropriately
+            m_Age = CalculateAgeOnPhase_Static(m_Mass, m_COCoreMass, m_Timescales[static_cast<int>(TIMESCALE::tHeMS)], m_GBParams);                                 // ... Set age appropriately
         EvolveOnPhase(0.0);
     }
 
 
     // member functions - alphabetically
     double      CalculateCriticalMassRatioClaeys14(const bool p_AccretorIsDegenerate) const ;
-    double      CalculateCriticalMassRatioHurleyHjellmingWebbink() const                                            { return 1.28; }                                    // From BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
+    double      CalculateCriticalMassRatioHurleyHjellmingWebbink() const                                            { return 1.28; }                                        // From BSE. Using the inverse owing to how qCrit is defined in COMPAS. See Hurley et al. 2002 sect. 2.6.1 for additional details.
 
     double      CalculateLuminosityOnPhase(const double p_CoreMass, const double p_GBPB, const double p_GBPD) const { return CalculateLuminosityOnPhase_Static(p_CoreMass, p_GBPB, p_GBPD); }
     double      CalculateLuminosityOnPhase() const                                                                  { return CalculateLuminosityOnPhase(m_CoreMass, m_GBParams[static_cast<int>(GBP::B)], m_GBParams[static_cast<int>(GBP::D)]); }
@@ -70,7 +70,7 @@ class HeGB: virtual public BaseStar, public HeHG {
     std::tuple <double, STELLAR_TYPE> CalculateRadiusAndStellarTypeOnPhase(const double p_Mass, const double p_Luminosity) const;
     std::tuple <double, STELLAR_TYPE> CalculateRadiusAndStellarTypeOnPhase() const                                  { return CalculateRadiusAndStellarTypeOnPhase(m_Mass, m_Luminosity); }
 
-    ENVELOPE    DetermineEnvelopeType() const                                                                       { return ENVELOPE::CONVECTIVE; }                    // Always CONVECTIVE
+    ENVELOPE    DetermineEnvelopeType() const                                                                       { return ENVELOPE::CONVECTIVE; }                        // Always CONVECTIVE    
 };
 
 #endif // __HeGB_h__

src/HeHG.h

@@ -128,7 +128,7 @@ class HeHG: virtual public BaseStar, public HeMS {
             bool            ShouldEvolveOnPhase() const;
             bool            ShouldSkipPhase() const                                                                 { return false; }                                                       // Never skip HeMS phase
 
-            void            UpdateAgeAfterMassLoss()                                                                { GiantBranch::UpdateAgeAfterMassLoss(); }                              // Skip HeMS
+            void            UpdateAgeAfterMassLoss()                                                                { GiantBranch::UpdateAgeAfterMassLoss(); }                              // No action for He giants
             void            UpdateInitialMass()                                                                     { GiantBranch::UpdateInitialMass(); }                                   // Skip HeMS
 };
 

src/HeMS.cpp

@@ -530,6 +530,8 @@ void HeMS::UpdateAgeAfterMassLoss() {
     double tHeMSprime = CalculateLifetimeOnPhase_Static(m_Mass);
 
     m_Age *= tHeMSprime / tHeMS;
+
+    CalculateTimescales(m_Mass, m_Timescales);
 }
 
 

src/MainSequence.cpp

@@ -747,6 +747,7 @@ void MainSequence::UpdateAgeAfterMassLoss() {
     double tMSprime  = MainSequence::CalculateLifetimeOnPhase(m_Mass, tBGBprime);
 
     m_Age *= tMSprime / tMS;
+    CalculateTimescales(m_Mass, m_Timescales);                                      // must update timescales
 }
 
 ///////////////////////////////////////////////////////////////////////////////////////

src/changelog.h

@@ -1391,8 +1391,10 @@
 //                                      - Added new critical mass ratio tables for He stars from Ge et al. team
 //                                      - Cleaned up the stability calculation for H-rich stars as well, specifically implementing nearest neighbor for extrapolation
 //                                      - Now all of their results from Papers I-V are included (including those requested in private comm.)
-// 03.08.03   VK - Nov 20, 2024     - Defect repairs
-//                                      - Fixed behavior for core spin to be retained after envelope loss.
-const std::string VERSION_STRING = "03.08.03";
+// 03.08.03   VK - Nov 20, 2024     - Defect repair:
+//                                      - Fixed behavior for core spin to be retained after envelope loss
+// 03.08.04   IM - Nov 25, 2024     - Defect repair:
+//                                      - Recalculate timescales when updating stellar age after mass loss (addresses issue #1231)
+const std::string VERSION_STRING = "03.08.04";
 
 # endif // __changelog_h__