TeamCOMPAS · ilyamandel · Dec 9, 2025 · Dec 4, 2025 · Dec 9, 2025
diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml
@@ -1,5 +1,5 @@
 ##~!!~## COMPAS option values
-##~!!~## File Created Mon Oct 27 18:34:22 2025 by COMPAS v03.27.01
+##~!!~## File Created Tue Dec  9 15:41:35 2025 by COMPAS v03.27.02
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -29,7 +29,7 @@ booleanChoices:
 #    --enhance-CHE-lifetimes-luminosities: True                            # Default: True
 #    --expel-convective-envelope-above-luminosity-threshold: False         # Default: False
 #    --natal-kick-for-PPISN: False                                         # Default: False
-#    --scale-mass-loss-with-surface-helium-abundance: True                 # Default: True
+#    --scale-CHE-mass-loss-with-surface-helium-abundance: True             # Default: True
 
     ### BINARY PROPERTIES
 #    --allow-touching-at-birth: False                                      # Default: False                                                                                                                                                                            # record binaries that have stars touching at birth in output files

diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst
@@ -1319,9 +1319,9 @@ Default = DECIN2023 |br|
 
 :ref:`Back to Top <options-props-top>`
 
-**--scale-mass-loss-with-surface-helium-abundance** |br|
-Scale mass loss for main sequence, including chemically homogeneously evolving (CHE), stars with the surface helium abundance. 
-Transition from OB/VMS to WR mass loss towards the end of the main sequence.
+**--scale-CHE-mass-loss-with-surface-helium-abundance** |br|
+Scale mass loss for chemically homogeneously evolving (CHE) stars with the surface helium abundance. 
+Transition from OB to WR mass loss towards the end of the main sequence.
 Default = TRUE
 
 **--scale-terminal-wind-velocity-with-metallicity-power** |br|
@@ -1535,7 +1535,7 @@ Go to :ref:`the top of this page <options-props-top>` for the full alphabetical
 
 --check-photon-tiring-limit, --cool-wind-mass-loss-multiplier, --luminous-blue-variable-prescription, --LBV-mass-loss-prescription
 --luminous-blue-variable-multiplier, --main-sequence-core-mass-prescription, --mass-loss-prescription, --overall-wind-mass-loss-multiplier, --wolf-rayet-multiplier, 
---expel-convective-envelope-above-luminosity-threshold, --luminosity-to-mass-threshold, --scale--mass-loss-with-surface-helium-abundance
+--expel-convective-envelope-above-luminosity-threshold, --luminosity-to-mass-threshold, --scale-CHE-mass-loss-with-surface-helium-abundance
 --OB-mass-loss, --OB-mass-loss-prescription, --RSG-mass-loss, --RSG-mass-loss-prescription, --VMS-mass-loss, --vms-mass-loss-prescription, --WR-mass-loss, --WR-mass-loss-prescription
 
 --chemically-homogeneous-evolution, --chemically-homogeneous-evolution-mode

diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp
@@ -2025,6 +2025,10 @@ void BaseBinaryStar::CalculateWindsMassLoss(double p_Dt) {
 
             m_aMassLossDiff = aWinds - m_SemiMajorAxisPrev;                                                                     // change to orbit (semi-major axis) due to winds mass loss
         }
+        else {                                                                                                                  // reset total mass loss rate to zero
+            m_Star1->UpdateTotalMassLossRate(0.0);
+            m_Star2->UpdateTotalMassLossRate(0.0);
+        }
     }
 }
 

diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp
@@ -2451,11 +2451,6 @@ double BaseStar::CalculateMassLossRateBelczynski2010() {
 
             otherWindsRate = CalculateMassLossRateOBVink2001();
             m_DominantMassLossRate = MASS_LOSS_TYPE::OB;
-
-            // If user wants to transition between OB and WR mass loss rates
-            if (OPTIONS->ScaleMassLossWithSurfaceHeliumAbundance()) {
-                otherWindsRate     = EnhanceWindsWithWolfRayetContribution(otherWindsRate, BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0), false);
-            }
         }
 
         if (utils::Compare(LBVRate, otherWindsRate) > 0) {                                                          // which is dominant?
@@ -2506,18 +2501,10 @@ double BaseStar::CalculateMassLossRateMerritt2025() {
         else if (utils::Compare(m_Mass, VMS_MASS_THRESHOLD) >= 0) {                                                 // mass at or above VMS winds threshold?
             otherWindsRate         = CalculateMassLossRateVMS(OPTIONS->VMSMassLossPrescription());                  // yes - use VMS mass loss rate
             m_DominantMassLossRate = MASS_LOSS_TYPE::VMS;                                                           // set dominant mass loss rate
-            // If user wants to transition between OB/VMS and WR mass loss rates
-            if (OPTIONS->ScaleMassLossWithSurfaceHeliumAbundance()) {
-                otherWindsRate     = EnhanceWindsWithWolfRayetContribution(otherWindsRate, 0.0, true);
-            }
         }
         else {                                                                                                      // otherwise...
             otherWindsRate         = CalculateMassLossRateOB(OPTIONS->OBMassLossPrescription());                    // use OB mass loss rate
             m_DominantMassLossRate = MASS_LOSS_TYPE::OB;                                                            // set dominant mass loss rate
-            // If user wants to transition between OB and WR mass loss rates
-            if (OPTIONS->ScaleMassLossWithSurfaceHeliumAbundance()) {
-                otherWindsRate     = EnhanceWindsWithWolfRayetContribution(otherWindsRate, 0.0, true);
-            }
         }
 
         if (utils::Compare(LBVRate, otherWindsRate) > 0) {                                                          // which is dominant?
@@ -2529,69 +2516,6 @@ double BaseStar::CalculateMassLossRateMerritt2025() {
 }
 
 
-/*
- * EnhanceWindsWithWolfRayetContribution
- *
- * Enhance the winds with a contribution due to WR finds (see CalculateMassLossFractionWR)
- *
- * double EnhanceWindsWithWolfRayetContribution (double p_OtherWindsRate, double p_WolfRayetRate, bool p_RecalculateWolfRayetRate)
- *
- * @param       p_OtherWindsRate                Wind mass loss rate due to OB or VMS winds to avoid recompution
- * @param       p_WolfRayetRate                 Wind mass loss rate due to WR winds if already computed
- * @param       p_RecalculateWolfRayetRate      If true, recompute the WR winds by cloning the star as a HeMS star
- * @return                                      Total wind mass loss rate
- */
-double BaseStar::EnhanceWindsWithWolfRayetContribution (double p_OtherWindsRate, double p_WolfRayetRate, bool p_RecalculateWolfRayetRate) {
-
-    double MdotWR = p_WolfRayetRate;
-    double fractionWR = CalculateMassLossFractionWR(m_HeliumAbundanceSurface);
-    double totalWindsRate = 0.0;
-
-    if (p_RecalculateWolfRayetRate && fractionWR > 0.0 ) {
-        MdotWR        = BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0);
-        // *Jeff* This used to clone the star as a HeMS star and query its CalculateMassLossRateMerritt2025(); eventually, let's switch to a static function to calculate the luminosity of the WR star
-    }
-
-    // Combine each of these prescriptions according to the OB wind fraction
-    totalWindsRate = ((1.0 - fractionWR) * p_OtherWindsRate) + (fractionWR * MdotWR);
-
-    if ( fractionWR * MdotWR > (1.0 - fractionWR) * p_OtherWindsRate) {
-        m_DominantMassLossRate =MASS_LOSS_TYPE::WR;
-    }
-
-    return totalWindsRate;
-}
-
-/*
- * CalculateMassLossFractionWR
- *
- * @brief
- * Calculate the fraction of mass loss attributable to WR 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.
- *
- *
- * double CalculateMassLossFractionWR(const double p_HeAbundanceSurface) const
- *
- * @param       p_HeAbundanceSurface            Helium abundance at the surface of the star
- * @return                                      Fraction of mass loss attributable to WR mass loss
- */
-double BaseStar::CalculateMassLossFractionWR(const double p_HeAbundanceSurface) const {
-
-    constexpr double limOB = 0.55;                                          // per Yoon et al. 2006
-    constexpr double limWR = 0.70;                                          // per Yoon et al. 2006
-
-    return std::min(1.0, std::max (0.0, (p_HeAbundanceSurface - limOB) / (limWR - limOB)));
-}
-
 /*
  * Calculate mass loss rate
  *

diff --git a/src/BaseStar.h b/src/BaseStar.h
@@ -615,8 +615,6 @@ class BaseStar {
             double              CalculateMassLossRateWolfRayetTemperatureCorrectionSander2023(const double p_Mdot) const;
             double              CalculateMassLossRateHeliumStarVink2017() const;
     virtual double              CalculateMassLossRateWolfRayetShenar2019() const;
-
-            double              CalculateMassLossFractionWR(const double p_HeAbundanceSurface) const;
 
     virtual double              CalculateMassTransferRejuvenationFactor()                                               { return 1.0; }
 
@@ -692,7 +690,6 @@ class BaseStar {
                                                     const double p_Rand,
                                                     const double p_EjectaMass,
                                                     const double p_RemnantMass);
-            double              EnhanceWindsWithWolfRayetContribution (double p_OtherWindsRate, double p_WolfRayetRate, bool p_RecalculateWolfRayetRate);
 
     virtual void                EvolveOneTimestepPreamble() { };                                                                                                                                    // Default is NO-OP
 

diff --git a/src/CH.cpp b/src/CH.cpp
@@ -285,6 +285,134 @@ void CH::UpdateAgeAfterMassLoss() {
 }
 
 
+/*
+ * CalculateMassLossFractionOB
+ *
+ * @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.
+ *
+ *
+ * double CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const
+ *
+ * @param       p_HeAbundanceSurface            Helium abundance at the surface of the star
+ * @return                                      Fraction of mass loss attributable to OB mass loss
+ */
+double CH::CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const {
+
+    constexpr double limOB = 0.55;                                          // per Yoon et al. 2006
+    constexpr double limWR = 0.70;                                          // per Yoon et al. 2006
+
+    return std::min(1.0, std::max (0.0, (limWR - p_HeAbundanceSurface) / (limWR - limOB)));
+}
+
+
+/*
+ * Calculate the dominant mass loss mechanism and associated rate for the star 
+ * at the current evolutionary phase
+ * 
+ * According to Belczynski et al. 2010 prescription - based on implementation in StarTrack
+ *
+ * Modifications for CH stars
+ * 
+ * double CalculateMassLossRateBelczynski2010()
+ *
+ * @return                                      Mass loss rate in Msol per year
+ */
+double CH::CalculateMassLossRateBelczynski2010() {
+
+    // Define variables
+    double Mdot   = 0.0;
+    double MdotOB = 0.0;
+    double MdotWR = 0.0;
+    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();
+
+    // If user wants to transition between OB and WR mass loss rates
+    if (OPTIONS->ScaleCHEMassLossWithSurfaceHeliumAbundance()) {
+
+        // Calculate WR mass loss rate
+        MdotWR = BaseStar::CalculateMassLossRateWolfRayetZDependent(0.0);
+
+        // Calculate fraction for combining these into total mass-loss rate
+        fractionOB = CalculateMassLossFractionOB(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 = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
+
+    // Enhance mass loss rate due to rotation
+    Mdot *= CalculateMassLossRateEnhancementRotation();
+
+    return Mdot;
+}
+
+
+/*
+ * Calculate the dominant mass loss mechanism and associated rate for the star 
+ * at the current evolutionary phase
+ * 
+ * According to Merritt et al. 2024 prescription
+ *
+ * Modifications for CH stars
+ * 
+ * double CalculateMassLossRateMerritt2025()
+ * 
+ * @return                                      Mass loss rate in Msol per year
+ */
+double CH::CalculateMassLossRateMerritt2025() {
+
+    // Define variables
+    double Mdot   = 0.0;
+    double MdotOB = 0.0;
+    double MdotWR = 0.0;
+    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 chosen prescription
+    MdotOB = BaseStar::CalculateMassLossRateOB(OPTIONS->OBMassLossPrescription());  
+
+    // If user wants to transition between OB and WR mass loss rates
+    if (OPTIONS->ScaleCHEMassLossWithSurfaceHeliumAbundance()) {
+
+        // Here we are going to pretend that this CH star is an HeMS star by
+        // cloning it, so that we can ask it what its mass loss rate would be if it were
+        // a HeMS star
+        HeMS *clone = HeMS::Clone((HeMS&)static_cast<const CH&>(*this), OBJECT_PERSISTENCE::EPHEMERAL, false);  // Do not initialise so that we can use same mass, luminosity, radius etc
+        MdotWR      = clone->CalculateMassLossRateMerritt2025();                                                // Calculate WR mass loss rate              
+        delete clone; clone = nullptr;                                                                          // return the memory allocated for the clone  
+
+        // Calculate weight for combining these into total mass-loss rate
+        fractionOB = CalculateMassLossFractionOB(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 = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
+
+    // Enhance mass loss rate due to rotation
+    Mdot *= CalculateMassLossRateEnhancementRotation();
+
+    return Mdot;
+}
+
+
 STELLAR_TYPE CH::EvolveToNextPhase() {
 
     STELLAR_TYPE stellarType = STELLAR_TYPE::MS_GT_07;

diff --git a/src/CH.h b/src/CH.h
@@ -8,6 +8,7 @@
 
 #include "BaseStar.h"
 #include "MS_gt_07.h"
+#include "HeMS.h"
 
 class BaseStar;
 class MS_gt_07;
@@ -68,11 +69,12 @@ class CH: virtual public BaseStar, public MS_gt_07 {
     double          CalculateLuminosityAtPhaseEnd() const               { return CalculateLuminosityAtPhaseEnd(m_Mass0); }                                                          // Use class member variables
 
     double          CalculateLuminosityOnPhase(const double p_Time, const double p_Mass, const double p_LZAMS) const;
-    double          CalculateLuminosityOnPhase() const                  { return m_Luminosity; }    
+    double          CalculateLuminosityOnPhase() const                  { return CalculateLuminosityOnPhase(m_Age, m_Mass0, m_LZAMS0); }    
 
     // Mass loss rate
-    double          CalculateMassLossRateBelczynski2010()               { return BaseStar::CalculateMassLossRateBelczynski2010() * CalculateMassLossRateEnhancementRotation(); }
-    double          CalculateMassLossRateMerritt2025()                  { return BaseStar::CalculateMassLossRateBelczynski2010() * CalculateMassLossRateEnhancementRotation(); }
+    double          CalculateMassLossRateBelczynski2010();
+    double          CalculateMassLossRateMerritt2025();
+    double          CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const;
 
     // Radius
     double          CalculateRadiusOnPhase() const                      { return m_RZAMS; }                                                                                         // Constant from birth

diff --git a/src/LogTypedefs.h b/src/LogTypedefs.h
@@ -1003,7 +1003,7 @@ enum class PROGRAM_OPTION: int {
     ROTATIONAL_FREQUENCY_1,
     ROTATIONAL_FREQUENCY_2,
 
-    SCALE_MASS_LOSS_SURF_HE_ABUNDANCE,
+    SCALE_CHE_MASS_LOSS_SURF_HE_ABUNDANCE,
     SCALE_TERMINAL_WIND_VEL_METALLICITY_POWER,
     SEMI_MAJOR_AXIS,
     SEMI_MAJOR_AXIS_DISTRIBUTION,
@@ -1233,7 +1233,7 @@ const COMPASUnorderedMap<PROGRAM_OPTION, std::string> PROGRAM_OPTION_LABEL = {
     { PROGRAM_OPTION::ROTATIONAL_FREQUENCY_1,                           "ROTATIONAL_FREQUENCY_1" },
     { PROGRAM_OPTION::ROTATIONAL_FREQUENCY_2,                           "ROTATIONAL_FREQUENCY_2" },
 
-    { PROGRAM_OPTION::SCALE_MASS_LOSS_SURF_HE_ABUNDANCE,                "SCALE_MASS_LOSS_SURF_HE_ABUNDANCE" },
+    { PROGRAM_OPTION::SCALE_CHE_MASS_LOSS_SURF_HE_ABUNDANCE,            "SCALE_CHE_MASS_LOSS_SURF_HE_ABUNDANCE" },
     { PROGRAM_OPTION::SCALE_TERMINAL_WIND_VEL_METALLICITY_POWER,        "SCALE_TERMINAL_WIND_VEL_METALLICITY_POWER" },
     { PROGRAM_OPTION::SEMI_MAJOR_AXIS,                                  "SEMI_MAJOR_AXIS" },
     { PROGRAM_OPTION::SEMI_MAJOR_AXIS_DISTRIBUTION,                     "SEMI_MAJOR_AXIS_DISTRIBUTION" },
@@ -1826,7 +1826,7 @@ const std::map<PROGRAM_OPTION, PROPERTY_DETAILS> PROGRAM_OPTION_DETAIL = {
     { PROGRAM_OPTION::ROTATIONAL_FREQUENCY_1,                                   { TYPENAME::DOUBLE,     "PO_Rotational_Frequency(1)",                "Hz",        24, 15}},
     { PROGRAM_OPTION::ROTATIONAL_FREQUENCY_2,                                   { TYPENAME::DOUBLE,     "PO_Rotational_Frequency(2)",                "Hz",        24, 15}},
 
-    { PROGRAM_OPTION::SCALE_MASS_LOSS_SURF_HE_ABUNDANCE,                        { TYPENAME::BOOL,       "PO_Scale_Mass_Loss_Surf_He_Abundance",  "flag",       0,  0}},
+    { PROGRAM_OPTION::SCALE_CHE_MASS_LOSS_SURF_HE_ABUNDANCE,                    { TYPENAME::BOOL,       "PO_Scale_CHE_Mass_Loss_Surf_He_Abundance",  "flag",       0,  0}},
     { PROGRAM_OPTION::SCALE_TERMINAL_WIND_VEL_METALLICITY_POWER,                { TYPENAME::DOUBLE,     "PO_Scale_Terminal_Wind_Vel_Metallicity_Power", "-",      24, 15}},
     { PROGRAM_OPTION::SEMI_MAJOR_AXIS,                                          { TYPENAME::DOUBLE,     "PO_Semi-Major_Axis",                        "AU",        24, 15}},
     { PROGRAM_OPTION::SEMI_MAJOR_AXIS_DISTRIBUTION,                             { TYPENAME::INT,        "PO_Semi-Major_Axis_Dstrbtn",                "-",          4, 1 }},