Partial commit

Author: Ilya Mandel

src/BaseStar.cpp

@@ -2448,8 +2448,14 @@ double BaseStar::CalculateMassLossRateBelczynski2010() {
             otherWindsRate = CalculateMassLossRateHurley() * OPTIONS->CoolWindMassLossMultiplier();                 // apply cool wind mass loss multiplier
         }
         else  {                                                                                                     // hot stars, add Vink et al. 2001 winds (ignoring bistability jump)
+            
             otherWindsRate = CalculateMassLossRateOBVink2001();
-            m_DominantMassLossRate = MASS_LOSS_TYPE::OB;                                                            // set dominant mass loss rate
+            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?
@@ -2499,11 +2505,19 @@ 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                             
+            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?
@@ -2515,6 +2529,72 @@ 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 have the code below instead; I'd prefer a static function to calculate the luminosity of the WR star.
+        //HeMS *clone = HeMS::Clone((HeMS&)static_cast<const MainSequence&>(*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
+    }
+    
+    // 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
  *
@@ -3849,7 +3929,8 @@ double BaseStar::DrawRemnantKickMullerMandel(const double p_COCoreMass,
     double rand = quantile0 + p_Rand * (1.0 - quantile0);
     remnantKick = muKick * (1.0 + gsl_cdf_gaussian_Pinv(rand, sigmaKick));
 
-    if (utils::SNEventType(m_SupernovaDetails.events.current) == SN_EVENT::USSN ) {                  // overwrite USSN kick prescription
+    // Mandel * Mueller 2020 call for USSN kicks to be treated in the same way as CCSN kicks; however, if this override flag is set, set the USSN kick to be equal to the user-provided magnitude
+    if (utils::SNEventType(m_SupernovaDetails.events.current) == SN_EVENT::USSN && OPTIONS->USSNKicksOverrideMandelMuller() ) {
         remnantKick = OPTIONS->KickMagnitudeDistributionSigmaForUSSN();
     }
 

src/BaseStar.h

@@ -615,6 +615,8 @@ 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; }
 
@@ -691,7 +693,8 @@ 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
 
             STELLAR_TYPE        EvolveOnPhase(const double p_DeltaTime);

src/CH.cpp

@@ -285,134 +285,6 @@ 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;

src/CH.h

@@ -8,7 +8,6 @@
 
 #include "BaseStar.h"
 #include "MS_gt_07.h"
-#include "HeMS.h"
 
 class BaseStar;
 class MS_gt_07;
@@ -72,9 +71,8 @@ class CH: virtual public BaseStar, public MS_gt_07 {
     double          CalculateLuminosityOnPhase() const                  { return m_Luminosity; }    
 
     // Mass loss rate
-    double          CalculateMassLossRateBelczynski2010();
-    double          CalculateMassLossRateMerritt2025();
-    double          CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const;
+    double          CalculateMassLossRateBelczynski2010()               { return BaseStar::CalculateMassLossRateBelczynski2010() * CalculateMassLossRateEnhancementRotation(); }
+    double          CalculateMassLossRateMerritt2025()                  { return BaseStar::CalculateMassLossRateBelczynski2010() * CalculateMassLossRateEnhancementRotation(); }
 
     // Radius
     double          CalculateRadiusOnPhase() const                      { return m_RZAMS; }                                                                                         // Constant from birth

src/LogTypedefs.h

@@ -1003,7 +1003,7 @@ enum class PROGRAM_OPTION: int {
     ROTATIONAL_FREQUENCY_1,
     ROTATIONAL_FREQUENCY_2,
 
-    SCALE_CHE_MASS_LOSS_SURF_HE_ABUNDANCE,
+    SCALE_MASS_LOSS_SURF_HE_ABUNDANCE,
     SCALE_TERMINAL_WIND_VEL_METALLICITY_POWER,
     SEMI_MAJOR_AXIS,
     SEMI_MAJOR_AXIS_DISTRIBUTION,
@@ -1013,6 +1013,8 @@ enum class PROGRAM_OPTION: int {
     STELLAR_ZETA_PRESCRIPTION,
 
     TIDES_PRESCRIPTION,
+    
+    USSN_KICKS_OVERRIDE_MANDEL_MULLER,
 
     WR_FACTOR,
 
@@ -1231,8 +1233,8 @@ 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_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::SCALE_MASS_LOSS_SURF_HE_ABUNDANCE,                "SCALE_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" },
     { PROGRAM_OPTION::SEMI_MAJOR_AXIS_DISTRIBUTION_MAX,                 "SEMI_MAJOR_AXIS_DISTRIBUTION_MAX" },
@@ -1241,6 +1243,8 @@ const COMPASUnorderedMap<PROGRAM_OPTION, std::string> PROGRAM_OPTION_LABEL = {
     { PROGRAM_OPTION::STELLAR_ZETA_PRESCRIPTION,                        "STELLAR_ZETA_PRESCRIPTION" },
 
     { PROGRAM_OPTION::TIDES_PRESCRIPTION,                               "TIDES_PRESCRIPTION" },
+    
+    { PROGRAM_OPTION::USSN_KICKS_OVERRIDE_MANDEL_MULLER,                "USSN_KICKS_OVERRIDE_MANDEL_MULLER" },
 
     { PROGRAM_OPTION::WR_FACTOR,                                        "WR_FACTOR" },
 
@@ -1822,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_CHE_MASS_LOSS_SURF_HE_ABUNDANCE,                    { TYPENAME::BOOL,       "PO_Scale_CHE_Mass_Loss_Surf_He_Abundance",  "flag",       0,  0}},
+    { PROGRAM_OPTION::SCALE_MASS_LOSS_SURF_HE_ABUNDANCE,                        { TYPENAME::BOOL,       "PO_Scale_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 }},
@@ -1832,6 +1836,8 @@ const std::map<PROGRAM_OPTION, PROPERTY_DETAILS> PROGRAM_OPTION_DETAIL = {
     { PROGRAM_OPTION::STELLAR_ZETA_PRESCRIPTION,                                { TYPENAME::INT,        "PO_Stellar_Zeta_Prscrptn",                  "-",          4, 1 }},
 
     { PROGRAM_OPTION::TIDES_PRESCRIPTION,                                       { TYPENAME::INT,        "PO_Tides_Prscrptn",                         "-",          4, 1 }},
+    
+    { PROGRAM_OPTION::USSN_KICKS_OVERRIDE_MANDEL_MULLER,                        { TYPENAME::BOOL,       "PO_USSN_Kicks_Override_Mandel_Muller",  "flag",       0,  0}},
 
     { PROGRAM_OPTION::WR_FACTOR,                                                { TYPENAME::DOUBLE,     "PO_WR_Factor",                              "-",         24, 15}},