diff --git a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst index a01d3f55e..c3e988bcb 100644 --- a/online-docs/pages/User guide/Program options/program-options-list-defaults.rst +++ b/online-docs/pages/User guide/Program options/program-options-list-defaults.rst @@ -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) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 13e4d7952..a03172f6c 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -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. diff --git a/src/BaseBinaryStar.cpp b/src/BaseBinaryStar.cpp index 87abf60d9..01e93a127 100644 --- a/src/BaseBinaryStar.cpp +++ b/src/BaseBinaryStar.cpp @@ -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; @@ -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(""); @@ -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 @@ -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? @@ -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 @@ -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., @@ -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 @@ -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 diff --git a/src/BaseBinaryStar.h b/src/BaseBinaryStar.h index 4f90e3f27..faa773bcd 100644 --- a/src/BaseBinaryStar.h +++ b/src/BaseBinaryStar.h @@ -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; @@ -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(); } @@ -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 @@ -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() diff --git a/src/BaseStar.cpp b/src/BaseStar.cpp index c1b8fe2ac..c9a032af8 100755 --- a/src/BaseStar.cpp +++ b/src/BaseStar.cpp @@ -200,6 +200,10 @@ BaseStar::BaseStar(const unsigned long int p_RandomSeed, m_SupernovaDetails.phi = p_KickParameters.phi; m_SupernovaDetails.meanAnomaly = p_KickParameters.meanAnomaly; + m_SupernovaDetails.rocketKickMagnitude = DEFAULT_INITIAL_DOUBLE_VALUE; + m_SupernovaDetails.rocketKickPhi = DEFAULT_INITIAL_DOUBLE_VALUE; + m_SupernovaDetails.rocketKickTheta = DEFAULT_INITIAL_DOUBLE_VALUE; + // Calculates the Baryonic mass for which the GravitationalRemnantMass will be equal to the maximumNeutronStarMass (inverse of SolveQuadratic()) // needed to decide whether to calculate Fryer+2012 for Neutron Star or Black Hole in GiantBranch::CalculateGravitationalRemnantMass() // calculate only once for entire simulation of N binaries in the future. @@ -361,6 +365,9 @@ COMPAS_VARIABLE BaseStar::StellarPropertyValue(const T_ANY_PROPERTY p_Property) case ANY_STAR_PROPERTY::RADIAL_EXPANSION_TIMESCALE: value = CalculateRadialExpansionTimescale(); break; case ANY_STAR_PROPERTY::RADIUS: value = Radius(); break; case ANY_STAR_PROPERTY::RANDOM_SEED: value = RandomSeed(); break; + case ANY_STAR_PROPERTY::ROCKET_KICK_MAGNITUDE: value = SN_RocketKickMagnitude(); break; + case ANY_STAR_PROPERTY::ROCKET_KICK_PHI: value = SN_RocketKickPhi(); break; + case ANY_STAR_PROPERTY::ROCKET_KICK_THETA: value = SN_RocketKickTheta(); break; case ANY_STAR_PROPERTY::RZAMS: value = RZAMS(); break; case ANY_STAR_PROPERTY::SN_TYPE: value = SN_Type(); break; case ANY_STAR_PROPERTY::SPEED: value = Speed(); break; diff --git a/src/BaseStar.h b/src/BaseStar.h index b5a5c7f1c..551e22da2 100644 --- a/src/BaseStar.h +++ b/src/BaseStar.h @@ -126,6 +126,9 @@ class BaseStar { double SN_KickMagnitude() const { return m_SupernovaDetails.kickMagnitude; } double SN_MeanAnomaly() const { return m_SupernovaDetails.meanAnomaly; } double SN_Phi() const { return m_SupernovaDetails.phi; } + double SN_RocketKickMagnitude() const { return m_SupernovaDetails.rocketKickMagnitude; } + double SN_RocketKickPhi() const { return m_SupernovaDetails.rocketKickPhi; } + double SN_RocketKickTheta() const { return m_SupernovaDetails.rocketKickTheta; } double SN_TotalMassAtCOFormation() const { return m_SupernovaDetails.totalMassAtCOFormation; } double SN_TrueAnomaly() const { return m_SupernovaDetails.trueAnomaly; } double SN_Theta() const { return m_SupernovaDetails.theta; } diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index 936d76f4d..842c73378 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -2074,6 +2074,9 @@ STELLAR_TYPE GiantBranch::ResolveSupernova() { } CalculateSNKickMagnitude(m_Mass, m_SupernovaDetails.totalMassAtCOFormation - m_Mass, stellarType); + if ( !utils::IsOneOf(stellarType, { STELLAR_TYPE::NEUTRON_STAR })) { + m_SupernovaDetails.rocketKickMagnitude = 0; // Only NSs can get rocket kicks + } // stash SN details for later printing to the SSE Supernova log // can't print it now because we may revert state (in Star::EvolveOneTimestep()) diff --git a/src/Options.cpp b/src/Options.cpp index 7800ae93e..46151626e 100644 --- a/src/Options.cpp +++ b/src/Options.cpp @@ -303,6 +303,13 @@ void Options::OptionValues::Initialise() { m_BlackHoleKicks.type = BLACK_HOLE_KICKS::FALLBACK; m_BlackHoleKicks.typeString = BLACK_HOLE_KICKS_LABEL.at(m_BlackHoleKicks.type); + // Rocket kicks + m_RocketKickMagnitude1 = 0.0; + m_RocketKickMagnitude2 = 0.0; + m_RocketKickPhi1 = 0.0; + m_RocketKickPhi2 = 0.0; + m_RocketKickTheta1 = 0.0; + m_RocketKickTheta2 = 0.0; // Chemically Homogeneous Evolution m_CheMode.type = CHE_MODE::PESSIMISTIC; @@ -1196,17 +1203,17 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt ( "kick-magnitude", po::value(&p_Options->m_KickMagnitude)->default_value(p_Options->m_KickMagnitude), - ("The magnitude of the kick velocity, in km/s, that the star receives during the a supernova (default = " + std::to_string(p_Options->m_KickMagnitude) + ")").c_str() + ("The magnitude of the kick velocity, in km/s, that the star receives during the supernova (default = " + std::to_string(p_Options->m_KickMagnitude) + ")").c_str() ) ( "kick-magnitude-1", po::value(&p_Options->m_KickMagnitude1)->default_value(p_Options->m_KickMagnitude1), - ("The magnitude of the kick velocity, in km/s, that the primary star receives during the a supernova (default = " + std::to_string(p_Options->m_KickMagnitude1) + ")").c_str() + ("The magnitude of the kick velocity, in km/s, that the primary star receives during the supernova (default = " + std::to_string(p_Options->m_KickMagnitude1) + ")").c_str() ) ( "kick-magnitude-2", po::value(&p_Options->m_KickMagnitude2)->default_value(p_Options->m_KickMagnitude2), - ("The magnitude of the kick velocity, in km/s, that the secondary star receives during the a supernova (default = " + std::to_string(p_Options->m_KickMagnitude2) + ")").c_str() + ("The magnitude of the kick velocity, in km/s, that the secondary star receives during the supernova (default = " + std::to_string(p_Options->m_KickMagnitude2) + ")").c_str() ) ( "kick-magnitude-max", @@ -1470,6 +1477,37 @@ bool Options::AddOptions(OptionValues *p_Options, po::options_description *p_Opt ("Minimum pulsar magnetic field, in log10(Gauss) (default = " + std::to_string(p_Options->m_PulsarLog10MinimumMagneticField) + ")").c_str() ) + ( + "rocket-kick-magnitude-1", + po::value(&p_Options->m_RocketKickMagnitude1)->default_value(p_Options->m_RocketKickMagnitude1), + ("The magnitude of the rocket kick velocity, in km/s, that primary neutron star receives following the supernova (default = " + std::to_string(p_Options->m_RocketKickMagnitude1) + ")").c_str() + ) + ( + "rocket-kick-magnitude-2", + po::value(&p_Options->m_RocketKickMagnitude2)->default_value(p_Options->m_RocketKickMagnitude2), + ("The magnitude of the rocket kick velocity, in km/s, that primary neutron star receives following the supernova (default = " + std::to_string(p_Options->m_RocketKickMagnitude2) + ")").c_str() + ) + ( + "rocket-kick-phi-1", + po::value(&p_Options->m_RocketKickPhi1)->default_value(p_Options->m_RocketKickPhi1), + ("The in-plane angle [0, 2pi) of the rocket kick velocity that primary neutron star receives following the supernova (default = " + std::to_string(p_Options->m_RocketKickPhi1) + ")").c_str() + ) + ( + "rocket-kick-phi-2", + po::value(&p_Options->m_RocketKickPhi2)->default_value(p_Options->m_RocketKickPhi2), + ("The in-plane angle [0, 2pi) of the rocket kick velocity that secondary neutron star receives following the supernova (default " + std::to_string(p_Options->m_RocketKickPhi2) + ")").c_str() + ) + ( + "rocket-kick-theta-1", + po::value(&p_Options->m_RocketKickTheta1)->default_value(p_Options->m_RocketKickTheta1), + ("The polar angle [0, pi] of the rocket kick velocity that primary neutron star receives following the supernova. 0 is aligned with orbital AM (default = " + std::to_string(p_Options->m_RocketKickTheta1) + ")").c_str() + ) + ( + "rocket-kick-theta-2", + po::value(&p_Options->m_RocketKickTheta2)->default_value(p_Options->m_RocketKickTheta2), + ("The polar angle [0, pi] of the rocket kick velocity that secondary neutron star receives following the supernova. 0 is aligned with orbital AM (default = " + std::to_string(p_Options->m_RocketKickTheta2) + ")").c_str() + ) + ( "rotational-frequency", po::value(&p_Options->m_RotationalFrequency)->default_value(p_Options->m_RotationalFrequency), @@ -2345,6 +2383,8 @@ std::string Options::OptionValues::CheckAndSetOptions() { COMPLAIN_IF(m_KickMagnitude < 0.0, "Kick magnitude (--kick-magnitude) must be >= 0"); COMPLAIN_IF(m_KickMagnitude1 < 0.0, "Kick magnitude (--kick-magnitude-1) must be >= 0"); COMPLAIN_IF(m_KickMagnitude2 < 0.0, "Kick magnitude (--kick-magnitude-2) must be >= 0"); + COMPLAIN_IF(m_RocketKickMagnitude1 < 0.0, "Rocket Kick magnitude 1 (--rocket-kick-magnitude-1) must be >= 0"); + COMPLAIN_IF(m_RocketKickMagnitude2 < 0.0, "Rocket Kick magnitude 2 (--rocket-kick-magnitude-2) must be >= 0"); COMPLAIN_IF(m_KickMagnitudeRandom < 0.0 || m_KickMagnitudeRandom >= 1.0, "Kick magnitude random (--kick-magnitude-random) must be >= 0 and < 1"); COMPLAIN_IF(m_KickMagnitudeRandom1 < 0.0 || m_KickMagnitudeRandom1 >= 1.0, "Kick magnitude random (--kick-magnitude-random-1) must be >= 0 and < 1"); @@ -4587,6 +4627,13 @@ COMPAS_VARIABLE Options::OptionValue(const T_ANY_PROPERTY p_Property) const { case PROGRAM_OPTION::REMNANT_MASS_PRESCRIPTION : value = static_cast(RemnantMassPrescription()); break; + case PROGRAM_OPTION::ROCKET_KICK_MAGNITUDE_1 : value = RocketKickMagnitude1(); break; + case PROGRAM_OPTION::ROCKET_KICK_MAGNITUDE_2 : value = RocketKickMagnitude2(); break; + case PROGRAM_OPTION::ROCKET_KICK_PHI_1 : value = RocketKick_Phi1(); break; + case PROGRAM_OPTION::ROCKET_KICK_PHI_2 : value = RocketKick_Phi2(); break; + case PROGRAM_OPTION::ROCKET_KICK_THETA_1 : value = RocketKick_Theta1(); break; + case PROGRAM_OPTION::ROCKET_KICK_THETA_2 : value = RocketKick_Theta2(); break; + case PROGRAM_OPTION::ROTATIONAL_VELOCITY_DISTRIBUTION : value = static_cast(RotationalVelocityDistribution()); break; case PROGRAM_OPTION::ROTATIONAL_FREQUENCY : value = RotationalFrequency(); break; diff --git a/src/Options.h b/src/Options.h index 582cdf4db..be6544204 100755 --- a/src/Options.h +++ b/src/Options.h @@ -388,6 +388,13 @@ class Options { "rotational-frequency-1", "rotational-frequency-2", + "rocket-kick-magnitude-1", + "rocket-kick-magnitude-2", + "rocket-kick-phi-1", + "rocket-kick-phi-2", + "rocket-kick-theta-1", + "rocket-kick-theta-2", + "semi-major-axis", "a", "semi-major-axis-distribution", "semi-major-axis-max", @@ -767,6 +774,14 @@ class Options { // Black hole kicks ENUM_OPT m_BlackHoleKicks; // Which black hole kicks mode + // Rocket kicks + double m_RocketKickMagnitude1; // Rocket kick magnitude primary - only for neutron stars + double m_RocketKickMagnitude2; // Rocket kick magnitude secondary - only for neutron stars + double m_RocketKickPhi1; // Rocket kick phi angle primary + double m_RocketKickPhi2; // Rocket kick phi angle secondary + double m_RocketKickTheta1; // Rocket kick theta angle primary + double m_RocketKickTheta2; // Rocket kick theta angle secondary + // CHE - Chemically Homogeneous Evolution ENUM_OPT m_CheMode; // Which Chemically Homogeneous Evolution mode @@ -1454,6 +1469,13 @@ class Options { bool RLOFPrinting() const { return m_CmdLine.optionValues.m_RlofPrinting; } + double RocketKickMagnitude1() const { return OPT_VALUE("rocket-kick-magnitude-1", m_RocketKickMagnitude1, true); } + double RocketKickMagnitude2() const { return OPT_VALUE("rocket-kick-magnitude-2", m_RocketKickMagnitude2, true); } + double RocketKick_Phi1() const { return OPT_VALUE("rocket-kick-phi-1", m_RocketKickPhi1, true); } + double RocketKick_Phi2() const { return OPT_VALUE("rocket-kick-phi-2", m_RocketKickPhi2, true); } + double RocketKick_Theta1() const { return OPT_VALUE("rocket-kick-theta-1", m_RocketKickTheta1, true); } + double RocketKick_Theta2() const { return OPT_VALUE("rocket-kick-theta-2", m_RocketKickTheta2, true); } + ROTATIONAL_VELOCITY_DISTRIBUTION RotationalVelocityDistribution() const { return OPT_VALUE("rotational-velocity-distribution", m_RotationalVelocityDistribution.type, true); } double RotationalFrequency() const { return OPT_VALUE("rotational-frequency", m_RotationalFrequency, true); } double RotationalFrequency1() const { return OPT_VALUE("rotational-frequency-1", m_RotationalFrequency1, true); } diff --git a/src/Star.h b/src/Star.h index fc1c79416..1670b807a 100755 --- a/src/Star.h +++ b/src/Star.h @@ -226,6 +226,10 @@ class Star { void SetSNPastEvent(const SN_EVENT p_SNEvent) { m_Star->SetSNPastEvent(p_SNEvent); } double SN_KickMagnitude() { return m_Star->SN_KickMagnitude() ; } + double SN_RocketKickMagnitude() { return m_Star->SN_RocketKickMagnitude(); } + double SN_RocketKickPhi() { return m_Star->SN_RocketKickPhi(); } + double SN_RocketKickTheta() { return m_Star->SN_RocketKickTheta(); } + STELLAR_TYPE SwitchTo(const STELLAR_TYPE p_StellarType, bool p_SetInitialType = false); diff --git a/src/changelog.h b/src/changelog.h index 262acc8bf..c2ba05dcb 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1113,7 +1113,9 @@ // - Defect repair : Added explicit definition `bool isUnstable = false` to avoid confusion in BaseBinaryStar.cpp // - Defect repair : Fixed erroneous core mass values in ResolveSNIa in WhiteDwarfs.cpp. Was previously 0 for all core masses. // - Enhancement: Added output parameter TZAMS for internal variable m_TZAMS +// 02.43.00 RTW - Mar 29, 2023 - Enhancement: +// - Added Hirai pulsar rocket kick, and related options -const std::string VERSION_STRING = "02.42.02"; +const std::string VERSION_STRING = "02.43.00"; # endif // __changelog_h__ diff --git a/src/compasConfigDefault.yaml b/src/compasConfigDefault.yaml index 0907aa81b..611e97819 100644 --- a/src/compasConfigDefault.yaml +++ b/src/compasConfigDefault.yaml @@ -1,5 +1,5 @@ ##~!!~## COMPAS option values -##~!!~## File Created Wed Nov 1 11:22:29 2023 by COMPAS v02.40.00 +##~!!~## File Created Tue Apr 2 15:10:06 2024 by COMPAS v02.43.00 ##~!!~## ##~!!~## The default COMPAS YAML file (``compasConfigDefault.yaml``), as distributed, has ##~!!~## all COMPAS option entries commented so that the COMPAS default value for the @@ -47,6 +47,9 @@ booleanChoices: # --common-envelope-lambda-nanjing-use-rejuvenated-mass: False # Default: False # --revised-energy-formalism-nandez-ivanova: False # Default: False + ### TIDES +# --enable-tides: False # Default: False + ### SUPERNOVAE, KICKS AND REMNANTS # --allow-non-stripped-ECSN: False # Default: False # --pair-instability-supernovae: True # Default: True @@ -187,6 +190,12 @@ numericalChoices: # --pisn-upper-limit: 135.00 # Default: 135.00 # Maximum core mass for PISN [Msol] # --ppi-lower-limit: 35.000000 # Default: 35.000000 # Minimum core mass for PPI [Msol] # --ppi-upper-limit: 60.000000 # Default: 60.000000 # Maximum core mass for PPI [Msol] +# --rocket-kick-magnitude-1: 0.000000 # Default: 0.000000 # (BSE) rocket kick magnitude for the primary should it undergo a supernova event [km/s] +# --rocket-kick-magnitude-2: 0.000000 # Default: 0.000000 # (BSE) rocket kick magnitude for the secondary should it undergo a supernova event [km/s] +# --rocket-kick-theta-1: 0.000000 # Default: 0.000000 # (BSE) angle between the orbital plane and the 'z' axis of the rocket kick for the primary star [radians] +# --rocket-kick-phi-1: 0.000000 # Default: 0.000000 # (BSE) angle between 'x' and 'y', both in the orbital plane of the rocket kick , for the primary star [radians] +# --rocket-kick-theta-2: 0.000000 # Default: 0.000000 # (BSE) angle between the orbital plane and the 'z' axis of the rocket kick for the secondary star [radians] +# --rocket-kick-phi-2: 0.000000 # Default: 0.000000 # (BSE) angle between 'x' and 'y', both in the orbital plane of the rocket kick , for the secondary star [radians] ### PULSAR PARAMETERS # --pulsar-birth-magnetic-field-distribution-min: 11.000000 # Default: 11.000000 # [log10(B/G)] @@ -207,6 +216,7 @@ stringChoices: # --notes: { } # Default: { } # --notes-hdrs: { } # Default: { } # --output-container: 'COMPAS_Output' # Default: 'COMPAS_Output' +# --timesteps-filename: '' # Default: '' ### STELLAR PROPERTIES # --chemically-homogeneous-evolution: 'PESSIMISTIC' # Default: 'PESSIMISTIC' # Options: ['PESSIMISTIC','OPTIMISTIC','NONE'] # chemically homogeneous evolution diff --git a/src/constants.h b/src/constants.h index 96f9122b9..8515a1643 100755 --- a/src/constants.h +++ b/src/constants.h @@ -1919,6 +1919,9 @@ const COMPASUnorderedMap PROPERTY_TYPE_LABEL = { RANDOM_SEED, \ RECYCLED_NEUTRON_STAR, \ RLOF_ONTO_NS, \ + ROCKET_KICK_MAGNITUDE, \ + ROCKET_KICK_PHI, \ + ROCKET_KICK_THETA, \ RZAMS, \ SN_TYPE, \ SPEED, \ @@ -2075,6 +2078,9 @@ const COMPASUnorderedMap STAR_PROPERTY_LABEL = { { STAR_PROPERTY::RANDOM_SEED, "RANDOM_SEED" }, { STAR_PROPERTY::RECYCLED_NEUTRON_STAR, "RECYCLED_NEUTRON_STAR" }, { STAR_PROPERTY::RLOF_ONTO_NS, "RLOF_ONTO_NS" }, + { STAR_PROPERTY::ROCKET_KICK_MAGNITUDE, "ROCKET_KICK_MAGNITUDE" }, + { STAR_PROPERTY::ROCKET_KICK_PHI, "ROCKET_KICK_PHI" }, + { STAR_PROPERTY::ROCKET_KICK_THETA, "ROCKET_KICK_THETA" }, { STAR_PROPERTY::RZAMS, "RZAMS" }, { STAR_PROPERTY::SN_TYPE, "SN_TYPE" }, { STAR_PROPERTY::SPEED, "SPEED" }, @@ -2586,6 +2592,13 @@ enum class PROGRAM_OPTION: int { REMNANT_MASS_PRESCRIPTION, + ROCKET_KICK_MAGNITUDE_1, + ROCKET_KICK_MAGNITUDE_2, + ROCKET_KICK_PHI_1, + ROCKET_KICK_PHI_2, + ROCKET_KICK_THETA_1, + ROCKET_KICK_THETA_2, + ROTATIONAL_VELOCITY_DISTRIBUTION, ROTATIONAL_FREQUENCY, ROTATIONAL_FREQUENCY_1, @@ -2796,6 +2809,13 @@ const COMPASUnorderedMap PROGRAM_OPTION_LABEL = { { PROGRAM_OPTION::REMNANT_MASS_PRESCRIPTION, "REMNANT_MASS_PRESCRIPTION" }, + { PROGRAM_OPTION::ROCKET_KICK_MAGNITUDE_1, "ROCKET_KICK_MAGNITUDE_1" }, + { PROGRAM_OPTION::ROCKET_KICK_MAGNITUDE_2, "ROCKET_KICK_MAGNITUDE_2" }, + { PROGRAM_OPTION::ROCKET_KICK_PHI_1, "ROCKET_KICK_PHI_1" }, + { PROGRAM_OPTION::ROCKET_KICK_PHI_2, "ROCKET_KICK_PHI_2" }, + { PROGRAM_OPTION::ROCKET_KICK_THETA_1, "ROCKET_KICK_THETA_1" }, + { PROGRAM_OPTION::ROCKET_KICK_THETA_2, "ROCKET_KICK_THETA_2" }, + { PROGRAM_OPTION::ROTATIONAL_VELOCITY_DISTRIBUTION, "ROTATIONAL_VELOCITY_DISTRIBUTION" }, { PROGRAM_OPTION::ROTATIONAL_FREQUENCY, "ROTATIONAL_FREQUENCY" }, { PROGRAM_OPTION::ROTATIONAL_FREQUENCY_1, "ROTATIONAL_FREQUENCY_1" }, @@ -2963,6 +2983,9 @@ const std::map ANY_STAR_PROPERTY_DETAIL = { { ANY_STAR_PROPERTY::RANDOM_SEED, { TYPENAME::ULONGINT, "SEED", "-", 12, 1 }}, { ANY_STAR_PROPERTY::RECYCLED_NEUTRON_STAR, { TYPENAME::BOOL, "Recycled_NS", "Event", 0, 0 }}, { ANY_STAR_PROPERTY::RLOF_ONTO_NS, { TYPENAME::BOOL, "RLOF->NS", "Event", 0, 0 }}, + { ANY_STAR_PROPERTY::ROCKET_KICK_MAGNITUDE, { TYPENAME::DOUBLE, "Rocket_Kick_Magnitude", "kms^-1", 14, 6 }}, + { ANY_STAR_PROPERTY::ROCKET_KICK_PHI, { TYPENAME::DOUBLE, "Rocket_Kick_Phi", "-", 14, 6 }}, + { ANY_STAR_PROPERTY::ROCKET_KICK_THETA, { TYPENAME::DOUBLE, "Rocket_Kick_Theta", "-", 14, 6 }}, { ANY_STAR_PROPERTY::RZAMS, { TYPENAME::DOUBLE, "Radius@ZAMS", "Rsol", 14, 6 }}, { ANY_STAR_PROPERTY::SN_TYPE, { TYPENAME::SN_EVENT, "SN_Type", "-", 4, 1 }}, { ANY_STAR_PROPERTY::SPEED, { TYPENAME::DOUBLE, "ComponentSpeed", "kms^-1", 14, 6 }}, @@ -3305,6 +3328,13 @@ const std::map PROGRAM_OPTION_DETAIL = { { PROGRAM_OPTION::REMNANT_MASS_PRESCRIPTION, { TYPENAME::INT, "Remnant_Mass_Prscrptn", "-", 4, 1 }}, + { PROGRAM_OPTION::ROCKET_KICK_MAGNITUDE_1, { TYPENAME::DOUBLE, "Rocket_Kick_Magnitude(1)", "kms^-1", 14, 6 }}, + { PROGRAM_OPTION::ROCKET_KICK_MAGNITUDE_2, { TYPENAME::DOUBLE, "Rocket_Kick_Magnitude(2)", "kms^-1", 14, 6 }}, + { PROGRAM_OPTION::ROCKET_KICK_PHI_1, { TYPENAME::DOUBLE, "Rocket_Kick_Phi(1)", "-", 14, 6 }}, + { PROGRAM_OPTION::ROCKET_KICK_PHI_2, { TYPENAME::DOUBLE, "Rocket_Kick_Phi(2)", "-", 14, 6 }}, + { PROGRAM_OPTION::ROCKET_KICK_THETA_1, { TYPENAME::DOUBLE, "Rocket_Kick_Theta(1)", "-", 14, 6 }}, + { PROGRAM_OPTION::ROCKET_KICK_THETA_2, { TYPENAME::DOUBLE, "Rocket_Kick_Theta(2)", "-", 14, 6 }}, + { PROGRAM_OPTION::ROTATIONAL_VELOCITY_DISTRIBUTION, { TYPENAME::INT, "Rotational_Velocity_Dstrbtn", "-", 4, 1 }}, { PROGRAM_OPTION::ROTATIONAL_FREQUENCY, { TYPENAME::DOUBLE, "Rotational_Frequency", "Hz", 14, 6 }}, { PROGRAM_OPTION::ROTATIONAL_FREQUENCY_1, { TYPENAME::DOUBLE, "Rotational_Frequency(1)", "Hz", 14, 6 }}, diff --git a/src/typedefs.h b/src/typedefs.h index 93b3274bb..22912e055 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -121,6 +121,9 @@ typedef struct SupernovaDetails { // Holds attributes, bool isHydrogenPoor; // Flag to indicate if exploding star is hydrogen-poor. We consider an H-rich star all SN progenitors that have an H envelope, otherwise H-poor double kickMagnitude; // Kick magnitude the system received during the supernova (km s^-1) double kickMagnitudeRandom; // Random number U(0,1) for choosing the supernova kick magnitude - drawn once at star creation + double rocketKickMagnitude; // Rocket kick magnitude the system received after the supernova (km s^-1) + double rocketKickPhi; // Rocket kick azimuthal angle phi the system received after the supernova + double rocketKickTheta; // Rocket kick polar angle theta the system received after the supernova double meanAnomaly; // Mean anomaly at instantaneous time of the SN - uniform in [0, 2pi] double phi; // Kick angle in the orbital plane, defined CCW from the radial vector pointed away from the Companion (rad) [0, 2pi) SN_STATE supernovaState; // Indicates which star (or stars) are undergoing / have undergone a supernova event diff --git a/src/vector3d.cpp b/src/vector3d.cpp index 42721f888..3d07e1eaf 100644 --- a/src/vector3d.cpp +++ b/src/vector3d.cpp @@ -9,7 +9,6 @@ // /////////////////////////////////////////////////////////// -// Default constructor Vector3d::Vector3d() { // Initialize member variables @@ -22,10 +21,10 @@ Vector3d::Vector3d() { } -// Initialize from values - -Vector3d::Vector3d(double p_x, double p_y, double p_z) { +Vector3d::Vector3d(double p_x, double p_y, double p_z) { // Initialize from values + // Initialize member variables + m_ObjectId = globalObjectId++; m_x = p_x; @@ -120,8 +119,9 @@ std::ostream &operator <<(std::ostream &os, Vector3d const p_Vec) { // Add in common vector calculations /* - * Calculate the magnitude of a velocity vector, the speed. - * + * Calculate the magnitude of a velocity vector. + * + * double Magnitude() const * * @return The magnitude of the velocity vector (speed) */ @@ -134,6 +134,7 @@ double Vector3d::Magnitude() const { /* * Convert the Vector3d to a DBL_VECTOR * + * DBL_VECTOR asDBL_VECTOR() * * @return The analogous DBL_VECTOR */ @@ -143,9 +144,66 @@ DBL_VECTOR Vector3d::asDBL_VECTOR() { DBL_VECTOR vFinal = { v[0], v[1], v[2] }; return vFinal; +} + + +////////////////////////////////// +// Add in rotation about an axis + +#define cTheta cos(p_Theta) +#define sTheta sin(p_Theta) +/* + * Rotate a vector about the X axis. + * + * Vector3d RotateVectorAboutX( const double p_Theta) + * + * @param [IN] p_Theta Rotation angle (rad) + * @return Vector after rotation + */ +Vector3d Vector3d::RotateVectorAboutX( const double p_Theta) { + std::vector RotationMatrix = { // Define the Rotation Matrix + { 1.0, 0.0, 0.0 }, + { 0.0, cTheta, -sTheta }, + { 0.0, sTheta, cTheta }}; + return linalg::matrixMult(RotationMatrix, (*this)); // Multiply RotationMatrix * p_oldVector } +/* + * Rotate a vector about the Y axis. + * + * Vector3d RotateVectorAboutY( const double p_Theta) + * + * @param [IN] p_Theta Rotation angle (rad) + * @return Vector after rotation + */ +Vector3d Vector3d::RotateVectorAboutY( const double p_Theta) { + std::vector RotationMatrix = { // Define the Rotation Matrix + { cTheta, 0.0, sTheta }, + { 0.0, 1.0, 0.0 }, + { -sTheta, 0.0, cTheta }}; + return linalg::matrixMult(RotationMatrix, (*this)); // Multiply RotationMatrix * p_oldVector +} + +/* + * Rotate a vector about the Z axis. + * + * Vector3d RotateVectorAboutZ( const double p_Theta) + * + * @param [IN] p_Theta Rotation angle (rad) + * @return Vector after rotation + */ +Vector3d Vector3d::RotateVectorAboutZ( const double p_Theta) { + std::vector RotationMatrix = { // Define the Rotation Matrix + { cTheta, -sTheta, 0.0 }, + { sTheta, cTheta, 0.0 }, + { 0.0, 0.0, 1.0 }}; + return linalg::matrixMult(RotationMatrix, (*this)); // Multiply RotationMatrix * p_oldVector +} + +#undef cTheta +#undef sTheta + /* * Redefines a vector from one coordinate basis to another using Euler Angles @@ -168,13 +226,14 @@ DBL_VECTOR Vector3d::asDBL_VECTOR() { * https://en.wikipedia.org/wiki/Euler_angles * https://en.wikipedia.org/wiki/Change_of_basis * + * Vector3d ChangeBasis(const double p_ThetaE, const double p_PhiE, const double p_PsiE) * * @param [IN] p_ThetaE Euler angle Theta (rad) * @param [IN] p_PhiE Euler angle Phi (rad) * @param [IN] p_PsiE Euler angle Psi (rad) * @return Vector in previous basis */ -Vector3d Vector3d::RotateVector(const double p_ThetaE, const double p_PhiE, const double p_PsiE) { +Vector3d Vector3d::ChangeBasis(const double p_ThetaE, const double p_PhiE, const double p_PsiE) { // Replace for convenience, undefine below #define cTheta cos(p_ThetaE) @@ -213,6 +272,7 @@ Vector3d Vector3d::RotateVector(const double p_ThetaE, const double p_PhiE, cons /* * Returns the unit vector in the direction of the given vector * + * Vector3d UnitVector() * * @return Unit vector of input */ @@ -235,6 +295,7 @@ namespace linalg { /* * Calculate the standard dot product of two vectors * + * double dot(const Vector3d& p_a, const Vector3d& p_b) * * @param [IN] a first vector * @param [IN] b second vector @@ -253,6 +314,7 @@ namespace linalg { /* * Calculate the standard cross product of two vectors * + * Vector3d cross(const Vector3d& p_a, const Vector3d& p_b) * * @param [IN] a first vector * @param [IN] b second vector @@ -272,14 +334,35 @@ namespace linalg { /* * Calculate the angle between two vectors. * + * double angleBetween(const Vector3d& p_a, const Vector3d& p_b) * - * @param [IN] a first vector - * @param [IN] b second vector + * @param [IN] p_a first vector + * @param [IN] p_b second vector * @return angle between them, in radians */ double angleBetween(const Vector3d& p_a, const Vector3d& p_b) { // Angle between 2 vectors, between [0, pi] return std::acos(linalg::dot(p_a, p_b) / (p_a.Magnitude() * p_b.Magnitude())); } + + /* + * Right multiply a matrix by a vector + * + * Vector3d matrixMult(const std::vector& p_matrix, const Vector3d& p_vector) + * + * @param [IN] p_matrix matrix + * @param [IN] p_vector vector + * @return angle between them, in radians + */ + Vector3d matrixMult(const std::vector& p_matrix, const Vector3d& p_vector) { + // Product of matrix and vector + Vector3d newVector = Vector3d(0.0, 0.0, 0.0); + for (size_t i=0; i< 3; i++) { + for (size_t j=0; j<3; j++) { + newVector[i] += p_matrix[i][j] * p_vector[j]; + } + } + return newVector; + } } diff --git a/src/vector3d.h b/src/vector3d.h index 78bb291c3..c083c88d7 100644 --- a/src/vector3d.h +++ b/src/vector3d.h @@ -39,7 +39,10 @@ class Vector3d { DBL_VECTOR asDBL_VECTOR(); // member functions - Vector3d RotateVector( const double p_ThetaE, + Vector3d RotateVectorAboutX( const double p_Theta); + Vector3d RotateVectorAboutY( const double p_Theta); + Vector3d RotateVectorAboutZ( const double p_Theta); + Vector3d ChangeBasis( const double p_ThetaE, const double p_PhiE, const double p_PsiE); Vector3d UnitVector(); @@ -107,6 +110,9 @@ namespace linalg { Vector3d cross(const Vector3d& p_a, const Vector3d& p_b); double angleBetween(const Vector3d& p_a, const Vector3d& p_b); + + Vector3d matrixMult(const std::vector& p_matrix, const Vector3d& p_vector); + } #endif // __vector3d_h__ diff --git a/src/yaml.h b/src/yaml.h index 69fd6ea76..d67cf4151 100644 --- a/src/yaml.h +++ b/src/yaml.h @@ -256,6 +256,12 @@ namespace yaml { " --pisn-upper-limit # Maximum core mass for PISN [Msol]", " --ppi-lower-limit # Minimum core mass for PPI [Msol]", " --ppi-upper-limit # Maximum core mass for PPI [Msol]", + " --rocket-kick-magnitude-1 # (BSE) rocket kick magnitude for the primary should it undergo a supernova event [km/s]", + " --rocket-kick-magnitude-2 # (BSE) rocket kick magnitude for the secondary should it undergo a supernova event [km/s]", + " --rocket-kick-theta-1 # (BSE) angle between the orbital plane and the 'z' axis of the rocket kick for the primary star [radians]", + " --rocket-kick-phi-1 # (BSE) angle between 'x' and 'y', both in the orbital plane of the rocket kick , for the primary star [radians]", + " --rocket-kick-theta-2 # (BSE) angle between the orbital plane and the 'z' axis of the rocket kick for the secondary star [radians]", + " --rocket-kick-phi-2 # (BSE) angle between 'x' and 'y', both in the orbital plane of the rocket kick , for the secondary star [radians]", "", " ### PULSAR PARAMETERS", " --pulsar-birth-magnetic-field-distribution-min # [log10(B/G)]",