Merge pull request #1062 from jeffriley/issue-1057

Fix for issue 1057

src/EAGB.cpp

@@ -772,6 +772,12 @@ double EAGB::CalculateRadiusOnPhase_Static(const double      p_Mass,
                                            const DBL_VECTOR &p_BnCoefficients) {
 #define b p_BnCoefficients  // for convenience and readability - undefined at end of function
 
+    // sanity check for mass and luminosity - just return 0.0 if mass or luminosity <= 0
+    // doing this will save some compute cycles - but it is not strictly what the equation in Hurley says
+    // (Hurley et al. 2000, eq 74 and immediately following) - there if mass is 0.0 but luminosity is
+    // non-zero we get a non-zero value for radius (shouldn't happen, but we need to code for all possibilities)
+    if (utils::Compare(p_Mass, 0.0) <= 0 || utils::Compare(p_Luminosity, 0.0) <= 0) return 0.0;
+
     // calculate radius constant A (Hurley et al. 2000, eq 74)
     // and coefficient b(50)
     double A;

src/HeGB.cpp

@@ -41,6 +41,9 @@ double HeGB::CalculateLuminosityOnPhase_Static(const double p_CoreMass, const do
  */
 std::tuple<double, double> HeGB::CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Luminosity) {
 
+    // sanity check for mass and luminosity - just return 0.0 if mass or luminosity <= 0
+    if (utils::Compare(p_Mass, 0.0) <= 0 || utils::Compare(p_Luminosity, 0.0) <= 0) return std::make_tuple(0.0, 0.0);
+
     double RZHe = HeMS::CalculateRadiusAtZAMS_Static(p_Mass);
     double LTHe = CalculateLuminosityAtPhaseEnd_Static(p_Mass);
 

src/HeHG.h

@@ -77,6 +77,8 @@ class HeHG: virtual public BaseStar, public HeMS {
 
             double          CalculateMassTransferRejuvenationFactor() const;
 
+            double          CalculateMomentOfInertia() const                                                        { return GiantBranch::CalculateMomentOfInertia(); }
+
             double          CalculatePerturbationMu() const;
             double          CalculatePerturbationMuAtPhaseEnd() const                                               { return m_Mu; }                                                        // NO-OP
 

src/HeMS.cpp

@@ -144,6 +144,9 @@ double HeMS::CalculateRadiusAtZAMS_Static(const double p_Mass) {
  */
 double HeMS::CalculateRadiusOnPhase_Static(const double p_Mass, const double p_Tau) {
 
+    // sanity check for mass - just return 0.0 if mass <= 0
+    if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0;
+
     double tau_6 = p_Tau * p_Tau * p_Tau * p_Tau * p_Tau * p_Tau;   // pow() is slow - use multiplication
     double beta  = std::max(0.0, 0.4 - 0.22 * log10(p_Mass));
 

src/HeMS.h

@@ -84,6 +84,8 @@ class HeMS: virtual public BaseStar, public TPAGB {
 
             double          CalculateMassTransferRejuvenationFactor() const;
 
+            double          CalculateMomentOfInertia() const                                                        { return MainSequence::CalculateMomentOfInertia(); }
+
             double          CalculatePerturbationMu() const                                                         { return 5.0; }                                                                 // Hurley et al. 2000, eqs 97 & 98
 
             double          CalculateRadialExtentConvectiveEnvelope() const                                         { return BaseStar::CalculateRadialExtentConvectiveEnvelope(); }                 // HeMS stars don't have a convective envelope

src/WhiteDwarfs.cpp

@@ -151,8 +151,13 @@ double WhiteDwarfs::CalculateLuminosityOnPhase_Static(const double p_Mass, const
  * @return                                      Radius of a White Dwarf in Rsol (since WD is ~ Earth sized, expect answer around 0.009)
  */
 double WhiteDwarfs::CalculateRadiusOnPhase_Static(const double p_Mass) {
+
+    // sanity check for mass - just return 0.0 if mass <= 0
+    if (utils::Compare(p_Mass, 0.0) <= 0) return 0.0;
+
     double MCH_Mass_one_third  = std::cbrt(MCH / p_Mass); 
     double MCH_Mass_two_thirds = MCH_Mass_one_third * MCH_Mass_one_third;
+
     return std::max(NEUTRON_STAR_RADIUS, 0.0115 * std::sqrt((MCH_Mass_two_thirds - 1.0 / MCH_Mass_two_thirds)));
 }
 

src/changelog.h

@@ -1090,8 +1090,13 @@
 //                                      - Fix for issue #1048
 // 02.41.05     YS - Jan 31, 2024    - Bug fix:
 //                                      - Fix for issue #1058: fixing calculation of pulsar spin period
+// 02.41.06     JR - Feb 10, 2024    - Defect repair:
+//                                      - Fix for issue #1057:
+//                                            HeMS::CalculateMomentOfInertia() falls back to MainSequence::CalculateMomentOfInertia()
+//                                            HeHG::CalculateMomentOfInertia() falls back to GiantBranch::CalculateMomentOfInertia()
+//                                      - Added sanity checks for mass and luminosity where necessary in variants of CalculateRadiusOnPhase_Static()
 
-const std::string VERSION_STRING = "02.41.05";
+const std::string VERSION_STRING = "02.41.06";
 
 
 # endif // __changelog_h__