Fixed random draws of SN kicks to avoid artificial pile-up at boundaries of distribution

Ilya Mandel · Ilya Mandel · commit f62596ee6c6c · 2025-07-17T11:02:23.000-06:00
diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml
@@ -1,5 +1,9 @@
 ##~!!~## COMPAS option values
+<<<<<<< HEAD
 ##~!!~## File Created Sun Jul  6 16:42:54 2025 by COMPAS v03.21.00
+=======
+##~!!~## File Created Thu Jul 17 09:38:18 2025 by COMPAS v03.21.00
+>>>>>>> dev
 ##~!!~##
 ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has
 ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the
@@ -244,7 +248,7 @@ stringChoices:
 #    --initial-mass-function: 'KROUPA'                                     # Default: 'KROUPA'                      # Options: ['KROUPA','UNIFORM','POWERLAW','SALPETER']
 #    --LBV-mass-loss-prescription: 'HURLEY_ADD'                            # Default: 'HURLEY_ADD'                  # Options: ['BELCZYNSKI','HURLEY','HURLEY_ADD','ZERO']
 #    --main-sequence-core-mass-prescription: 'MANDEL'                      # Default: 'MANDEL'                      # Options: ['BRCEK','MANDEL','ZERO']
-#    --mass-loss-prescription: 'MERRITT2024'                               # Default: 'MERRITT2024'                 # Options: ['MERRITT2024','BELCZYNSKI2010','HURLEY','ZERO']
+#    --mass-loss-prescription: 'MERRITT2025'                               # Default: 'MERRITT2025'                 # Options: ['MERRITT2025','BELCZYNSKI2010','HURLEY','ZERO']
 #    --OB-mass-loss-prescription: 'VINK2021'                               # Default: 'VINK2021'                    # Options: ['KRTICKA2018','BJORKLUND2022','VINK2021','VINK2001','ZERO']
 #    --RSG-mass-loss-prescription: 'DECIN2023'                             # Default: 'DECIN2023'                   # Options: ['NJ90','KEE2021','YANG2023','DECIN2023','BEASOR2020','VINKSABHAHIT2023','ZERO']
 #    --VMS-mass-loss-prescription: 'SABHAHIT2023'                          # Default: 'SABHAHIT2023'                # Options: ['SABHAHIT2023','BESTENLEHNER2020','VINK2011','ZERO']
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
@@ -826,12 +826,12 @@ Default = 0.0
 
 **--mass-loss-prescription** |br|
 Mass loss prescription. |br|
-Options: { ZERO, HURLEY, BELCZYNSKI2010, MERRITT2024 } |br|
+Options: { ZERO, HURLEY, BELCZYNSKI2010, MERRITT2025 } |br|
 ``ZERO``          : No wind mass loss |br|
 ``HURLEY``        : Mass loss according to Hurley (2000) |br|
 ``BELCZYNSKI2010``: Mass loss as per Belczynski (2010), and the default prior to 2023 |br|
-``MERRITT2024``   : Flexible mass loss with phase specific options: (OB, RSG, WR, VMS) |br|
-Default = MERRITT2024 |br|
+``MERRITT2025``   : Flexible mass loss with phase specific options: (OB, RSG, WR, VMS) |br|
+Default = MERRITT2025 |br|
 
 **--mass-ratio [ -q ]** |br|
 Mass ratio :math:`\frac{m2}{m1}` used to determine secondary mass if not specified via ``--initial-mass-2``. |br|
diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst
@@ -3,10 +3,17 @@ What's new
 
 Following is a brief list of important updates to the COMPAS code.  A complete record of changes can be found in the file ``changelog.h``.
 
+<<<<<<< HEAD
 **03.21.00 July 6, 2025**
 * Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
 * Added options to set beta and gamma prescription for second stage of 2-stage CE (``--common-envelope-second-stage-beta``, ``--common-envelope-second-stage-gamma-prescription``)
 * Fixed a bug in the calculation of zeta_equilibrium, which impacts when mass transfer is declared to proceed on a nuclear timescale (and hence how conservative it is)
