diff --git a/compas_python_utils/preprocessing/compasConfigDefault.yaml b/compas_python_utils/preprocessing/compasConfigDefault.yaml index a7ba0e29d..fd582524c 100644 --- a/compas_python_utils/preprocessing/compasConfigDefault.yaml +++ b/compas_python_utils/preprocessing/compasConfigDefault.yaml @@ -1,5 +1,5 @@ ##~!!~## COMPAS option values -##~!!~## File Created Thu Oct 23 15:46:49 2025 by COMPAS v03.26.01 +##~!!~## File Created Mon Oct 27 09:51:18 2025 by COMPAS v03.26.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-CHE-mass-loss-with-surface-helium-abundance: True # Default: True +# --scale-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 @@ -57,6 +57,7 @@ booleanChoices: # --allow-non-stripped-ECSN: False # Default: False # --pair-instability-supernovae: True # Default: True # --pulsational-pair-instability: True # Default: True +# --USSN-kicks-override-Mandel-Muller: False # Default: False ### PULSAR PARAMETERS # --evolve-pulsars: False # Default: False 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 index 97ddb81a8..0397aebb0 100644 --- 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 ` -**--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. +**--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. Default = TRUE **--scale-terminal-wind-velocity-with-metallicity-power** |br| @@ -1431,6 +1431,10 @@ This option is primarily intended for debugging/testing of convergence issues ra Enable mass transfer. |br| Default = TRUE +**--USSN-kicks-override-mandel-muller** |br| +Use user-defined USSN kicks (as a fixed value) in lieu of the Mandel & Muller kick prescription for USSNe. |br| +Default = FALSE + .. _options-props-V: :ref:`Back to Top ` @@ -1531,7 +1535,7 @@ Go to :ref:`the top of this page ` 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, +--expel-convective-envelope-above-luminosity-threshold, --luminosity-to-mass-threshold, --scale--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 @@ -1583,8 +1587,8 @@ Go to :ref:`the top of this page ` for the full alphabetical --kick-magnitude-distribution, --kick-magnitude-sigma-CCSN-BH, --kick-magnitude-sigma-CCSN-NS, --kick-magnitude-sigma-ECSN, --kick-magnitude-sigma-USSN, --black-hole-kicks, --black-hole-kicks-mode, --fix-dimensionless-kick-magnitude, --kick-magnitude, --kick-magnitude-1, --kick-magnitude-2, --kick-magnitude-min, --kick-magnitude-max, --kick-magnitude-random, --kick-magnitude-random-1, --kick-magnitude-random-2, --kick-scaling-factor, -muller-mandel-kick-multiplier-BH, ---muller-mandel-kick-multiplier-NS, --muller-mandel-sigma-kick-BH, --muller-mandel-sigma-kick-NS, --kick-direction, ---kick-direction-distribution, --kick-direction-power, --kick-mean-anomaly-1, --kick-mean-anomaly-2, --kick-phi-1, --kick-phi-2, --kick-theta-1, --kick-theta-2 +--muller-mandel-kick-multiplier-NS, --muller-mandel-sigma-kick-BH, --muller-mandel-sigma-kick-NS, --USSN-kicks-override-mandel-muller, +--kick-direction, --kick-direction-distribution, --kick-direction-power, --kick-mean-anomaly-1, --kick-mean-anomaly-2, --kick-phi-1, --kick-phi-2, --kick-theta-1, --kick-theta-2 :ref:`Back to Top ` diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 221a28393..5669dc356 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,11 @@ 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``. +**03.26.02 October 27, 2025** + +* Added option --USSN-kicks-override-mandel-muller ; if set to true, use user-defined USSN kicks (as a fixed value) in lieu of the Mandel & Muller kick prescription for USSNe +* Replaced --scale-CHE-mass-loss-with-surface-helium-abundance with the more general --scale-mass-loss-with-surface-helium-abundance (applies to all MS stars, not just CHE stars) + **03.26.00 September 2, 2025** * Added HAMSTARS mass transfer efficiency prescription diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index 4fc20e175..dc8de97a7 100755 --- a/src/BaseStar.cpp +++ b/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,69 @@ 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 * @@ -2967,28 +3044,14 @@ double BaseStar::CalculateOStarRotationalVelocityAnalyticCDF_Static(const double boost::math::inverse_gamma_distribution<> gammaComponent(alpha, beta); // (shape, scale) = (alpha, beta) boost::math::normal_distribution<> normalComponent(mu, sigma); + + // Compute CDF at zero rotational velocity -- the CDF should be relative to this quantity + double CDFzero = (iGamma * boost::math::cdf(gammaComponent, 0.0)) + ((1.0 - iGamma) * boost::math::cdf(normalComponent, 0.0)); + + double CDFunnormalised = (iGamma * boost::math::cdf(gammaComponent, p_Ve)) + ((1.0 - iGamma) * boost::math::cdf(normalComponent, p_Ve)); + + return ((CDFunnormalised-CDFzero) / (1.0 - CDFzero)); - return (iGamma * boost::math::cdf(gammaComponent, p_Ve)) + ((1.0 - iGamma) * boost::math::cdf(normalComponent, p_Ve)); -} - - -/* - * Calculate the inverse of the analytic cumulative distribution function (CDF) for the - * equatorial rotational velocity of single O stars. - * - * (i.e. calculate the inverse of CalculateOStarRotationalVelocityAnalyticCDF_Static()) - * - * - * double CalculateOStarRotationalVelocityAnalyticCDFInverse_Static(const double p_Ve, const void *p_Params) - * - * @param [IN] p_vE Rotational velocity (in km s^-1) - value of the kick vk which we want to find - * @param [IN] p_Params Pointer to RotationalVelocityParams structure containing y, the CDF draw U(0,1) - * @return Inverse CDF - * Should be zero when p_Ve = vk, the value of the kick to draw - */ -double BaseStar::CalculateOStarRotationalVelocityAnalyticCDFInverse_Static(double p_Ve, void* p_Params) { - RotationalVelocityParams* params = (RotationalVelocityParams*) p_Params; - return CalculateOStarRotationalVelocityAnalyticCDF_Static(p_Ve) - params->u; } @@ -3001,68 +3064,93 @@ double BaseStar::CalculateOStarRotationalVelocityAnalyticCDFInverse_Static(doubl * Ramirez-Agudelo et al. 2013 https://arxiv.org/abs/1309.2929 * * - * double CalculateOStarRotationalVelocity_Static(const double p_Xmin, const double p_Xmax) + * double CalculateOStarRotationalVelocity * - * @param [IN] p_Xmin Minimum value for root - * @param [IN] p_Xmax Maximum value for root * @return Rotational velocity in km s^-1 */ -double BaseStar::CalculateOStarRotationalVelocity_Static(const double p_Xmin, const double p_Xmax) { - - double xMin = p_Xmin; - double xMax = p_Xmax; - - double result = xMin; - - double maximumInverse = CalculateOStarRotationalVelocityAnalyticCDF_Static(xMax); - double minimumInverse = CalculateOStarRotationalVelocityAnalyticCDF_Static(xMin); - - double rand = RAND->Random(); - - while (utils::Compare(rand, maximumInverse) > 0) { - xMax *= 2.0; - maximumInverse = CalculateOStarRotationalVelocityAnalyticCDF_Static(xMax); - } - - if (utils::Compare(rand, minimumInverse) >= 0) { - - const gsl_root_fsolver_type *T; - gsl_root_fsolver *s; - gsl_function F; - - RotationalVelocityParams params = {rand}; - - F.function = &CalculateOStarRotationalVelocityAnalyticCDFInverse_Static; - F.params = ¶ms; - - // gsl_root_fsolver_brent - // gsl_root_fsolver_bisection - T = gsl_root_fsolver_brent; - s = gsl_root_fsolver_alloc(T); - - gsl_root_fsolver_set(s, &F, xMin, xMax); - - int status = GSL_CONTINUE; - int iter = 0; - int maxIter = 100; - - while (status == GSL_CONTINUE && iter < maxIter) { - iter++; - status = gsl_root_fsolver_iterate(s); - result = gsl_root_fsolver_root(s); - xMin = gsl_root_fsolver_x_lower(s); - xMax = gsl_root_fsolver_x_upper(s); - status = gsl_root_test_interval(xMin, xMax, 0, 0.001); +double BaseStar::CalculateOStarRotationalVelocity() { + + double desiredCDF = RAND->Random(); // Random desired CDF + + const boost::uintmax_t maxit = ADAPTIVE_RV_MAX_ITERATIONS; // Limit to maximum iterations. + boost::uintmax_t it = maxit; // Initially our chosen max iterations, but updated with actual. + + // find root + // we use an iterative algorithm to find the root here: + // - if the root finder throws an exception, we stop and return a negative value for the root (indicating no root found) + // - if the root finder reaches the maximum number of (internal) iterations, we stop and return a negative value for the root (indicating no root found) + // - if the root finder returns a solution, we check that func(solution) = 0.0 +/ ROOT_ABS_TOLERANCE + // - if the solution is acceptable, we stop and return the solution + // - if the solution is not acceptable, we reduce the search step size and try again + // - if we reach the maximum number of search step reduction iterations, or the search step factor reduces to 1.0 (so search step size = 0.0), + // we stop and return a negative value for the root (indicating no root found) + + double guess = 100.0; // guess at 100 km s^-1 (arbitrary initial guess) + + double factorFrac = ADAPTIVE_RV_SEARCH_FACTOR_FRAC; // search step size factor fractional part + double factor = 1.0 + factorFrac; // factor to determine search step size (size = guess * factor) + + std::pair root(-1.0, -1.0); // initialise root - default return + std::size_t tries = 0; // number of tries + bool done = false; // finished (found root or exceed maximum tries)? + ERROR error = ERROR::NONE; + OStarRotationVelocityFunctor func = OStarRotationVelocityFunctor(desiredCDF); + while (!done) { // while no error and acceptable root found + + bool isRising = true; //guess for direction of search; CDF increases monotonically + + // run the root finder + // regardless of any exceptions or errors, display any problems as a warning, then + // check if the root returned is within tolerance - so even if the root finder + // bumped up against the maximum iterations, or couldn't bracket the root, use + // whatever value it ended with and check if it's good enough for us - not finding + // an acceptable root should be the exception rather than the rule, so this strategy + // shouldn't cause undue performance issues. + try { + error = ERROR::NONE; + root = boost::math::tools::bracket_and_solve_root(func, guess, factor, isRising, utils::BracketTolerance, it); // find root + // root finder returned without raising an exception + if (error != ERROR::NONE) { SHOW_WARN(error); } // root finder encountered an error + else if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_RV_ITERATIONS); } // too many root finder iterations + } + catch(std::exception& e) { // catch generic boost root finding error + // root finder exception + // could be too many iterations, or unable to bracket root - it may not + // be a hard error - so no matter what the reason is that we are here, + // we'll just emit a warning and keep trying + if (it >= maxit) { SHOW_WARN(ERROR::TOO_MANY_RV_ITERATIONS); } // too many root finder iterations + else { SHOW_WARN(ERROR::ROOT_FINDER_FAILED, e.what()); } // some other problem - show it as a warning } - // JR: should we issue a warning, or throw an error, if the root finder didn't actually find the roor here (i.e. we stopped because pf maxIter)? - // To be consistent, should we use the Boost root solver here? - // **Ilya** both questions above -- IM: yes to both, TBC - - gsl_root_fsolver_free(s); // de-allocate memory for root solver + // we have a solution from the root finder - it may not be an acceptable solution + // so we check if it is within our preferred tolerance + if (fabs(func(root.first + (root.second - root.first) / 2.0)) <= ROOT_ABS_TOLERANCE) { // solution within tolerance? + done = true; // yes - we're done + } + else if (fabs(func(root.first)) <= ROOT_ABS_TOLERANCE) { // solution within tolerance at endpoint 1? + root.second=root.first; + done = true; // yes - we're done + } + else if (fabs(func(root.second)) <= ROOT_ABS_TOLERANCE) { // solution within tolerance at endpoint 2? + root.first=root.second; + done = true; // yes - we're done + } + else { // no - try again + // we don't have an acceptable solution - reduce search step size and try again + factorFrac /= 2.0; // reduce fractional part of factor + factor = 1.0 + factorFrac; // new search step size + tries++; // increment number of tries + if (tries > ADAPTIVE_RV_MAX_TRIES || fabs(factor - 1.0) <= ROOT_ABS_TOLERANCE) { // too many tries, or step size 0.0? + // we've tried as much as we can - fail here with -ve return value + root.first = -1.0; // yes - set error return + root.second = -1.0; + SHOW_WARN(ERROR::TOO_MANY_RV_TRIES); // show warning + done = true; // we're done + } + } } - - return result; + + return root.first + (root.second - root.first) / 2.0; // Midway between brackets is our result, if necessary we could return the result as an interval here. } @@ -3096,16 +3184,18 @@ double BaseStar::CalculateRotationalVelocity(double p_MZAMS) { case ROTATIONAL_VELOCITY_DISTRIBUTION::VLTFLAMES: // VLTFLAMES // Rotational velocity based on VLT-FLAMES survey. - // For O-stars use results of Ramirez-Agudelo et al. (2013) https://arxiv.org/abs/1309.2929 (single stars) + // For O-stars (taken to be above 16 Msol), use results + // of Ramirez-Agudelo et al. (2013) https://arxiv.org/abs/1309.2929 (single stars) // and Ramirez-Agudelo et al. (2015) https://arxiv.org/abs/1507.02286 (spectroscopic binaries) - // For B-stars use results of Dufton et al. (2013) https://arxiv.org/abs/1212.2424 - // For lower mass stars, I don't know what updated results there are so default back to - // Hurley et al. 2000 distribution for now + // For B-stars (taken to be between 2 and 16 Msol) use results + // of Dufton et al. (2013) https://arxiv.org/abs/1212.2424 + // For lower mass stars, default back to Hurley et al. 2000 distribution for now - if (utils::Compare(p_MZAMS, 16.0) >= 0) { // JR: what does 16.0 represent? Not another mass threshold that should be in constants.h ...? /*ilya*/ - vRot = CalculateOStarRotationalVelocity_Static(0.0, 800.0); + if (utils::Compare(p_MZAMS, 16.0) >= 0) { + vRot = CalculateOStarRotationalVelocity(); + vRot = max(vRot, 0.0); // Set to no rotation if no positive solution found; warning already raised } - else if (utils::Compare(p_MZAMS, 2.0) >= 0) { // JR: what does 2.0 represent? Not another mass threshold that should be in constants.h ...? **Ilya** + else if (utils::Compare(p_MZAMS, 2.0) >= 0) { vRot = utils::InverseSampleFromTabulatedCDF(RAND->Random(), BStarRotationalVelocityCDFTable); } else { @@ -3848,7 +3938,12 @@ double BaseStar::DrawRemnantKickMullerMandel(const double p_COCoreMass, double quantile0 = gsl_cdf_gaussian_P(-1.0, sigmaKick); //quantile of -1 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, sigmaKick)); - + + // 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(); + } + return remnantKick; } diff --git a/src/BaseStar.h b/src/BaseStar.h index c1dd4aea7..19f95a37a 100644 --- a/src/BaseStar.h +++ b/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; } @@ -625,8 +627,7 @@ class BaseStar { static double CalculateOpacity_Static(const double p_HeliumAbundanceSurface); static double CalculateOStarRotationalVelocityAnalyticCDF_Static(const double p_Ve); - static double CalculateOStarRotationalVelocityAnalyticCDFInverse_Static(double p_Ve, void *p_Params); - static double CalculateOStarRotationalVelocity_Static(const double p_Xmin, const double p_Xmax); + double CalculateOStarRotationalVelocity(); double CalculatePerturbationB(const double p_Mass) const; double CalculatePerturbationC(double p_Mass) const; @@ -691,7 +692,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); @@ -740,6 +742,34 @@ class BaseStar { void UpdateAttributesAndAgeOneTimestepPreamble(const double p_DeltaMass, const double p_DeltaMass0, const double p_DeltaTime); + + /* + * Functor for CalculateOStarRotationalVelocity() + * + * + * Constructor: initialise the class + * template OStarRotationVelocityFunctor(double p_CDF, ERROR *p_Error) + * + * @param [IN] p_CDF Desired CDF value + * + * Function: calculate the CDF of the O star rotational velocity and compare to desired value + * T OStarRotationVelocityFunctor(double const& p_Ve) + * + * @param [IN] p_Ve Rotational velocity, km s^-1 + * @return Difference between star's Roche Lobe radius and radius after mass loss + */ + template + struct OStarRotationVelocityFunctor { + OStarRotationVelocityFunctor(double p_CDF) { + m_CDF = p_CDF; + } + T operator()(double const& p_Ve) { + + return (CalculateOStarRotationalVelocityAnalyticCDF_Static(p_Ve) - m_CDF); + } + private: + double m_CDF; + }; }; #endif // __BaseStar_h__ diff --git a/src/CH.cpp b/src/CH.cpp index ab2b59b6c..2896b46e2 100755 --- a/src/CH.cpp +++ b/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(*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 index 9a8ff7300..22d83fedb 100755 --- a/src/CH.h +++ b/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 diff --git a/src/ErrorCatalog.h b/src/ErrorCatalog.h index fea6e5ce4..8a28a7c06 100644 --- a/src/ErrorCatalog.h +++ b/src/ErrorCatalog.h @@ -133,6 +133,8 @@ enum class ERROR: int { TOO_MANY_RETRIES, // generic too many retries TOO_MANY_RLOF_ITERATIONS, // too many iterations in RLOF root finder TOO_MANY_RLOF_TRIES, // too many tries in RLOF root finder + TOO_MANY_RV_ITERATIONS, // too many iterations in rotational velocity root finder + TOO_MANY_RV_TRIES, // too many tries in rotational velocity root finder TOO_MANY_TIMESTEPS_IN_TIMESTEPS_FILE, // too many timesteps in timesteps file (exceeds maximum) UNABLE_TO_CREATE_DIRECTORY, // unable to create directory UNABLE_TO_REMOVE_DIRECTORY, // unable to remove directory @@ -310,6 +312,8 @@ const COMPASUnorderedMap> ERROR_CATA { ERROR::TOO_MANY_RETRIES, { ERROR_SCOPE::ALWAYS, "Too many retries" }}, { ERROR::TOO_MANY_RLOF_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of iterations when fitting star inside Roche Lobe in RLOF" }}, { ERROR::TOO_MANY_RLOF_TRIES, { ERROR_SCOPE::ALWAYS, "Reached maximum number of tries when fitting star inside Roche Lobe in RLOF" }}, + { ERROR::TOO_MANY_RV_ITERATIONS, { ERROR_SCOPE::ALWAYS, "Reached maximum number of tries when inverting the CDF of rotational velocitiesF" }}, + { ERROR::TOO_MANY_RV_TRIES, { ERROR_SCOPE::ALWAYS, "Reached maximum number of tries when inverting the CDF of rotational velocities" }}, { ERROR::TOO_MANY_TIMESTEPS_IN_TIMESTEPS_FILE, { ERROR_SCOPE::ALWAYS, "Number of timesteps in timestpes file exceeds maximum timesteps" }}, { ERROR::UNABLE_TO_CREATE_DIRECTORY, { ERROR_SCOPE::ALWAYS, "Unable to create directory" }}, { ERROR::UNABLE_TO_REMOVE_DIRECTORY, { ERROR_SCOPE::ALWAYS, "Unable to remove directory" }}, diff --git a/src/HeHG.cpp b/src/HeHG.cpp index 3a6235eae..de947ab54 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -270,7 +270,7 @@ double HeHG::CalculateRemnantRadius() const { * * This implementation adapted from the STARTRACK implementation (STARTRACK courtesy Chris Belczynski) * - * This function good for HeHG and HeGB stars (for Helium stars: always use Natasha's fit) + * This function is for HeHG and HeGB stars (for Helium stars: always use Natasha Ivanova's fit) * * * double CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) @@ -282,11 +282,11 @@ double HeHG::CalculateRemnantRadius() const { */ double HeHG::CalculateLambdaNanjingStarTrack(const double p_Mass, const double p_Metallicity) const { - double rMin = 0.25; // minimum considered radius: Natasha JR: should this be in constants.h? Maybe not... Who is Natasha? **Ilya** - double rMax = 120.0; // maximum considered radius: Natasha JR: should this be in constants.h? Maybe not... Who is Natasha? **Ilya** + double rMin = 0.25; // minimum considered radius: Natasha + double rMax = 120.0; // maximum considered radius: Natasha - double rMinLambda = 0.3 * PPOW(rMin, -0.8); // JR: todo: should this be in constants.h? JR: should this be in constants.h? Maybe not... **Ilya** - double rMaxLambda = 0.3 * PPOW(rMax, -0.8); // JR: todo: should this be in constants.h? JR: should this be in constants.h? Maybe not... **Ilya** + double rMinLambda = 0.3 * PPOW(rMin, -0.8); + double rMaxLambda = 0.3 * PPOW(rMax, -0.8); return m_Radius < rMin ? rMinLambda : (m_Radius > rMax ? rMaxLambda : 0.3 * PPOW(m_Radius, -0.8)); } diff --git a/src/LogTypedefs.h b/src/LogTypedefs.h index 69ff56539..b3f579b77 100644 --- a/src/LogTypedefs.h +++ b/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_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_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_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_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}}, diff --git a/src/Options.cpp b/src/Options.cpp index bc2d3c2a1..9aa26b105 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -322,6 +322,7 @@ void Options::OptionValues::Initialise() { m_KickMagnitudeDistributionSigmaForECSN = 30.0; m_KickMagnitudeDistributionSigmaForUSSN = 30.0; m_KickScalingFactor = 1.0; + m_USSNKicksOverrideMandelMuller = false; // Default from Mandel&Mueller 2020 is to treat USSNe as normal CCSNe for kicks // Kick direction option m_KickDirectionDistribution.type = KICK_DIRECTION_DISTRIBUTION::ISOTROPIC; @@ -375,7 +376,7 @@ void Options::OptionValues::Initialise() { m_CheMode.type = CHE_MODE::PESSIMISTIC; m_CheMode.typeString = CHE_MODE_LABEL.at(m_CheMode.type); m_EnhanceCHELifetimesLuminosities = true; // default is to enhance - m_ScaleCHEMassLossWithSurfaceHeliumAbundance = true; // default is to scale the mass loss + m_ScaleMassLossWithSurfaceHeliumAbundance = true; // default is to scale the mass loss // Supernova remnant mass prescription options m_RemnantMassPrescription.type = REMNANT_MASS_PRESCRIPTION::MULLERMANDEL; @@ -1018,15 +1019,20 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt ("Print switch log to file (default = " + std::string(p_Options->m_SwitchLog ? "TRUE" : "FALSE") + ")").c_str() ) ( - "scale-CHE-mass-loss-with-surface-helium-abundance", - po::value(&p_Options->m_ScaleCHEMassLossWithSurfaceHeliumAbundance)->default_value(p_Options->m_ScaleCHEMassLossWithSurfaceHeliumAbundance)->implicit_value(true), - ("Whether to transition mass loss rates for chemically homogeneously evolving (CHE) stars between OB mass loss rates and Wolf-Rayet (WR) mass loss rates as a function of the surface helium abundance (Ys) as described by Yoon et al. 2006 (default = " + std::string(p_Options->m_ScaleCHEMassLossWithSurfaceHeliumAbundance ? "TRUE" : "FALSE") + ")").c_str() + "scale-mass-loss-with-surface-helium-abundance", + po::value(&p_Options->m_ScaleMassLossWithSurfaceHeliumAbundance)->default_value(p_Options->m_ScaleMassLossWithSurfaceHeliumAbundance)->implicit_value(true), + ("Whether to transition mass loss rates for stars between OB/VMS mass loss rates and Wolf-Rayet (WR) mass loss rates as a function of the surface helium abundance (Ys) as described by Yoon et al. 2006 (default = " + std::string(p_Options->m_ScaleMassLossWithSurfaceHeliumAbundance ? "TRUE" : "FALSE") + ")").c_str() ) ( "use-mass-transfer", po::value(&p_Options->m_UseMassTransfer)->default_value(p_Options->m_UseMassTransfer)->implicit_value(true), ("Enable mass transfer (default = " + std::string(p_Options->m_UseMassTransfer ? "TRUE" : "FALSE") + ")").c_str() ) + ( + "USSN-kicks-override-Mandel-Muller", + po::value(&p_Options->m_USSNKicksOverrideMandelMuller)->default_value(p_Options->m_USSNKicksOverrideMandelMuller)->implicit_value(true), + ("Whether to use user-defined USSN kicks (as a fixed value) in lieu of the Mandel & Muller kick prescription for USSNe (default = " + std::string(p_Options->m_USSNKicksOverrideMandelMuller ? "TRUE" : "FALSE") + ")").c_str() + ) // numerical options - alphabetically grouped by type @@ -5113,7 +5119,7 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const { case PROGRAM_OPTION::ROTATIONAL_FREQUENCY_1 : value = RotationalFrequency1(); break; case PROGRAM_OPTION::ROTATIONAL_FREQUENCY_2 : value = RotationalFrequency2(); break; - case PROGRAM_OPTION::SCALE_CHE_MASS_LOSS_SURF_HE_ABUNDANCE : value = ScaleCHEMassLossWithSurfaceHeliumAbundance(); break; + case PROGRAM_OPTION::SCALE_MASS_LOSS_SURF_HE_ABUNDANCE : value = ScaleMassLossWithSurfaceHeliumAbundance(); break; case PROGRAM_OPTION::SCALE_TERMINAL_WIND_VEL_METALLICITY_POWER : value = ScaleTerminalWindVelocityWithMetallicityPower(); break; case PROGRAM_OPTION::SEMI_MAJOR_AXIS : value = SemiMajorAxis(); break; @@ -5124,6 +5130,7 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const { case PROGRAM_OPTION::STELLAR_ZETA_PRESCRIPTION : value = static_cast(StellarZetaPrescription()); break; case PROGRAM_OPTION::TIDES_PRESCRIPTION : value = static_cast(TidesPrescription()); break; + case PROGRAM_OPTION::USSN_KICKS_OVERRIDE_MANDEL_MULLER : value = static_cast(USSNKicksOverrideMandelMuller()); break; case PROGRAM_OPTION::WR_FACTOR : value = WolfRayetFactor(); break; diff --git a/src/Options.h b/src/Options.h index b3e93c5e1..4f03d33fa 100755 --- a/src/Options.h +++ b/src/Options.h @@ -710,6 +710,7 @@ class Options { "use-mass-loss", "use-mass-transfer", + "USSN-kicks-override-mandel-muller", "VMW-mass-loss-prescription", "version", "v", @@ -968,6 +969,8 @@ class Options { double m_MullerMandelKickNS; // Multiplier for NS kicks per Mandel and Mueller, 2020 double m_MullerMandelSigmaKickBH; // Scatter for BH kicks per Mandel and Mueller, 2020 double m_MullerMandelSigmaKickNS; // Scatter for NS kicks per Mandel and Mueller, 2020 + bool m_USSNKicksOverrideMandelMuller; // Use user-defined USSN kicks (as a fixed value) in lieu of the Mandel & Muller kick prescription for USSNe + // Black hole kicks ENUM_OPT m_BlackHoleKicksMode; // Which black hole kicks mode @@ -983,7 +986,7 @@ class Options { // CHE - Chemically Homogeneous Evolution ENUM_OPT m_CheMode; // Which Chemically Homogeneous Evolution mode bool m_EnhanceCHELifetimesLuminosities; // Whether to enhance the lifetimes and luminosities of CHE stars relative to SSE MS stars - bool m_ScaleCHEMassLossWithSurfaceHeliumAbundance; // Whether to transition between OB and WR mass loss rates for CHE stars on the MS + bool m_ScaleMassLossWithSurfaceHeliumAbundance; // Whether to transition between OB/VMS and WR mass loss rates for stars on the MS/CHE // Supernova remnant mass ENUM_OPT m_RemnantMassPrescription; // Which remnant mass prescription @@ -1735,7 +1738,7 @@ class Options { double RotationalFrequency2() const { return OPT_VALUE("rotational-frequency-2", m_RotationalFrequency2, true); } RSG_MASS_LOSS_PRESCRIPTION RSGMassLossPrescription() const { return OPT_VALUE("RSG-mass-loss-prescription", m_RSGMassLossPrescription.type, true); } - bool ScaleCHEMassLossWithSurfaceHeliumAbundance() const { return OPT_VALUE("scale-CHE-mass-loss-with-surface-helium-abundance", m_ScaleCHEMassLossWithSurfaceHeliumAbundance, false); } + bool ScaleMassLossWithSurfaceHeliumAbundance() const { return OPT_VALUE("scale-mass-loss-with-surface-helium-abundance", m_ScaleMassLossWithSurfaceHeliumAbundance, false); } double ScaleTerminalWindVelocityWithMetallicityPower() const { return OPT_VALUE("scale-terminal-wind-velocity-with-metallicity-power", m_ScaleTerminalWindVelocityWithMetallicityPower, true);} double SemiMajorAxis() const { return OPT_VALUE("semi-major-axis", m_SemiMajorAxis, true); } SEMI_MAJOR_AXIS_DISTRIBUTION SemiMajorAxisDistribution() const { return OPT_VALUE("semi-major-axis-distribution", m_SemiMajorAxisDistribution.type, true); } @@ -1771,6 +1774,7 @@ class Options { bool UseMassTransfer() const { return OPT_VALUE("use-mass-transfer", m_UseMassTransfer, true); } bool UsePairInstabilitySupernovae() const { return OPT_VALUE("pair-instability-supernovae", m_UsePairInstabilitySupernovae, true); } bool UsePulsationalPairInstability() const { return OPT_VALUE("pulsational-pair-instability", m_UsePulsationalPairInstability, true); } + bool USSNKicksOverrideMandelMuller() const { return OPT_VALUE("USSN-kicks-override-Mandel-Muller", m_USSNKicksOverrideMandelMuller, true); } VMS_MASS_LOSS_PRESCRIPTION VMSMassLossPrescription() const { return OPT_VALUE("VMS-mass-loss-prescription", m_VMSMassLossPrescription.type, true); } double WolfRayetFactor() const { return OPT_VALUE("wolf-rayet-multiplier", m_WolfRayetFactor, true); } diff --git a/src/changelog.h b/src/changelog.h index cd88383f9..0abe179db 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1669,6 +1669,10 @@ // - Corrected errors in combining OB and WR winds in CH::CalculateMassLossRateBelczynski2010(), CalculateMassLossRateMerritt2025() and CH::CalculateMassLossFractionOB() [previously CalculateMassLossRateWeightOB()] // 03.26.01 AB - October 24, 2025 - Option name change: // - Main sequence core mass prescription ZERO renamed to HURLEY; deprecated ZERO +// 03.26.02 IM - October 27, 2025 - Enhancements +// - Added option --USSN-kicks-override-mandel-muller ; if set to true, use user-defined USSN kicks (as a fixed value) in lieu of the Mandel & Muller kick prescription for USSNe +// - Replaced --scale-CHE-mass-loss-with-surface-helium-abundance with the more general --scale-mass-loss-with-surface-helium-abundance (applies to all MS stars, not just CHE stars) +// - Updated rotational velocity solver to use boost root finder // // Version string format is MM.mm.rr, where // @@ -1679,7 +1683,7 @@ // 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.26.01"; +const std::string VERSION_STRING = "03.26.02"; # endif // __changelog_h__ diff --git a/src/constants.h b/src/constants.h index a77b18a1d..090e0f764 100755 --- a/src/constants.h +++ b/src/constants.h @@ -279,6 +279,10 @@ constexpr int ADAPTIVE_RLOF_MAX_TRIES = 30; constexpr int ADAPTIVE_RLOF_MAX_ITERATIONS = 50; // Maximum number of root finder iterations in BaseBinaryStar::MassLossToFitInsideRocheLobe() constexpr double ADAPTIVE_RLOF_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in BaseBinaryStar::MassLossToFitInsideRocheLobe() (added to 1.0) +constexpr int ADAPTIVE_RV_MAX_TRIES = 30; // Maximum number of tries in BaseStar::CalculateOStarRotationalVelocity_Static() +constexpr int ADAPTIVE_RV_MAX_ITERATIONS = 50; // Maximum number of root finder iterations in BaseStar::CalculateOStarRotationalVelocity_Static() +constexpr double ADAPTIVE_RV_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in BaseStar::CalculateOStarRotationalVelocity_Static() (added to 1.0) + constexpr int ADAPTIVE_MASS0_MAX_TRIES = 30; // Maximum number of tries in HG::Mass0ToMatchDesiredCoreMass() constexpr int ADAPTIVE_MASS0_MAX_ITERATIONS = 50; // Maximum number of iterations in HG::Mass0ToMatchDesiredCoreMass() constexpr double ADAPTIVE_MASS0_SEARCH_FACTOR_FRAC = 1.0; // Search size factor (fractional part) in HG::Mass0ToMatchDesiredCoreMass() (added to 1.0) diff --git a/src/yaml.h b/src/yaml.h index 5689b1e8d..82f2df4d4 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -95,7 +95,7 @@ namespace yaml { " --enhance-CHE-lifetimes-luminosities", " --expel-convective-envelope-above-luminosity-threshold", " --natal-kick-for-PPISN", - " --scale-CHE-mass-loss-with-surface-helium-abundance", + " --scale-mass-loss-with-surface-helium-abundance", "", " ### BINARY PROPERTIES", " --allow-touching-at-birth # record binaries that have stars touching at birth in output files", @@ -123,6 +123,7 @@ namespace yaml { " --allow-non-stripped-ECSN", " --pair-instability-supernovae", " --pulsational-pair-instability", + " --USSN-kicks-override-Mandel-Muller", "", " ### PULSAR PARAMETERS", " --evolve-pulsars",