inf (divide-by-zero) in BaseBinaryStar::ResolveSupernova() when kick magnitude is 0.0

[jeffriley]

When kick magnitude is 0.0 (e.g. --kick-magnitude-distribution ZERO), this conditional statement in BaseBinaryStar::ResolveSupernova():

1338 m_PsiE = utils::Compare( dot(eccentricityVector, orbitalAngularMomentumVectorPrev), 0.0) >= 0 // are eccentricityVector and orbitalAngularMomentumVectorPrev in the same hemisphere?
1339     ? angleBetween(eccentricityVector, orbitalPivotAxis)                                      // yes - psi in [0,pi)
1340     : -angleBetween(eccentricityVector, orbitalPivotAxis);                                    // no  - psi in [-pi,0)

results in a divide-by-zero in Vector3d::AngleBetween(), because orbitalPivotAxis is a zero vector.

orbitalPivotAxis is calculated here in BaseBinaryStar::ResolveSupernova():

1323 Vector3d orbitalPivotAxis = cross(orbitalAngularMomentumVectorPrev, orbitalAngularMomentumVector); // cross product of the orbit normals

orbitalAngularMomentumVectorPrev and orbitalAngularMomentumVector are calculated here in BaseBinaryStar::ResolveSupernova():

1219 Vector3d orbitalAngularMomentumVectorPrev = cross(separationVectorPrev, relativeVelocityVectorPrev); // specific orbital angular momentum vector (km^2 s^-1)
...
1249 Vector3d orbitalAngularMomentumVector = cross(separationVectorPrev, relativeVelocityVector);              // post-SN specific orbital angular momentum vector (km^2 s^-1)

with the difference in the calculations being whether relativeVelocityVectorPrev or relativeVelocityVector is used.

relativeVelocityVectorPrev is calculated here in BaseBinaryStar::ResolveSupernova():

1218 Vector3d relativeVelocityVectorPrev = Vector3d(-fact1 * sinEccAnomaly, fact1 * cosEccAnomaly * sqrt1MinusEccPrevSquared, 0.0); // relative velocity vector, in the m1Prev rest frame (km/s)

and relativeVelocityVector is calculated here:

1247 Vector3d relativeVelocityVector = relativeVelocityVectorPrev + (natalKickVector - companionRecoilVector); // post-SN relative velocity vector (km/s)

natalKickVector is a zero vector (kick is zero), and companionRecoilVector is defined here as a zero vector:

1235 Vector3d companionRecoilVector = Vector3d(0.0, 0.0, 0.0); // km/s - The recoil of the companion due to ablation

and is never changed - so relativeVelocityVectorPrev and relativeVelocityVector are the same vector, which results in orbitalAngularMomentumVectorPrev and orbitalAngularMomentumVector being the same vector, so orbitalPivotAxis is just the cross product of a vector with itself, resulting in a zero vector, hence the divide-by-zero when calculating m_PsiE

Also note that companionRecoilVector, which is set to zero and never changed, is used in the calculation of centerOfMassVelocity here:

1244 Vector3d centerOfMassVelocity   = (-m2Prev * dm1 / fact2 + m1Prev * dm2 / fact2) * relativeVelocityVectorPrev + 
1245                                   (m1 / totalMass) * natalKickVector + (m2 / totalMass) * companionRecoilVector; // post-SN center of mass velocity vector (km/s)

For clarity, or future updates?

COMPAS version: v03.07.00
OS version: Ubuntu 22.04

I imagine that there is some code in BaseBinayStar::ResolveSupernova() that we can skip, or at least simplify, if the kick magnitude is 0.0, but I'll leave that to someone more familiar with the subject (@reinhold-willcox?). However we resolve this issue, we should not be dividing by zero.

To Reproduce
Run COMPAS with:

./COMPAS --critical-mass-ratio-prescription NONE --kick-magnitude-distribution ZERO --random-seed 1724542090 --common-envelope-alpha 1.6265542767303343 --wolf-rayet-multiplier 0.8084912681292984 --number-of-systems 1 --initial-mass-function KROUPA --mass-ratio-distribution DUQUENNOYMAYOR1991 --semi-major-axis-distribution SANA2012 --minimum-secondary-mass 5 --metallicity-distribution LOGUNIFORM --tides-prescription PERFECT --remnant-mass-prescription FRYER2022 --common-envelope-allow-radiative-envelope-survive False --common-envelope-allow-immediate-RLOF-post-CE-survive False

Expected behavior
Production code does not cause inf

**Versioning

OS: Ubuntu 22.04
COMPAS v03.07.01

[ilyamandel]

@reinhold-willcox ,
There's clearly a logic issue here. You do try to check for the situation when the current and previous angular momentum are aligned or anti-aligned (m_ThetaE is either 0 or pi) -- but you only seem to care if that's true and the orbit was and remains eccentric, and ignore the situation of a circular orbit whose plane has not changed...