@@ -285,6 +285,134 @@ void CH::UpdateAgeAfterMassLoss() {
285285}
286286
287287
288+ /*
289+ * CalculateMassLossFractionOB
290+ *
291+ * @brief
292+ * Calculate the fraction of mass loss attributable to OB mass loss, per Yoon et al. 2006
293+ *
294+ * The model described in Yoon et al. 2006 (also Szecsi et al. 2015) uses OB mass loss while the
295+ * He surface abundance is below 0.55, WR mass loss when the surface He abundance is above 0.7,
296+ * and linearly interpolate when the He surface abundance is between those limits.
297+ *
298+ * This function calculates the fraction of mass loss attributable to OB mass loss, based on
299+ * the He surface abundance and the abundance limits described in Yoon et al. 2006. The value
300+ * returned will be 1.0 if 100% of the mass loss is attributable to OB mass lass, 0.0 if 100% of
301+ * the mass loss is attributable to WR mass loss, and in the range (0.0, 1.0) if the mass loss is
302+ * a mix of OB and WR.
303+ *
304+ *
305+ * double CalculateMassLossFractionOB(const double p_HeAbundanceSurface) const
306+ *
307+ * @param p_HeAbundanceSurface Helium abundance at the surface of the star
308+ * @return Fraction of mass loss attributable to OB mass loss
309+ */
310+ double CH::CalculateMassLossFractionOB (const double p_HeAbundanceSurface) const {
311+
312+ constexpr double limOB = 0.55 ; // per Yoon et al. 2006
313+ constexpr double limWR = 0.70 ; // per Yoon et al. 2006
314+
315+ return std::min (1.0 , std::max (0.0 , (limWR - p_HeAbundanceSurface) / (limWR - limOB)));
316+ }
317+
318+
319+ /*
320+ * Calculate the dominant mass loss mechanism and associated rate for the star
321+ * at the current evolutionary phase
322+ *
323+ * According to Belczynski et al. 2010 prescription - based on implementation in StarTrack
324+ *
325+ * Modifications for CH stars
326+ *
327+ * double CalculateMassLossRateBelczynski2010()
328+ *
329+ * @return Mass loss rate in Msol per year
330+ */
331+ double CH::CalculateMassLossRateBelczynski2010 () {
332+
333+ // Define variables
334+ double Mdot = 0.0 ;
335+ double MdotOB = 0.0 ;
336+ double MdotWR = 0.0 ;
337+ double fractionOB = 1.0 ; // Initialised to 1.0 to allow us to use the OB mass loss rate by default
338+
339+ // Calculate OB mass loss rate according to the Vink et al. formalism
340+ MdotOB = BaseStar::CalculateMassLossRateBelczynski2010 ();
341+
342+ // If user wants to transition between OB and WR mass loss rates
343+ if (OPTIONS->ScaleCHEMassLossWithSurfaceHeliumAbundance ()) {
344+
345+ // Calculate WR mass loss rate
346+ MdotWR = BaseStar::CalculateMassLossRateWolfRayetZDependent (0.0 );
347+
348+ // Calculate fraction for combining these into total mass-loss rate
349+ fractionOB = CalculateMassLossFractionOB (m_HeliumAbundanceSurface);
350+
351+ }
352+
353+ // Finally, combine each of these prescriptions according to the OB wind fraction
354+ Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR);
355+
356+ // Set dominant mass loss rate
357+ m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
358+
359+ // Enhance mass loss rate due to rotation
360+ Mdot *= CalculateMassLossRateEnhancementRotation ();
361+
362+ return Mdot;
363+ }
364+
365+
366+ /*
367+ * Calculate the dominant mass loss mechanism and associated rate for the star
368+ * at the current evolutionary phase
369+ *
370+ * According to Merritt et al. 2024 prescription
371+ *
372+ * Modifications for CH stars
373+ *
374+ * double CalculateMassLossRateMerritt2025()
375+ *
376+ * @return Mass loss rate in Msol per year
377+ */
378+ double CH::CalculateMassLossRateMerritt2025 () {
379+
380+ // Define variables
381+ double Mdot = 0.0 ;
382+ double MdotOB = 0.0 ;
383+ double MdotWR = 0.0 ;
384+ double fractionOB = 1.0 ; // Initialised to 1.0 to allow us to use the OB mass loss rate by default
385+
386+ // Calculate OB mass loss rate according to the chosen prescription
387+ MdotOB = BaseStar::CalculateMassLossRateOB (OPTIONS->OBMassLossPrescription ());
388+
389+ // If user wants to transition between OB and WR mass loss rates
390+ if (OPTIONS->ScaleCHEMassLossWithSurfaceHeliumAbundance ()) {
391+
392+ // Here we are going to pretend that this CH star is an HeMS star by
393+ // cloning it, so that we can ask it what its mass loss rate would be if it were
394+ // a HeMS star
395+ 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
396+ MdotWR = clone->CalculateMassLossRateMerritt2025 (); // Calculate WR mass loss rate
397+ delete clone; clone = nullptr ; // return the memory allocated for the clone
398+
399+ // Calculate weight for combining these into total mass-loss rate
400+ fractionOB = CalculateMassLossFractionOB (m_HeliumAbundanceSurface);
401+ }
402+
403+ // Finally, combine each of these prescriptions according to the OB wind fraction
404+ Mdot = (fractionOB * MdotOB) + ((1.0 - fractionOB) * MdotWR);
405+
406+ // Set dominant mass loss rate
407+ m_DominantMassLossRate = (fractionOB * MdotOB) > ((1.0 - fractionOB) * MdotWR) ? MASS_LOSS_TYPE::OB : MASS_LOSS_TYPE::WR;
408+
409+ // Enhance mass loss rate due to rotation
410+ Mdot *= CalculateMassLossRateEnhancementRotation ();
411+
412+ return Mdot;
413+ }
414+
415+
288416STELLAR_TYPE CH::EvolveToNextPhase () {
289417
290418 STELLAR_TYPE stellarType = STELLAR_TYPE::MS_GT_07;
0 commit comments