+=======
+**03.20=1.00 July 17, 2025**
+
+* Deprecated mass loss prescription MERRITT2024 in favour of MERRITT2025
+* Added version strings for gsl, boost, and HDF5 to COMPAS splashscreen
+>>>>>>> dev
 
 **03.20.06 June 25, 2025**
 
diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp
@@ -2469,11 +2469,11 @@ double BaseStar::CalculateMassLossRateBelczynski2010() {
  * or are added to other wind mass loss if LBV_MASS_LOSS_PRESCRIPTION::HURLEY_ADD is used.
  * 
  *
- * double CalculateMassLossRateMerritt2024()
+ * double CalculateMassLossRateMerritt2025()
  * 
  * @return                  Mass loss rate in Msol per year
  */
-double BaseStar::CalculateMassLossRateMerritt2024() {
+double BaseStar::CalculateMassLossRateMerritt2025() {
 
     m_DominantMassLossRate = MASS_LOSS_TYPE::NONE;
 
@@ -2550,8 +2550,8 @@ double BaseStar::CalculateMassLossRate() {
                 mDot = CalculateMassLossRateBelczynski2010();
                 break;
 
-            case MASS_LOSS_PRESCRIPTION::MERRITT2024:
-                mDot = CalculateMassLossRateMerritt2024();
+            case MASS_LOSS_PRESCRIPTION::MERRITT2025:
+                mDot = CalculateMassLossRateMerritt2025();
                 break;
 
             default:                                                                                                // unknown prescription
@@ -3824,9 +3824,8 @@ double BaseStar::DrawRemnantKickMuller(const double p_COCoreMass) const {
 double BaseStar::DrawRemnantKickMullerMandel(const double p_COCoreMass, 
                                              const double p_Rand,
                                              const double p_RemnantMass) const {					
-	double remnantKick = -1.0;
+	double remnantKick;
 	double muKick      = 0.0;
-    double rand        = p_Rand;    // makes it possible to adjust if p_Rand is too low, to avoid getting stuck
 
 	if (utils::Compare(p_RemnantMass, OPTIONS->MaximumNeutronStarMass()) <  0) {
 		muKick = max(OPTIONS->MullerMandelKickMultiplierNS() * (p_COCoreMass - p_RemnantMass) / p_RemnantMass, 0.0);
@@ -3835,10 +3834,9 @@ double BaseStar::DrawRemnantKickMullerMandel(const double p_COCoreMass,
 		muKick = max(OPTIONS->MullerMandelKickMultiplierBH() * (p_COCoreMass - p_RemnantMass) / p_RemnantMass, 0.0);
 	}
 
-	while (remnantKick < 0.0) {
-		remnantKick = muKick * (1.0 + gsl_cdf_gaussian_Pinv(rand, OPTIONS->MullerMandelSigmaKick()));
-		rand        = min(rand + p_Rand + 0.0001, 1.0);
-	}
+    double quantile0 = gsl_cdf_gaussian_P(0.0, sigma);  //quantile of 0 in the Gaussian CDF; the goal is to draw from the cut-off Gaussian since the kick must exceed 0
+    double rand = quantile0 + p_Rand * (1.0 - quantile0);
+    remnantKick = muKick * (1.0 + gsl_cdf_gaussian_Pinv(rand, OPTIONS->MullerMandelSigmaKick()));
 
 	return remnantKick;
 }
diff --git a/src/BaseStar.h b/src/BaseStar.h
@@ -608,7 +608,7 @@ class BaseStar {
             double              CalculateMassLossRateVMSSabhahit2023();
             double              CalculateMassLossRateVMSVink2011() const;
     virtual double              CalculateMassLossRateBelczynski2010();
-    virtual double              CalculateMassLossRateMerritt2024();
+    virtual double              CalculateMassLossRateMerritt2025();
             double              CalculateMassLossRateWolfRayetZDependent(const double p_Mu) const;
             double              CalculateMassLossRateWolfRayet(const double p_Mu) const;
             double              CalculateMassLossRateWolfRayetSanderVink2020(const double p_Mu) const;
diff --git a/src/CH.cpp b/src/CH.cpp
@@ -374,11 +374,11 @@ double CH::CalculateMassLossRateBelczynski2010() {
  *
  * Modifications for CH stars
  * 
- * double CalculateMassLossRateMerritt2024()
+ * double CalculateMassLossRateMerritt2025()
  * 
  * @return                                      Mass loss rate in Msol per year
  */
-double CH::CalculateMassLossRateMerritt2024() {
+double CH::CalculateMassLossRateMerritt2025() {
     
     // Define variables
     double Mdot   = 0.0;
@@ -396,7 +396,7 @@ double CH::CalculateMassLossRateMerritt2024() {
         // 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->CalculateMassLossRateMerritt2024();                                                // Calculate WR mass loss rate              
+        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
diff --git a/src/CH.h b/src/CH.h
@@ -73,7 +73,7 @@ class CH: virtual public BaseStar, public MS_gt_07 {
 
     // Mass loss rate
     double          CalculateMassLossRateBelczynski2010();
-    double          CalculateMassLossRateMerritt2024();
+    double          CalculateMassLossRateMerritt2025();
     double          CalculateMassLossRateWeightOB(const double p_HeliumAbundanceSurface);
     
     // Radius
diff --git a/src/HeMS.cpp b/src/HeMS.cpp
@@ -399,11 +399,11 @@ double HeMS::CalculateMassLossRateWolfRayetShenar2019() const {
  * Calculate the mass loss rate for helium stars in the updated prescription
  * Uses Sander & Vink 2020 for Wolf--Rayet stars
  * 
- * double CalculateMassLossRateMerritt2024()
+ * double CalculateMassLossRateMerritt2025()
  *
  * @return                                      Mass loss rate in Msol per year
  */
-double HeMS::CalculateMassLossRateMerritt2024() {
+double HeMS::CalculateMassLossRateMerritt2025() {
 
     double MdotWR = 0.0;
 
diff --git a/src/HeMS.h b/src/HeMS.h
@@ -47,7 +47,7 @@ class HeMS: virtual public BaseStar, public TPAGB {
     static double   CalculateLuminosityAtPhaseEnd_Static(const double p_Mass);
 
            double   CalculateMassLossRateBelczynski2010();
-           double   CalculateMassLossRateMerritt2024();
+           double   CalculateMassLossRateMerritt2025();
 
     static DBL_DBL  CalculateRadiusAtPhaseEnd_Static(const double p_Mass, const double p_Luminosity);
     static double   CalculateRadiusAtZAMS_Static(const double p_Mass);
diff --git a/src/Options.cpp b/src/Options.cpp
@@ -424,7 +424,7 @@ void Options::OptionValues::Initialise() {
     m_ExpelConvectiveEnvelopeAboveLuminosityThreshold               = false;
     m_LuminosityToMassThreshold                                     = 4.2;      // Podsiadlowski, private communication
 
-    m_MassLossPrescription.type                                     = MASS_LOSS_PRESCRIPTION::MERRITT2024;
+    m_MassLossPrescription.type                                     = MASS_LOSS_PRESCRIPTION::MERRITT2025;
     m_MassLossPrescription.typeString                               = MASS_LOSS_PRESCRIPTION_LABEL.at(m_MassLossPrescription.type);
 
     m_LBVMassLossPrescription.type                                  = LBV_MASS_LOSS_PRESCRIPTION::HURLEY_ADD;
diff --git a/src/Options.h b/src/Options.h
@@ -192,6 +192,9 @@ class Options {
     //       be "false" in the vector, and will be set true if and when the deprecation notice for that
     //       option is shown the first time in a COMPAS run (a deprecation notice for a deprecated option
     //       is only shown once per COMPAS run).
+    //     - datestring indicating the date the option string was deprecated.  Deprecated option strings
+    //       should be manually removed from the "deprecatedOptionStrings" vector 12 months (too long?)
+    //       after the deprecation date.  Datestring format is yyyymmdd (e.g.20251107 indicates November 07, 2025).
     // 
     //
     // "deprecatedOptionValues" vector
@@ -216,23 +219,27 @@ class Options {
     //       be "false" in the vector, and will be set true if and when the deprecation notice for that option
     //       value is shown the first time in a COMPAS run (a deprecation notice for a deprecated option value
     //       is only shown once per COMPAS run).
+    //     - datestring indicating the date the option string was deprecated.  Deprecated option values should
+    //       be manually removed from the "deprecatedOptionValues" vector 12 months (too long?) after the
+    //       deprecation date.  Datestring format is yyyymmdd (e.g.20251107 indicates November 07, 2025).
 
-    std::vector<std::tuple<std::string, std::string, bool>> deprecatedOptionStrings = {
-        { "retain-core-mass-during-caseA-mass-transfer", "", false }
+    std::vector<std::tuple<std::string, std::string, bool, std::string>> deprecatedOptionStrings = {
+        { "retain-core-mass-during-caseA-mass-transfer", "", false, "20250116" }
     };
 
-    std::vector<std::tuple<std::string, std::string, std::string, bool>> deprecatedOptionValues = {
-        { "critical-mass-ratio-prescription",          "GE20",      "GE",        false },
-        { "critical-mass-ratio-prescription",          "GE20_IC",   "GE_IC",     false },
-        { "pulsational-pair-instability-prescription", "COMPAS",    "WOOSLEY",   false },
-	    { "pulsar-birth-spin-period-distribution",     "ZERO",      "NOSPIN",    false },
-        { "tides-prescription",                        "KAPIL2024", "KAPIL2025", false }
+    std::vector<std::tuple<std::string, std::string, std::string, bool, std::string>> deprecatedOptionValues = {
+        { "critical-mass-ratio-prescription",          "GE20",        "GE",          false, "20241118" },
+        { "critical-mass-ratio-prescription",          "GE20_IC",     "GE_IC",       false, "20241118" },
+        { "pulsational-pair-instability-prescription", "COMPAS",      "WOOSLEY",     false, "20250208" },
+	    { "pulsar-birth-spin-period-distribution",     "ZERO",        "NOSPIN",      false, "20250303" },
+        { "tides-prescription",                        "KAPIL2024",   "KAPIL2025",   false, "20250525" },
+        { "mass-loss-prescription",                    "MERRITT2024", "MERRITT2025", false, "20250717" }
     };
 
     // the following vector is used to replace deprecated options in the logfile-definitions file
-    std::vector<std::tuple<std::string, std::string, bool>> deprecatedOptionProperties = {
-        { "black_hole_kicks", "black_hole_kicks_mode",      false },
-        { "lbv_prescription", "LBV-mass-loss-prescription", false }
+    std::vector<std::tuple<std::string, std::string, bool, std::string>> deprecatedOptionProperties = {
+        { "black_hole_kicks", "black_hole_kicks_mode",      false, "20241030" },
+        { "lbv_prescription", "LBV_mass_loss_prescription", false, "20241030" }
     };
 
 
diff --git a/src/Remnants.h b/src/Remnants.h
@@ -74,7 +74,7 @@ class Remnants: virtual public BaseStar, public HeGB {
 
     double          CalculateMassLossRateHurley()                                                               { m_DominantMassLossRate = MASS_LOSS_TYPE::NONE ; return 0.0; }
     double          CalculateMassLossRateBelczynski2010()                                                       { m_DominantMassLossRate = MASS_LOSS_TYPE::NONE ; return 0.0; }
-    double          CalculateMassLossRateMerritt2024()                                                          { m_DominantMassLossRate = MASS_LOSS_TYPE::NONE ; return 0.0; }                                                         // 
+    double          CalculateMassLossRateMerritt2025()                                                          { m_DominantMassLossRate = MASS_LOSS_TYPE::NONE ; return 0.0; }                                                         // 
 
     double          CalculatePerturbationMuOnPhase() const                                                      { return m_Mu; }                                                        // NO-OP
 
diff --git a/src/changelog.h b/src/changelog.h
@@ -1599,14 +1599,28 @@
 //                                          - Fixed error in MainSequence::CalculateInitialMainSequenceCoreMass()
 //  03.20.09  RTW - Jun 30, 2025        - Enhancement:
 //                                          - Added individual velocity components for stars to the LogTypedefs file so they can be included in the output (as ANY_STAR_PROPERTY::VELOCITY_X, or Y, Z)
-//  03.21.00   IM - July 6, 2025        - Enhancements, defect repair:
+//  03.21.00   JR - Jul 17, 2025        - Enhancement:
+//                                          - Changed mass loss prescription MERRITT2024 to MERRIT2025; deprecated MERRIT2024
+//                                          - Added datestrings to vectors supporting deprecation in Options.h to allow timely removal
+//                                          - Added version strings for gsl, boost, and HDF5 to COMPAS splashscreen (for now, gsl & hdf5 are installed versions, boost is version compiled with COMPAS)
+//  03.22.00   IM - July 17, 2025        - Enhancements, defect repair:
 //                                          - Changed default values of --enhance-CHE-lifetimes-luminosities and --scale-CHE-mass-loss-with-surface-helium-abundance to true
 //                                          - Added options to set beta and gamma prescription for second stage of 2-stage CE (--common-envelope-second-stage-beta, --common-envelope-second-stage-gamma-prescription)
 //                                          - Fixed a bug in CalculateZetaEquilibrium(), which impacted when mass transfer is declared nuclear (and how conservative it is)
 //                                          - Added missing virtual declaration to ShouldEnvelopeBeExpelledByPulsations
 //                                          - Now calculate mass accretion rate for nuclear timescale mass transfer on the fly to match with donor mass loss rate set by donor mass loss (required to fit into Roche lobe) divided by time step
+//                                          - Fixed random draws of SN kicks to avoid artificial pile-up at boundaries of distribution
 
 
+// Version string format is MM.mm.rr, where
+//
+// MM is the MAJOR release number: this should be incremented whenever major new functionality is introduced
+// mm is the MINOR release number: this should be incremented whenever minor new functionality (e.g. small enhancemet) is introduced
+// rr is the fix RELEASE number:   this should be incremented whenever a defect repair is made
+//
+// if MM is incremented, set mm and rr to 00, even if defect repairs and minor enhancements were also made
+// if mm is incremented, set rr to 00, even if defect repairs were also made
+
 const std::string VERSION_STRING = "03.21.00";
 
 
diff --git a/src/typedefs.h b/src/typedefs.h
@@ -598,12 +598,12 @@ enum class MASS_CUTOFF: int {
 };
 
 // mass loss prescriptions
-enum class MASS_LOSS_PRESCRIPTION: int { ZERO, HURLEY, BELCZYNSKI2010, MERRITT2024 };
+enum class MASS_LOSS_PRESCRIPTION: int { ZERO, HURLEY, BELCZYNSKI2010, MERRITT2025 };
 const COMPASUnorderedMap<MASS_LOSS_PRESCRIPTION, std::string> MASS_LOSS_PRESCRIPTION_LABEL = {
     { MASS_LOSS_PRESCRIPTION::ZERO,           "ZERO" },
     { MASS_LOSS_PRESCRIPTION::HURLEY,         "HURLEY" },
     { MASS_LOSS_PRESCRIPTION::BELCZYNSKI2010, "BELCZYNSKI2010" },
-    { MASS_LOSS_PRESCRIPTION::MERRITT2024,    "MERRITT2024" }
+    { MASS_LOSS_PRESCRIPTION::MERRITT2025,    "MERRITT2025" }
 };
 
 // symbolic names for mass loss rate type
diff --git a/src/utils.cpp b/src/utils.cpp
diff --git a/src/utils.h b/src/utils.h

Original file line number	Diff line number	Diff line change
`@@ -399,11 +399,11 @@ double HeMS::CalculateMassLossRateWolfRayetShenar2019() const {`
`399`	`399`	`* Calculate the mass loss rate for helium stars in the updated prescription`
`400`	`400`	`* Uses Sander & Vink 2020 for Wolf--Rayet stars`
`401`	`401`	`*`
`402`		`- * double CalculateMassLossRateMerritt2024()`
	`402`	`+ * double CalculateMassLossRateMerritt2025()`
`403`	`403`	`*`
`404`	`404`	`* @return Mass loss rate in Msol per year`
`405`	`405`	`*/`
`406`		`-double HeMS::CalculateMassLossRateMerritt2024() {`
	`406`	`+double HeMS::CalculateMassLossRateMerritt2025() {`
`407`	`407`
`408`	`408`	`double MdotWR = 0.0;`
`409`	`409`