Skip to content

Commit

Permalink
Merge branch 'dev' into dynamical_tides
Browse files Browse the repository at this point in the history
  • Loading branch information
@veome22
veome22 committed Apr 4, 2024
2 parents 3d02ece + bd10dda commit 05f3097
Show file tree
Hide file tree
Showing 17 changed files with 352 additions and 39 deletions.
25 changes: 25 additions & 0 deletions online-docs/pages/User guide/Program options/program-options-list-defaults.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1093,6 +1093,31 @@ Default = FALSE
Print RLOF events to logfile. |br|
Default = TRUE


**--rocket-kick-magnitude-1** |br|
Magnitude of post-SN pulsar rocket kick for the primary, in km/s. |br|
Default = 0.0

**--rocket-kick-magnitude-2** |br|
Magnitude of post-SN pulsar rocket kick for the secondary, in km/s. |br|
Default = 0.0

**--rocket-kick-phi-1** |br|
The in-plane angle [0, 2pi) of the rocket kick velocity that primary neutron star receives following the supernova. |br|
Default = 0.0

**--rocket-kick-phi-2** |br|
The in-plane angle [0, 2pi) of the rocket kick velocity that secondary neutron star receives following the supernova. |br|
Default = 0.0

**--rocket-kick-theta-1** |br|
The polar angle [0, pi] of the rocket kick velocity that primary neutron star receives following the supernova. 0 is aligned with orbital AM. |br|
Default = 0.0

**--rocket-kick-theta-2** |br|
The polar angle [0, pi] of the rocket kick velocity that secondary neutron star receives following the supernova. 0 is aligned with orbital AM. |br|
Default = 0.0

**--rotational-frequency** |br|
Initial rotational frequency of the star for SSE (Hz). |br|
Default = 0.0 (``--rotational-velocity-distribution`` used if ``--rotational-frequency`` not specified)
Expand Down
4 changes: 4 additions & 0 deletions online-docs/pages/whats-new.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,6 +6,10 @@ Following is a brief list of important updates to the COMPAS code. A complete r

**LATEST RELEASE** |br|

**02.43.00 Mar 29, 2024**

* Implementation of the neutrino rocket kick.

**02.42.00 Jan 04, 2023**

