Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010() and CH::CalculateMassLossRateFractionOB() [previously CalculateMassLossRateWeightOB()]

Author: Ilya Mandel

src/CH.cpp

@@ -286,36 +286,33 @@ void CH::UpdateAgeAfterMassLoss() {
 
 
 /*
- * Calculate the weight for the OB and WR mass loss prescriptions
+ * CalculateMassLossFractionOB
  *
- * According to the prescription described in Yoon et al. 2006 (and Szecsi et al. 2015)
- * Use OB mass loss rates until Y1 = 0.55, WR mass loss rates above Y2 = 0.7, and 
- * linearly interpolate for Y1 < Ys < Y2 where Ys is the surface helium fraction.
- * 
- * Since we don't have Ys by default in our models, and detailed models show Ys
- * rises ~linearly from 0.2 to 1.0 across the main-sequence, here we simply 
- * use tau = t / tMS as a proxy for Ys.
+ * @brief
+ * Calculate the fraction of mass loss attributable to OB mass loss, per Yoon et al. 2006
+ *
+ * The model described in Yoon et al. 2006 (also Szecsi et al. 2015) uses OB mass loss while the
+ * He surface abundance is below 0.55, WR mass loss when the surface He abundance is above 0.7,
+ * and linearly interpolate when the He surface abundance is between those limits.
+ *
+ * This function calculates the fraction of mass loss attributable to OB mass loss, based on
+ * the He surface abundance and the abundance limits described in Yoon et al. 2006.  The value
+ * returned will be 1.0 if 100% of the mass loss is attributable to OB mass lass, 0.0 if 100% of
+ * the mass loss is attributable to WR mass loss, and in the range (0.0, 1.0) if the mass loss is
+ * a mix of OB and WR.
  *
- * CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface)
  *
- * @param   [IN]        p_HeliumAbundanceSurface           Surface helium abundance Ys
+ * double CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const
  *
- * @return                                   weight for OB mass loss rate
+ * @param       p_HeAbundanceSurface            Helium abundance at the surface of the star
+ * @return                                      Fraction of mass loss attributable to OB mass loss
  */
-double CH::CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface) {
+double CH::CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const {
 
-    const double y1 = 0.55;
-    const double y2 = 0.70;
-    double       Ys = p_HeliumAbundanceSurface; // Get surface helium abundance
-    double weight   = 0.0;
+    constexpr double limOB = 0.55;                                          // per Yoon et al. 2006
+    constexpr double limWR = 0.70;                                          // per Yoon et al. 2006
 
-
-    // Calculate the weight as a function of Ys according to the prescription from Yoon et al. 2006
-         if (Ys < y1) weight = 1.0;
-    else if (Ys > y2) weight = 0.0;
-    else              weight = (y2 - Ys) / y2 - y1;
-
-    return weight;
+    return std::min(1.0, std::max (0.0, (limWR - p_HeAbundanceSurface) / (limWR - limOB)));
 }
 
 
@@ -337,7 +334,7 @@ double CH::CalculateMassLossRateBelczynski2010() {
     double Mdot   = 0.0;
     double MdotOB = 0.0;
     double MdotWR = 0.0;
-    double weight = 1.0;    // Initialised to 1.0 to allow us to use the OB mass loss rate by default 
+    double fractionOB = 1.0;    // Initialised to 1.0 to allow us to use the OB mass loss rate by default
 
     // Calculate OB mass loss rate according to the Vink et al. formalism
     MdotOB = BaseStar::CalculateMassLossRateBelczynski2010();
@@ -348,16 +345,16 @@ double CH::CalculateMassLossRateBelczynski2010() {
         // Calculate WR mass loss rate
         MdotWR = BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0);
 
-        // Calculate weight for combining these into total mass-loss rate
-        weight = CalculateMassLossRateWeightOB(m_HeliumAbundanceSurface);
+        // Calculate fraction for combining these into total mass-loss rate
+        fractionOB = CalculateMassLossRateFractionOB(m_HeliumAbundanceSurface);
 
     }
 
+    // Finally, combine each of these prescriptions according to the OB wind fraction
+    Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR);
+    
     // Set dominant mass loss rate
-    m_DominantMassLossRate = weight == 0 ? MASS_LOSS_TYPE::WR : MASS_LOSS_TYPE::OB;   
-
-    // Finally, combine each of these prescriptions according to the weight
-    Mdot = (weight * MdotOB) + ((1.0 - weight) * MdotWR);
+    m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
 
     // Enhance mass loss rate due to rotation
     Mdot *= CalculateMassLossRateEnhancementRotation();

src/CH.h

@@ -74,7 +74,7 @@ class CH: virtual public BaseStar, public MS_gt_07 {
     // Mass loss rate
     double          CalculateMassLossRateBelczynski2010();
     double          CalculateMassLossRateMerritt2025();
-    double          CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface);
+    double          CalculateMassLossRateFractionOB(const double p_HeliumAbundanceSurface);
 
     // Radius
     double          CalculateRadiusOnPhase() const                      { return m_RZAMS; }                                                                                         // Constant from birth

src/changelog.h

@@ -1664,8 +1664,9 @@
 //                                          - 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)
-//  03.26.00  IM - September 2, 2025    - Enhancement:
+//  03.26.00  IM - September 2, 2025    - Enhancement, defect repairs:
 //                                          - First (simplified) implementation of the Lau+ (2024) Hamstars thermally limited accretion prescription
+//                                          - Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010() and CH::CalculateMassLossRateFractionOB() [previously CalculateMassLossRateWeightOB()]
 //
 // Version string format is MM.mm.rr, where
 //