Skip to content

Fixed bug #1247, div by zero in SN Euler angles #1314

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
ilyamandel merged 2 commits into from
Dec 18, 2024
Merged

Fixed bug #1247, div by zero in SN Euler angles #1314

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
9b1f62e
Fixed supernova Euler angles code. Incomplete logic was leading to a …
reinhold-willcox Dec 18, 2024
fcdb390
changelog
reinhold-willcox Dec 18, 2024
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
30 changes: 17 additions & 13 deletions src/BaseBinaryStar.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
4 changes: 3 additions & 1 deletion src/changelog.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1421,7 +1421,9 @@
// - exactly preserve the product of semi-major axis * total mass on wind mass loss
// 03.10.03 JR - Dec 16, 2024 - Defect repair:
// - fix for issue #1310 - run terminates prematurely if error in grid file
// 03.10.04 RTW - Nov 27, 2024 - Defect repair:
// - fix for issue #1247 - SN Euler angles had incomplete logic, leading to a div by zero in some cases

const std::string VERSION_STRING = "03.10.03";
const std::string VERSION_STRING = "03.10.04";

# endif // __changelog_h__
Loading