SN kick fixes

src/BaseBinaryStar.cpp

@@ -1201,6 +1201,10 @@ void BaseBinaryStar::ResolveSupernova() {
         m_Supernova->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;

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__

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,

src/vector3d.cpp

@@ -81,23 +81,24 @@ 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<DBL_VECTOR> rotationMatrix = {
         { cPhi * cPsi - sPhi * cTheta * sPsi ,  -cPhi * sPsi - sPhi * cTheta * cPsi ,  sTheta * sPhi },
         { sPhi * cPsi + cPhi * cTheta * sPsi ,  -sPhi * sPsi + cPhi * cTheta * cPsi , -sTheta * cPhi },
         { 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

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;
     }