Adding fix for issue #1257

TeamCOMPAS · Jan 7, 2025 · a879403 · a879403
1 parent c938c8f
commit a879403
Show file tree

Hide file tree

Showing 2 changed files with 32 additions and 7 deletions.
diff --git a/src/NS.cpp b/src/NS.cpp
@@ -330,15 +330,28 @@ void NS::CalculateAndSetPulsarParameters() {
 
     m_PulsarDetails.magneticField     = PPOW(10.0, CalculateBirthMagneticField()) * GAUSS_TO_TESLA;         // magnetic field in Gauss -> convert to Tesla
     m_PulsarDetails.spinPeriod        = CalculateBirthSpinPeriod();                                         // spin period in ms
-    m_PulsarDetails.spinFrequency     = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS);
-    m_PulsarDetails.birthPeriod       = m_PulsarDetails.spinPeriod * SECONDS_IN_MS;                         // convert from ms to s 
-
     m_MomentOfInertia_CGS             = CalculateMomentOfInertiaCGS();                                      // in CGS g cm^2
 
-    // Note we convert neutronStarMomentOfInertia from CGS to SI here
-    m_PulsarDetails.spinDownRate      = CalculateSpinDownRate(m_PulsarDetails.spinFrequency, m_MomentOfInertia_CGS, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM);  
-    m_PulsarDetails.birthSpinDownRate = m_PulsarDetails.spinDownRate; 
-    m_AngularMomentum_CGS             = m_MomentOfInertia_CGS * m_PulsarDetails.spinFrequency;              // in CGS g cm^2 s^-1
+    if ((utils::Compare(m_PulsarDetails.spinPeriod, 0.0) == 0) || 
+        (utils::Compare(m_PulsarDetails.magneticField, 0.0) == 0)) {
+            // if initial spin period or magnetic field is set to ZERO
+            // then assume spin or magnetic field evolution is not calculated
+            // and use 0.0 as placeholders instead of going through calculations with 0s 
+            m_PulsarDetails.spinDownRate = 0.0;
+            m_PulsarDetails.spinFrequency = 0.0;
+            m_PulsarDetails.spinPeriod = 0.0; 
+            m_PulsarDetails.magneticField = 0.0;
+            m_AngularMomentum_CGS = 0.0;
+    }
+    else {
+        m_PulsarDetails.spinFrequency     = _2_PI / (m_PulsarDetails.spinPeriod * SECONDS_IN_MS);
+        m_PulsarDetails.birthPeriod       = m_PulsarDetails.spinPeriod * SECONDS_IN_MS;                         // convert from ms to s 
+
+        // Note we convert neutronStarMomentOfInertia from CGS to SI here
+        m_PulsarDetails.spinDownRate      = CalculateSpinDownRate(m_PulsarDetails.spinFrequency, m_MomentOfInertia_CGS, m_PulsarDetails.magneticField, m_Radius * RSOL_TO_KM);  
+        m_PulsarDetails.birthSpinDownRate = m_PulsarDetails.spinDownRate; 
+        m_AngularMomentum_CGS             = m_MomentOfInertia_CGS * m_PulsarDetails.spinFrequency;             // in CGS g cm^2 s^-1
+    }
 }
 
 
@@ -496,6 +509,17 @@ DBL_DBL_DBL NS::DeltaAngularMomentumByPulsarAccretion_Static(const double p_Mass
  */
 void NS::UpdateMagneticFieldAndSpin(const bool p_CommonEnvelope, const bool p_RecycledNS, const double p_Stepsize, const double p_MassGainPerTimeStep, const double p_Epsilon) {
 
+    if ((utils::Compare(m_PulsarDetails.spinPeriod, 0.0) == 0) || 
+        (utils::Compare(m_PulsarDetails.magneticField, 0.0) == 0)) {
+            // if spin period or magnetic field is zero, set all pulsar parameters to 0 
+            // instead of doing calculations with 0s. 
+            m_PulsarDetails.spinDownRate = 0.0;
+            m_PulsarDetails.spinFrequency = 0.0;
+            m_PulsarDetails.spinPeriod = 0.0; 
+            m_PulsarDetails.magneticField = 0.0;
+        return; 
+    }
+
     double magFieldLowerLimit = PPOW(10.0, OPTIONS->PulsarLog10MinimumMagneticField()) * GAUSS_TO_TESLA;    
     double magFieldLowerLimit_G = magFieldLowerLimit * TESLA_TO_GAUSS ;
     double kappa              = OPTIONS->PulsarMagneticFieldDecayMassscale() * MSOL_TO_KG;     

diff --git a/src/changelog.h b/src/changelog.h
@@ -1429,6 +1429,7 @@
 //                                      - Update to NS::UpdateMagneticFieldAndSpin(). For NS accreting in RLOF, adding BOOST ODE integrator to calculate the spin-up/recycling. 
 //                                      - Added program option "--neutron-star-accretion-in-ce" to account for how a neutron star would accrete mass during a common envelope event. 
 //                                      - Fix for issue #1002
+//                                      - Fixing for issue #1257
 
 const std::string VERSION_STRING = "03.11.00";