M*a should be constant on wind mass loss

TeamCOMPAS · Dec 5, 2024 · 3f33a58 · 3f33a58
1 parent 1f2e382
commit 3f33a58
Showing 1 changed file with 4 additions and 1 deletion.
diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp
@@ -1925,8 +1925,11 @@ void BaseBinaryStar::CalculateWindsMassLoss() {
             double mWinds1 = m_Star1->CalculateMassLossValues(true);                                                            // calculate new values assuming mass loss applied
             double mWinds2 = m_Star2->CalculateMassLossValues(true);                                                            // calculate new values assuming mass loss applied
 
-            double aWinds  = m_SemiMajorAxisPrev / (2.0 - ((m_Star1->MassPrev() + m_Star2->MassPrev()) / (mWinds1 + mWinds2))); // new semi-major axis for circularlised orbit
+            double aWinds  = m_SemiMajorAxisPrev * (m_Star1->Mass() + m_Star2->Mass()) / (mWinds1 + mWinds2);                   // new semi-major axis after wind mass loss, integrated to ensure a*M conservation
 
+            std::cout<<"Types"<<(int)m_Star1->StellarType()<<" "<<(int)m_Star2->StellarType()<<"prev masses"<<m_Star1->Mass()<<" "<<m_Star2->Mass()<<std::endl;
+            std::cout<<" new masses"<<mWinds1<<" "<<mWinds2<<"old a "<<m_SemiMajorAxisPrev<<"new a "<<aWinds<<std::endl;
+
             m_Star1->SetMassLossDiff(mWinds1 - m_Star1->Mass());                                                                // JR: todo: find a better way?
             m_Star2->SetMassLossDiff(mWinds2 - m_Star2->Mass());                                                                // JR: todo: find a better way?