Skip to content

Commit

Permalink
Merge pull request #1090 from TeamCOMPAS/issue-1084
Browse files Browse the repository at this point in the history 
Defect Repair: Issue 1084
  • Loading branch information
@ilyamandel
ilyamandel authored Apr 20, 2024
2 parents 5aa2d21 + 4efb884 commit 3fc98f6
Show file tree
Hide file tree
Showing 41 changed files with 816 additions and 781 deletions.
19 changes: 12 additions & 7 deletions src/BH.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,16 +15,22 @@ class BH: virtual public BaseStar, public Remnants {

public:

BH(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), Remnants(p_BaseStar, false) {
if (p_Initialise) Initialise();
BH(const BaseStar &p_BaseStar, const bool p_Initialise = true) : BaseStar(p_BaseStar), Remnants(p_BaseStar) {
m_StellarType = STELLAR_TYPE::BLACK_HOLE; // Set stellar type
if (p_Initialise) Initialise(); // Initialise if required
}

BH& operator = (const BaseStar &p_BaseStar) {
static_cast<BaseStar&>(*this) = p_BaseStar;
Initialise();
return *this;
BH* Clone(const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
BH* clone = new BH(*this, p_Initialise);
clone->SetPersistence(p_Persistence);
return clone;
}

static BH* Clone(BH& p_Star, const OBJECT_PERSISTENCE p_Persistence, const bool p_Initialise = true) {
BH* clone = new BH(p_Star, p_Initialise);
clone->SetPersistence(p_Persistence);
return clone;
}

// member functions - alphabetically
static DBL_DBL_DBL CalculateCoreCollapseSNParams_Static(const double p_Mass);
Expand All @@ -39,7 +45,6 @@ class BH: virtual public BaseStar, public Remnants {
protected:

void Initialise() {
m_StellarType = STELLAR_TYPE::BLACK_HOLE; // Set stellar type
CalculateTimescales(); // Initialise timescales
m_Age = 0.0; // Set age appropriately
}
Expand Down
88 changes: 42 additions & 46 deletions src/BaseBinaryStar.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -233,13 +233,14 @@ void BaseBinaryStar::SetInitialValues(const unsigned long int p_Seed, const long

m_Error = ERROR::NONE;

m_ObjectId = globalObjectId++;
m_ObjectType = OBJECT_TYPE::BASE_BINARY_STAR;
m_StellarType = STELLAR_TYPE::BINARY_STAR;
m_RandomSeed = p_Seed;
m_Id = p_Id;
m_ObjectId = globalObjectId++;
m_ObjectType = OBJECT_TYPE::BASE_BINARY_STAR;
m_ObjectPersistence = OBJECT_PERSISTENCE::PERMANENT;
m_StellarType = STELLAR_TYPE::BINARY_STAR;
m_RandomSeed = p_Seed;
m_Id = p_Id;

m_EvolutionStatus = EVOLUTION_STATUS::CONTINUE;
m_EvolutionStatus = EVOLUTION_STATUS::CONTINUE;

if (OPTIONS->PopulationDataPrinting()) { // user wants to see details of binary?
SAY("Using supplied random seed " << m_RandomSeed << " for Binary Star id = " << m_ObjectId); // yes - show them
Expand Down Expand Up @@ -1161,6 +1162,13 @@ void BaseBinaryStar::ResolveCoalescence() {
* @return True if a supernova event occurred, otherwise false
*/
bool BaseBinaryStar::ResolveSupernova() {
// Functions defined in vector3d.h
// Defined here for convenience - undefined later
#define cross(x,y) Vector3d::Cross(x, y)
#define dot(x,y) Vector3d::Dot(x, y)
#define angleBetween(x,y) Vector3d::AngleBetween(x, y)
#define mag Magnitude()
#define hat UnitVector()

if (!m_Supernova->IsSNevent()) {
SHOW_WARN(ERROR::RESOLVE_SUPERNOVA_IMPROPERLY_CALLED);
Expand Down Expand Up @@ -1201,14 +1209,6 @@ bool BaseBinaryStar::ResolveSupernova() {
else { // no - evaluate orbital changes and calculate velocities
// Evolve SN out of binary

// Functions defined in vector3d.h
// Defined here for convenience - undefined later
#define cross(x,y) linalg::cross(x,y)
#define dot(x,y) linalg::dot(x,y)
#define angleBetween(x,y) linalg::angleBetween(x,y)
#define mag Magnitude()
#define hat UnitVector()

// Pre-SN parameters
double semiMajorAxisPrev_km = m_SemiMajorAxis * AU_TO_KM; // km - Semi-Major axis
double eccentricityPrev = m_Eccentricity; // -- - Eccentricity, written with a prev to distinguish from later use
Expand All @@ -1224,29 +1224,23 @@ bool BaseBinaryStar::ResolveSupernova() {
double sinEccAnomaly = sin(m_Supernova->SN_EccentricAnomaly());

// Derived quantities
double aPrev = semiMajorAxisPrev_km;
double aPrev_2 = aPrev * aPrev;
double aPrev_3 = aPrev_2 * aPrev;

double omega = std::sqrt(G_km_Msol_s * totalMassPrev / aPrev_3); // rad/s - Keplerian orbital frequency
double aPrev = semiMajorAxisPrev_km;
double aPrev_2 = aPrev * aPrev;
double aPrev_3 = aPrev_2 * aPrev;

Vector3d separationVectorPrev = Vector3d(aPrev * (cosEccAnomaly - eccentricityPrev),
aPrev * (sinEccAnomaly) * sqrt1MinusEccPrevSquared,
0.0); // km - Relative position vector, from m1Prev to m2Prev
double separationPrev = separationVectorPrev.mag; // km - Instantaneous Separation
double fact1 = aPrev_2 * omega / separationPrev;
double omega = std::sqrt(G_km_Msol_s * totalMassPrev / aPrev_3); // rad/s - Keplerian orbital frequency

Vector3d relativeVelocityVectorPrev = Vector3d(-fact1 * sinEccAnomaly,
fact1 * cosEccAnomaly * sqrt1MinusEccPrevSquared,
0.0); // km/s - Relative velocity vector, in the m1Prev rest frame
Vector3d separationVectorPrev = Vector3d(aPrev * (cosEccAnomaly - eccentricityPrev), aPrev * (sinEccAnomaly) * sqrt1MinusEccPrevSquared, 0.0); // km - Relative position vector, from m1Prev to m2Prev
double separationPrev = separationVectorPrev.mag; // km - Instantaneous Separation
double fact1 = aPrev_2 * omega / separationPrev;

Vector3d relativeVelocityVectorPrev = Vector3d(-fact1 * sinEccAnomaly, fact1 * cosEccAnomaly * sqrt1MinusEccPrevSquared, 0.0); // km/s - Relative velocity vector, in the m1Prev rest frame
Vector3d orbitalAngularMomentumVectorPrev = cross(separationVectorPrev, relativeVelocityVectorPrev); // km^2 s^-1 - Specific orbital angular momentum vector

Vector3d eccentricityVectorPrev = cross(relativeVelocityVectorPrev, orbitalAngularMomentumVectorPrev) /
(G_km_Msol_s * totalMassPrev) - separationVectorPrev.hat; // -- - Laplace-Runge-Lenz vector (magnitude = eccentricity)

m_OrbitalVelocityPreSN = relativeVelocityVectorPrev.mag; // km/s - Set the Pre-SN orbital velocity and
m_uK = m_Supernova->SN_KickMagnitude() / m_OrbitalVelocityPreSN; // -- - Dimensionless kick magnitude
m_OrbitalVelocityPreSN = relativeVelocityVectorPrev.mag; // km/s - Set the Pre-SN orbital velocity and
m_uK = m_Supernova->SN_KickMagnitude() / m_OrbitalVelocityPreSN; // -- - Dimensionless kick magnitude

// Note: In the following,
// orbitalAngularMomentumVectorPrev defines the Z-axis,
Expand Down Expand Up @@ -1380,19 +1374,15 @@ bool BaseBinaryStar::ResolveSupernova() {
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
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.0 * m_OrbitalEnergy/reducedMass); // AU/yr - average orbital velocity post-SN
double kGrav = 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 / kGrav); // unitless!
double theta_rotation = 3.0 * rocketKickVector.mag * KM_TO_AU * SECONDS_IN_YEAR / (2.0 * averageOrbitalVelocityPreRocket); // rad - need to convert velocities to same units

// Apply hPlus and hMinus support vectors
Vector3d hPlusVector = amVectorNormalizedByCircularAmPreRocket + eccentricityVectorPreRocket;
Expand All @@ -1403,8 +1393,8 @@ bool BaseBinaryStar::ResolveSupernova() {
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 );
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);
Expand All @@ -1415,7 +1405,7 @@ bool BaseBinaryStar::ResolveSupernova() {
Vector3d eccentricityVectorPostRocket = 0.5 * (hPlusVector_prime - hMinusVector_prime);

m_NormalizedOrbitalAngularMomentumVector = normalizedAngularMomentumVectorPostRocket ;
m_Eccentricity = eccentricityVectorPostRocket.mag;
m_Eccentricity = eccentricityVectorPostRocket.mag;

UpdateSystemicVelocity(rocketKickVector.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
m_Supernova->UpdateComponentVelocity(rocketKickVector.ChangeBasis(m_ThetaE, m_PhiE, m_PsiE));
Expand All @@ -1440,6 +1430,12 @@ bool BaseBinaryStar::ResolveSupernova() {
m_Supernova->ClearCurrentSNEvent();

return true;

#undef hat
#undef mag
#undef angleBetween
#undef dot
#undef cross
}


Expand Down
Loading

0 comments on commit 3fc98f6

Please sign in to comment.