From 9b1f62e53330adf1dbe7e10001034b32f9b031e2 Mon Sep 17 00:00:00 2001 From: Reinhold Willcox Date: Wed, 18 Dec 2024 10:05:52 +0100 Subject: [PATCH] Fixed supernova Euler angles code. Incomplete logic was leading to a div by zero. --- src/BaseBinaryStar.cpp | 30 +++++++++++++++++------------- 1 file changed, 17 insertions(+), 13 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 609e0bf5..93f4a033 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1312,19 +1312,23 @@ void BaseBinaryStar::ResolveSupernova() { // then the cross product is not well-defined, and we need to account for degeneracy between eccentricity vectors. // Also, if either eccentricity is 0.0, then the eccentricity vector is not well defined. - if ((utils::Compare(m_ThetaE, 0.0) == 0) && // orbitalAngularMomentumVectorPrev parallel to orbitalAngularMomentumVector? - ((utils::Compare(eccentricityPrev, 0.0) > 0) && (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ... and both eccentricityVectorPrev and eccentricityVector well-defined? - - double psiPlusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi + phi is constant - m_PhiE = _2_PI * RAND->Random(); - m_PsiE = psiPlusPhi - m_PhiE; - } - else if ((utils::Compare(m_ThetaE, M_PI) == 0) && // orbitalAngularMomentumVectorPrev anti-parallel to orbitalAngularMomentumVector? - ((utils::Compare(eccentricityPrev, 0.0) > 0) && (utils::Compare(m_Eccentricity, 0.0) > 0))) { // ... and both eccentricityVectorPrev and eccentricityVector well-defined? - - double psiMinusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi - phi is constant - m_PhiE = _2_PI * RAND->Random(); - m_PsiE = psiMinusPhi + m_PhiE; + if ((utils::Compare(m_ThetaE, 0.0) == 0) || (utils::Compare(m_ThetaE, M_PI) == 0)) { // orbitalAngularMomentumVectorPrev parallel or anti-parallel to orbitalAngularMomentumVector + if ((utils::Compare(eccentricityPrev, 0.0) == 0) || (utils::Compare(m_Eccentricity, 0.0) == 0)) { // either e_prev or e_now is 0, so eccentricity vector is not well-defined + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = _2_PI * RAND->Random(); + } + else { // both eccentricityVectorPrev and eccentricityVector well-defined + if (utils::Compare(m_ThetaE, 0.0) == 0){ // Orbital AM is parallel ? + double psiPlusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // yes - then psi + phi is constant + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = psiPlusPhi - m_PhiE; + } + else { + double psiMinusPhi = angleBetween(eccentricityVector, eccentricityVectorPrev); // no - then psi - phi is constant + m_PhiE = _2_PI * RAND->Random(); + m_PsiE = psiMinusPhi + m_PhiE; + } + } } else { // neither - the cross product of the orbit normals is well-defined Vector3d orbitalPivotAxis = cross(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // cross product of the orbit normals