* Timesteps are now quantised to an integral multiple of 1e-12Myr.
Expand Down
93 changes: 74 additions & 19 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -387,11 +387,11 @@ void BaseBinaryStar::SetRemainingValues() {
m_Flags.mergesInHubbleTime = false;
m_Unbound = false;

m_SystemicVelocity = Vector3d();
m_OrbitalAngularMomentumVector = Vector3d();
m_ThetaE = DEFAULT_INITIAL_DOUBLE_VALUE;
m_PhiE = DEFAULT_INITIAL_DOUBLE_VALUE;
m_PsiE = DEFAULT_INITIAL_DOUBLE_VALUE;
m_SystemicVelocity = Vector3d();
m_NormalizedOrbitalAngularMomentumVector = Vector3d();
m_ThetaE = DEFAULT_INITIAL_DOUBLE_VALUE;
m_PhiE = DEFAULT_INITIAL_DOUBLE_VALUE;
m_PsiE = DEFAULT_INITIAL_DOUBLE_VALUE;

m_SynchronizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE;
m_CircularizationTimescale = DEFAULT_INITIAL_DOUBLE_VALUE;
Expand Down Expand Up @@ -1181,10 +1181,18 @@ bool BaseBinaryStar::ResolveSupernova() {
Vector3d natalKickVector = m_Supernova->SN_KickMagnitude() *Vector3d(cos(theta) * cos(phi),
cos(theta) * sin(phi),
sin(theta));

// Define the rocket kick vector - will be 0 if unused.
double rocketTheta = m_Supernova->SN_RocketKickTheta(); // Azimuthal angle
double rocketPhi = m_Supernova->SN_RocketKickPhi(); // Polar angle
Vector3d rocketKickVector = m_Supernova->SN_RocketKickMagnitude() * Vector3d( sin(rocketTheta)*cos(rocketPhi),
sin(rocketTheta)*sin(rocketPhi),
cos(rocketTheta)); // The rocket is aligned with the NS spin axis, which by default is aligned with the pre-SN orbit (0.0, 0.0, 1.0). Defined here in case the system is already unbound.

// Check if the system is already unbound
if (IsUnbound()) { // is system already unbound?

m_Supernova->UpdateComponentVelocity( natalKickVector.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // yes - only need to update the velocity of the star undergoing SN
m_Supernova->UpdateComponentVelocity( (natalKickVector+rocketKickVector).ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); // yes - only need to update the velocity of the star undergoing SN

// The quantities below are meaningless in this context, so they are set to nan to avoid misuse
m_OrbitalVelocityPreSN = -nan("");
Expand Down Expand Up @@ -1265,8 +1273,8 @@ bool BaseBinaryStar::ResolveSupernova() {
Vector3d relativeVelocityVector = relativeVelocityVectorPrev + (natalKickVector - companionRecoilVector); // km/s - PostSN relative velocity vector

Vector3d orbitalAngularMomentumVector = cross(separationVectorPrev, relativeVelocityVector); // km^2 s^-1 - PostSN specific orbital angular momentum vector
double orbitalAngularMomentum = orbitalAngularMomentumVector.mag; // km^2 s^-1 - PostSN specific orbital angular momentum
m_OrbitalAngularMomentumVector = orbitalAngularMomentumVector / orbitalAngularMomentum; // set unit vector here to make printing out the inclination vector easier
double orbitalAngularMomentum = orbitalAngularMomentumVector.mag; // km^2 s^-1 - PostSN specific orbital angular momentum
m_NormalizedOrbitalAngularMomentumVector = orbitalAngularMomentumVector/orbitalAngularMomentum; // set unit vector here to make printing out the inclination vector easier

Vector3d eccentricityVector = cross(relativeVelocityVector, orbitalAngularMomentumVector) /
(G_km_Msol_s * totalMass) - separationVectorPrev / separationPrev; // PostSN Laplace-Runge-Lenz vector
Expand All @@ -1282,7 +1290,9 @@ bool BaseBinaryStar::ResolveSupernova() {
// eccentricityVector defines the X'-axis, and
// (orbitalAngularMomentumVector x eccentricityVector) defines the Y'-axis

UpdateSystemicVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE)); // update the system velocity with the new center of mass velocity
UpdateSystemicVelocity(centerOfMassVelocity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); // Update the system velocity with the new center of mass velocity
double reducedMass = m_Supernova->Mass() * m_Companion->Mass() / totalMass; // Reduced Mass
m_Supernova->SetOrbitalEnergyPostSN(CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis)); // Orbital energy

// Split off and evaluate depending on whether the binary is now bound or unbound
if (utils::Compare(m_Eccentricity, 1.0) >= 0) { // unbound?
Expand All @@ -1300,8 +1310,8 @@ bool BaseBinaryStar::ResolveSupernova() {
Vector3d component2VelocityVectorAtInfinity = -(m1 / totalMass) * relativeVelocityVectorAtInfinity + centerOfMassVelocity;

// Update the component velocities
m_Supernova->UpdateComponentVelocity(component1VelocityVectorAtInfinity.RotateVector(m_ThetaE, m_PhiE, m_PsiE));
m_Companion->UpdateComponentVelocity(component2VelocityVectorAtInfinity.RotateVector(m_ThetaE, m_PhiE, m_PsiE));
m_Supernova->UpdateComponentVelocity(component1VelocityVectorAtInfinity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
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
Expand All @@ -1312,10 +1322,10 @@ bool BaseBinaryStar::ResolveSupernova() {

// Set the component velocites to the system velocity. System velocity was already correctly set above.

m_Supernova->UpdateComponentVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE));
m_Companion->UpdateComponentVelocity(centerOfMassVelocity.RotateVector(m_ThetaE, m_PhiE, m_PsiE));
m_Supernova->UpdateComponentVelocity(centerOfMassVelocity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
m_Companion->UpdateComponentVelocity(centerOfMassVelocity.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));

// Calculate Euler angles - see RotateVector() in vector.cpp for details
// Calculate Euler angles - see ChangeBasis() in vector.cpp for details
m_ThetaE = angleBetween(orbitalAngularMomentumVector, orbitalAngularMomentumVectorPrev); // angle between the angular momentum unit vectors, always well defined

// If the new orbital A.M. is parallel or anti-parallel to the previous orbital A.M.,
Expand Down Expand Up @@ -1365,6 +1375,53 @@ bool BaseBinaryStar::ResolveSupernova() {
// the angle of periapsis around the new orbital angular momentum, (i.e, Psi) - RTW 15/05/20
m_PsiE = _2_PI * RAND->Random();
}

// Account for possible neutrino rocket - see Hirai+ 2024
if (ShouldResolveNeutrinoRocketMechanism()) {

if (IsUnbound()) { // Is system unbound?
m_Supernova->UpdateComponentVelocity(rocketKickVector.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE)); // yes - simply update the component velocity
} else { // no - need to update the eccentricity and system velocity

Vector3d eccentricityVectorPreRocket = eccentricityVector; // defined earlier
double averageOrbitalVelocityPreRocket = std::sqrt( -2*m_OrbitalEnergy/reducedMass); // AU/yr - average orbital velocity post-SN
double k_grav = averageOrbitalVelocityPreRocket*averageOrbitalVelocityPreRocket
* reducedMass * m_SemiMajorAxis; // AU^3 * Msol / yr^2
Vector3d totalAmVectorPreRocket = orbitalAngularMomentumVector * reducedMass
* KM_TO_AU * KM_TO_AU * SECONDS_IN_YEAR; // Msol * AU^2 / yr (orbitalAngularMomentumVector is the specific orbital AM)
Vector3d amVectorNormalizedByCircularAmPreRocket = totalAmVectorPreRocket
*(averageOrbitalVelocityPreRocket / k_grav) ; // unitless!
double theta_rotation = 3*rocketKickVector.mag * KM_TO_AU * SECONDS_IN_YEAR
/ (2*averageOrbitalVelocityPreRocket); // rad - need to convert velocities to same units

// Apply hPlus and hMinus support vectors
Vector3d hPlusVector = amVectorNormalizedByCircularAmPreRocket + eccentricityVectorPreRocket;
Vector3d hMinusVector = amVectorNormalizedByCircularAmPreRocket - eccentricityVectorPreRocket;

// Rotate hPlus and hMinus vectors so that the thrust is parallel to the z-axis, in order to apply the rotation below
hPlusVector = hPlusVector.RotateVectorAboutZ( -rocketPhi).RotateVectorAboutY(-rocketTheta);
hMinusVector = hMinusVector.RotateVectorAboutZ(-rocketPhi).RotateVectorAboutY(-rocketTheta);

// Rotate vectors about the new "z-axis" - parallel to the rocket thrust
Vector3d hPlusVector_prime = hPlusVector.RotateVectorAboutZ( theta_rotation );
Vector3d hMinusVector_prime = hMinusVector.RotateVectorAboutZ( -theta_rotation );

// Rotate new hPlus and hMinus vectors back to the original frame
hPlusVector = hPlusVector.RotateVectorAboutY( rocketTheta).RotateVectorAboutZ(rocketPhi);
hMinusVector = hMinusVector.RotateVectorAboutY(rocketTheta).RotateVectorAboutZ(rocketPhi);

// Calculate post-rocket values
Vector3d normalizedAngularMomentumVectorPostRocket = 0.5 * (hPlusVector_prime + hMinusVector_prime);
Vector3d eccentricityVectorPostRocket = 0.5 * (hPlusVector_prime - hMinusVector_prime);

m_NormalizedOrbitalAngularMomentumVector = normalizedAngularMomentumVectorPostRocket ;
m_Eccentricity = eccentricityVectorPostRocket.mag;

UpdateSystemicVelocity(rocketKickVector.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
m_Supernova->UpdateComponentVelocity(rocketKickVector.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
m_Companion->UpdateComponentVelocity(rocketKickVector.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
}
}

#undef hat
#undef mag
Expand All @@ -1373,12 +1430,10 @@ bool BaseBinaryStar::ResolveSupernova() {
#undef cross
}

// Set remaining post-SN values
double totalMass = m_Supernova->Mass() + m_Companion->Mass(); // total Mass
double reducedMass = m_Supernova->Mass() * m_Companion->Mass() / totalMass; // reduced Mass
m_Supernova->SetOrbitalEnergyPostSN(CalculateOrbitalEnergy(reducedMass, totalMass, m_SemiMajorAxis)); // orbital energy
//////////////////////////
// Do for all systems

m_IPrime = m_ThetaE; // inclination angle between preSN and postSN orbital planes
m_IPrime = m_ThetaE; // Inclination angle between preSN and postSN orbital planes
m_CosIPrime = cos(m_IPrime);

(void)PrintSupernovaDetails(); // log record to supernovae logfile
Expand Down
13 changes: 8 additions & 5 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -101,7 +101,7 @@ class BaseBinaryStar {
m_SynchronizationTimescale = p_Star.m_SynchronizationTimescale;

m_SystemicVelocity = p_Star.m_SystemicVelocity;
m_OrbitalAngularMomentumVector = p_Star.m_OrbitalAngularMomentumVector;
m_NormalizedOrbitalAngularMomentumVector = p_Star.m_NormalizedOrbitalAngularMomentumVector;
m_ThetaE = p_Star.m_ThetaE;
m_PhiE = p_Star.m_PhiE;
m_PsiE = p_Star.m_PsiE;
Expand Down Expand Up @@ -257,9 +257,9 @@ class BaseBinaryStar {
STELLAR_TYPE StellarType2PostCEE() const { return m_Star2->StellarTypePostCEE(); }
STELLAR_TYPE StellarType2PreCEE() const { return m_Star2->StellarTypePreCEE(); }
double SN_OrbitInclinationAngle() const { return m_ThetaE; }
double SN_OrbitInclinationVectorX() const { return m_OrbitalAngularMomentumVector.xValue(); }
double SN_OrbitInclinationVectorY() const { return m_OrbitalAngularMomentumVector.yValue(); }
double SN_OrbitInclinationVectorZ() const { return m_OrbitalAngularMomentumVector.zValue(); }
double SN_OrbitInclinationVectorX() const { return m_NormalizedOrbitalAngularMomentumVector.xValue(); }
double SN_OrbitInclinationVectorY() const { return m_NormalizedOrbitalAngularMomentumVector.yValue(); }
double SN_OrbitInclinationVectorZ() const { return m_NormalizedOrbitalAngularMomentumVector.zValue(); }
SN_STATE SN_State() const { return m_SupernovaState; }
double SynchronizationTimescale() const { return m_SynchronizationTimescale; }
double SystemicSpeed() const { return m_SystemicVelocity.Magnitude(); }
Expand Down Expand Up @@ -367,7 +367,7 @@ class BaseBinaryStar {
double m_SynchronizationTimescale;

Vector3d m_SystemicVelocity; // Systemic velocity vector, relative to ZAMS Center of Mass
Vector3d m_OrbitalAngularMomentumVector; // Orbital AM vector postSN, in preSN frame
Vector3d m_NormalizedOrbitalAngularMomentumVector; // Orbital AM vector postSN, in preSN frame
double m_ThetaE; // Euler Theta
double m_PhiE; // Euler Phi
double m_PsiE; // Euler Psi
Expand Down Expand Up @@ -517,6 +517,9 @@ class BaseBinaryStar {
return LOGGING->LogBSESupernovaDetails(this, p_RecordType);
}

bool ShouldResolveNeutrinoRocketMechanism() const {
return (OPTIONS->RocketKickMagnitude1() > 0) || (OPTIONS->RocketKickMagnitude2() > 0);
}

/*
* Functor for MassLossToFitInsideRocheLobe()
Expand Down
Loading

0 comments on commit 05f3097

Please sign in to comment.