Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/check_if_latex_present' into add…
Browse files Browse the repository at this point in the history 
…_plot_in_pr_after_successful_workflow
  • Loading branch information
@avivajpeyi
avivajpeyi committed Apr 9, 2024
2 parents 10f12e0 + ce1597e commit 0a64173
Show file tree
Hide file tree
Showing 22 changed files with 407 additions and 50 deletions.
16 changes: 10 additions & 6 deletions .github/workflows/compas-compile-ci.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,12 +30,16 @@ jobs:
python-version: '3.9'
cache: 'pip' # caching pip dependencies

- name: Install dependencies on ubuntu
if: startsWith(matrix.os, 'ubuntu-20')
run: |
cd misc/cicd-scripts
chmod 755 linux-dependencies
./linux-dependencies
- name: Install TeXLive
uses: teatimeguest/setup-texlive-action@v3


- name: Install dependencies on ubuntu
if: startsWith(matrix.os, 'ubuntu-20')
run: |
cd misc/cicd-scripts
chmod 755 linux-dependencies
./linux-dependencies
- name: Build Compas
run: |
Expand Down
3 changes: 2 additions & 1 deletion compas_python_utils/detailed_evolution_plotter/plot_detailed_evolution.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,6 +5,7 @@
###################################################################

import os, sys
import shutil
import numpy as np
import h5py as h5
import matplotlib.pyplot as plt
Expand Down Expand Up @@ -45,7 +46,7 @@ def run_main_plotter(data_path, outdir='.', show=True):

fontparams = {
"font.serif": "Times New Roman",
"text.usetex": "True",
"text.usetex": str(shutil.which("latex") is not None),
"axes.grid": "True",
"grid.color": "gray",
"grid.linestyle": ":",
Expand Down
16 changes: 16 additions & 0 deletions ...ges/User guide/COMPAS output/standard-logfiles-record-specification-stellar.rst
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2050,6 +2050,22 @@ but not both. If both are printed then the file will contain two columns with th
* - Header Strings:
- Teff<CE(1), Teff<CE(2), Teff<CE(SN), Teff<CE(CP)

.. flat-table::
:widths: 25 75 1 1
:header-rows: 0
:class: aligned-text

* - :cspan:`2` **TZAMS**
-
* - Data type:
- DOUBLE
* - COMPAS variable:
- BaseStar::m_TZAMS
* - Description:
- ZAMS Effective Temperature (K).
* - Header Strings:
- Teff@\ ZAMS, Teff@ZAMS(1), Teff@ZAMS(2), Teff@ZAMS(SN), Teff@ZAMS(CP)

.. flat-table::
:widths: 25 75 1 1
:header-rows: 0
Expand Down
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
1 change: 0 additions & 1 deletion py_tests/test_plot_detailed_evolution.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,7 +5,6 @@
from compas_python_utils.detailed_evolution_plotter import plot_detailed_evolution


@pytest.mark.skip(reason="RuntimeError: Failed to process string with tex because latex could not be found")
def test_plotter(example_compas_output_path, capsys, test_archive_dir):
data_path = example_compas_output_path
bse_detailed_out_path = os.path.join(
Expand Down
95 changes: 75 additions & 20 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 Expand Up @@ -1879,7 +1934,7 @@ void BaseBinaryStar::CalculateMassTransfer(const double p_Dt) {
bool accretorIsWD = m_Accretor->IsOneOf(WHITE_DWARFS);

// Determine stability
bool isUnstable;
bool isUnstable = false;
if (donorIsHeHGorHeGB && (caseBBAlwaysStable || caseBBAlwaysUnstable || (caseBBAlwaysUnstableOntoNSBH && accretorIsNSorBH))) { // Determine stability based on case BB
isUnstable = (caseBBAlwaysUnstable || (caseBBAlwaysUnstableOntoNSBH && accretorIsNSorBH)); // Already established that donor is HeHG or HeGB - need to check if new case BB prescriptions are added
}
Expand Down
Loading

0 comments on commit 0a64173

Please sign in to comment.