From 22a0b25d0456b57b4626d570438d74542396b95b Mon Sep 17 00:00:00 2001 From: Reinhold Willcox Date: Wed, 27 Nov 2024 16:59:36 +0100 Subject: [PATCH 1/3] fixed issues with vector3d indexing and rotation --- src/BaseBinaryStar.cpp | 32 ++++++++++++++++++++++++++++++-- src/typedefs.h | 7 +++++++ src/vector3d.cpp | 9 +++++---- src/vector3d.h | 4 ++-- 4 files changed, 44 insertions(+), 8 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 4ef7a6487..ff95f28fe 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1158,6 +1158,8 @@ void BaseBinaryStar::ResolveSupernova() { #define mag Magnitude() #define hat UnitVector() + Vector3d thisVect; + // set relevant preSN parameters m_EccentricityPreSN = m_Eccentricity; m_SemiMajorAxisPreSN = m_SemiMajorAxis; @@ -1170,6 +1172,7 @@ void BaseBinaryStar::ResolveSupernova() { double theta = m_Supernova->SN_Theta(); // angle out of the binary plane double phi = m_Supernova->SN_Phi(); // angle in the binary plane Vector3d natalKickVector = m_Supernova->SN_KickMagnitude() * Vector3d(cos(theta) * cos(phi), cos(theta) * sin(phi), sin(theta)); + thisVect = natalKickVector;std::cout<<"kick vec : "<CalculateSNAnomalies(eccentricityPrev); double cosEccAnomaly = cos(m_Supernova->SN_EccentricAnomaly()); double sinEccAnomaly = sin(m_Supernova->SN_EccentricAnomaly()); + if ((utils::Compare(eccentricityPrev, 0.0) == 0) && m_Companion->IsOneOf(SN_REMNANTS)) { // If circular and first SN, fix eccentric anomaly to 0 + cosEccAnomaly = 1; + sinEccAnomaly = 0; + } // Derived quantities double aPrev = semiMajorAxisPrev_km; @@ -1210,13 +1217,16 @@ void BaseBinaryStar::ResolveSupernova() { double omega = std::sqrt(G_km_Msol_s * totalMassPrev / aPrev_3); // Keplerian orbital frequency (rad/s) Vector3d separationVectorPrev = Vector3d(aPrev * (cosEccAnomaly - eccentricityPrev), aPrev * (sinEccAnomaly) * sqrt1MinusEccPrevSquared, 0.0); // relative position vector, from m1Prev to m2Prev (km) + thisVect = separationVectorPrev;std::cout<<"sep initial vec : "<UpdateComponentVelocity(component1VelocityVectorAtInfinity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); - m_Companion->UpdateComponentVelocity(component2VelocityVectorAtInfinity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); + m_Companion->UpdateComponentVelocity(tmpvec); // Set Euler Angles m_ThetaE = angleBetween(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // angle between the angular momentum unit vectors, always well defined diff --git a/src/typedefs.h b/src/typedefs.h index 23dc6e6cd..a60d38dbc 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -246,6 +246,13 @@ const STELLAR_TYPE_LIST COMPACT_OBJECTS = { STELLAR_TYPE::MASSLESS_REMNANT }; +// (convenience) initializer list for SN REMNANTS +const STELLAR_TYPE_LIST SN_REMNANTS = { + STELLAR_TYPE::NEUTRON_STAR, + STELLAR_TYPE::BLACK_HOLE, + STELLAR_TYPE::MASSLESS_REMNANT +}; + // (convenience) initializer list for GIANTS const STELLAR_TYPE_LIST GIANTS = { STELLAR_TYPE::FIRST_GIANT_BRANCH, diff --git a/src/vector3d.cpp b/src/vector3d.cpp index d7dadd57f..86539b2b5 100644 --- a/src/vector3d.cpp +++ b/src/vector3d.cpp @@ -81,8 +81,6 @@ Vector3d Vector3d::ChangeBasis(const double p_ThetaE, const double p_PhiE, const #define sPhi sin(p_PhiE) #define sPsi sin(p_PsiE) - Vector3d result = *this; // default return is this vector - // Define the Rotation Matrix std::vector rotationMatrix = { { cPhi * cPsi - sPhi * cTheta * sPsi , -cPhi * sPsi - sPhi * cTheta * cPsi , sTheta * sPhi }, @@ -90,14 +88,17 @@ Vector3d Vector3d::ChangeBasis(const double p_ThetaE, const double p_PhiE, const { sTheta * sPsi , sTheta * cPsi , cTheta } }; + Vector3d oldVector = *this; // input + Vector3d newVector = Vector3d(0,0,0); // output + // Apply rotation for (size_t row = 0; row < 3; row++) { for (size_t col = 0; col < 3; col++) { - result[row] += result[col] * rotationMatrix[row][col]; + newVector[row] += oldVector[col] * rotationMatrix[row][col]; } } - return result; + return newVector; #undef cTheta #undef cPhi diff --git a/src/vector3d.h b/src/vector3d.h index ed3d86d5d..d450c4834 100644 --- a/src/vector3d.h +++ b/src/vector3d.h @@ -41,8 +41,8 @@ class Vector3d { THROW_ERROR_IF(p_i < 0 || p_i > 2, ERROR::INDEX_OUT_OF_RANGE); // this is a code defect - if (p_i == 1) return m_x; - else if (p_i == 2) return m_y; + if (p_i == 0) return m_x; + else if (p_i == 1) return m_y; else return m_z; } From 6d61e2686541966fb37fab87db66dee51f1e5dfd Mon Sep 17 00:00:00 2001 From: Reinhold Willcox Date: Wed, 27 Nov 2024 17:03:54 +0100 Subject: [PATCH 2/3] Added tweak that sets the Ecc Anomaly for circular systems at first SN (so sets x-basis). Also cleaned up std::couts --- src/BaseBinaryStar.cpp | 28 ++-------------------------- 1 file changed, 2 insertions(+), 26 deletions(-) diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index ff95f28fe..0025303b4 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -1158,8 +1158,6 @@ void BaseBinaryStar::ResolveSupernova() { #define mag Magnitude() #define hat UnitVector() - Vector3d thisVect; - // set relevant preSN parameters m_EccentricityPreSN = m_Eccentricity; m_SemiMajorAxisPreSN = m_SemiMajorAxis; @@ -1172,7 +1170,6 @@ void BaseBinaryStar::ResolveSupernova() { double theta = m_Supernova->SN_Theta(); // angle out of the binary plane double phi = m_Supernova->SN_Phi(); // angle in the binary plane Vector3d natalKickVector = m_Supernova->SN_KickMagnitude() * Vector3d(cos(theta) * cos(phi), cos(theta) * sin(phi), sin(theta)); - thisVect = natalKickVector;std::cout<<"kick vec : "<UpdateComponentVelocity(component1VelocityVectorAtInfinity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); - m_Companion->UpdateComponentVelocity(tmpvec); + m_Companion->UpdateComponentVelocity(component2VelocityVectorAtInfinity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); // Set Euler Angles m_ThetaE = angleBetween(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // angle between the angular momentum unit vectors, always well defined From 10cbbe5c2cb82a1151eaf0920ab8da957a041d3f Mon Sep 17 00:00:00 2001 From: Reinhold Willcox Date: Wed, 27 Nov 2024 17:06:03 +0100 Subject: [PATCH 3/3] changelog --- src/changelog.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/changelog.h b/src/changelog.h index a658ad9b8..3743c765a 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1401,6 +1401,9 @@ // - Fix issue (likely introduced in 03.08.00) with the accretor not gaining mass appropriately // 03.09.01 RTW - Nov 27, 2024 - Enhancement: // - Added systemic velocity components x, y, and z to the output -const std::string VERSION_STRING = "03.09.01"; +// 03.09.02 RTW - Nov 27, 2024 - Defect repair, enhancement: +// - Fixed bugs in vector3d related to indexing and rotation +// - Added tweak for circular systems at first SN, to fix the x-axis along the separation vector +const std::string VERSION_STRING = "03.09.02"; # endif // __changelog_h__