Skip to content

Commit

Permalink
Merge pull request #1089 from jeffriley/issue-1084
Browse files Browse the repository at this point in the history 
Defect repair: Issue 1084
  • Loading branch information
@jeffriley
jeffriley authored Apr 14, 2024
2 parents 1c1c42c + 3f0af9f commit 5c962ac
Show file tree
Hide file tree
Showing 41 changed files with 610 additions and 680 deletions.
6 changes: 0 additions & 6 deletions src/BH.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,12 +19,6 @@ class BH: virtual public BaseStar, public Remnants {
if (p_Initialise) Initialise();
}

BH& operator = (const BaseStar &p_BaseStar) {
static_cast<BaseStar&>(*this) = p_BaseStar;
Initialise();
return *this;
}


// member functions - alphabetically
static DBL_DBL_DBL CalculateCoreCollapseSNParams_Static(const double p_Mass);
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
34 changes: 21 additions & 13 deletions src/BaseBinaryStar.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -142,9 +142,10 @@ class BaseBinaryStar {
// Copy constructor
BaseBinaryStar(const BaseBinaryStar& p_Star) {

m_ObjectId = globalObjectId++; // get unique object id (don't copy source)
m_ObjectType = OBJECT_TYPE::BASE_BINARY_STAR; // can only copy from BASE_BINARY_STAR
m_StellarType = STELLAR_TYPE::BINARY_STAR; // always
m_ObjectId = globalObjectId++; // get unique object id (don't copy source)
m_ObjectType = OBJECT_TYPE::BASE_BINARY_STAR; // can only copy from BASE_BINARY_STAR
m_ObjectPersistence = OBJECT_PERSISTENCE::PERMANENT; // permanent - not an ephemeral clone
m_StellarType = STELLAR_TYPE::BINARY_STAR; // always

CopyMemberVariables(p_Star); // copy member variables
}
Expand All @@ -155,9 +156,10 @@ class BaseBinaryStar {

if (this != &p_Star) { // make sure we're not not copying ourselves...

m_ObjectId = globalObjectId++; // get unique object id (don't copy source)
m_ObjectType = OBJECT_TYPE::BASE_BINARY_STAR; // can only copy from BASE_BINARY_STAR
m_StellarType = STELLAR_TYPE::BINARY_STAR; // always
m_ObjectId = globalObjectId++; // get unique object id (don't copy source)
m_ObjectType = OBJECT_TYPE::BASE_BINARY_STAR; // can only copy from BASE_BINARY_STAR
m_ObjectPersistence = OBJECT_PERSISTENCE::PERMANENT; // permanent - not an ephemeral clone
m_StellarType = STELLAR_TYPE::BINARY_STAR; // always

CopyMemberVariables(p_Star); // copy member variables
}
Expand All @@ -171,6 +173,7 @@ class BaseBinaryStar {
// object identifiers - all classes have these
OBJECT_ID ObjectId() const { return m_ObjectId; }
OBJECT_TYPE ObjectType() const { return m_ObjectType; }
OBJECT_PERSISTENCE ObjectPersistence() const { return m_ObjectPersistence; }
STELLAR_TYPE StellarType() const { return m_StellarType; }
long int Id() const { return m_Id; }

Expand Down Expand Up @@ -271,6 +274,9 @@ class BaseBinaryStar {
double ZetaLobe() const { return m_ZetaLobe; }
double ZetaStar() const { return m_ZetaStar; }

// setters
void SetObjectId(const OBJECT_ID p_ObjectId) { m_ObjectId = p_ObjectId; }
void SetPersistence(const OBJECT_PERSISTENCE p_Persistence) { m_ObjectPersistence = p_Persistence; }

// member functions - alphabetically
COMPAS_VARIABLE BinaryPropertyValue(const T_ANY_PROPERTY p_Property) const;
Expand All @@ -279,7 +285,7 @@ class BaseBinaryStar {

EVOLUTION_STATUS Evolve();

bool PrintSwitchLog(const bool p_PrimarySwitching) { return OPTIONS->SwitchLog() ? LOGGING->LogBSESwitchLog(this, p_PrimarySwitching) : true; }
bool PrintSwitchLog(const bool p_PrimarySwitching) { return OPTIONS->SwitchLog() ? (LOGGING->ObjectSwitchingPersistence() == OBJECT_PERSISTENCE::PERMANENT ? LOGGING->LogBSESwitchLog(this, p_PrimarySwitching) : true) : true; }

COMPAS_VARIABLE PropertyValue(const T_ANY_PROPERTY p_Property) const;

Expand All @@ -290,12 +296,13 @@ class BaseBinaryStar {

BaseBinaryStar() { }

OBJECT_ID m_ObjectId; // Instantiated object's unique object id
OBJECT_TYPE m_ObjectType; // Instantiated object's object type
STELLAR_TYPE m_StellarType; // Stellar type defined in Hurley et al. 2000
long int m_Id; // Id used to name detailed output file - uses p_Id as passed (usually the index number of multiple binaries being produced)
OBJECT_ID m_ObjectId; // Instantiated object's unique object id
OBJECT_TYPE m_ObjectType; // Instantiated object's object type
OBJECT_PERSISTENCE m_ObjectPersistence; // Instantiated object's persistence (permanent or ephemeral)
STELLAR_TYPE m_StellarType; // Stellar type defined in Hurley et al. 2000
long int m_Id; // Id used to name detailed output file - uses p_Id as passed (usually the index number of multiple binaries being produced)

ERROR m_Error; // Records most recent error encountered for this binary
ERROR m_Error; // Records most recent error encountered for this binary

// member variables - alphabetical in groups (sort of...)

Expand Down Expand Up @@ -559,7 +566,8 @@ class BaseBinaryStar {
double donorMass = m_Donor->Mass();
double accretorMass = m_Accretor->Mass();

BinaryConstituentStar* donorCopy = new BinaryConstituentStar(*m_Donor);
BinaryConstituentStar *donorCopy = BinaryConstituentStar::Clone(*m_Donor, OBJECT_PERSISTENCE::EPHEMERAL);

double semiMajorAxis = m_Binary->CalculateMassTransferOrbit(donorCopy->Mass(), -p_dM , *m_Accretor, m_FractionAccreted);
double RLRadius = semiMajorAxis * (1.0 - m_Binary->Eccentricity()) * CalculateRocheLobeRadius_Static(donorMass - p_dM, accretorMass + (m_Binary->FractionAccreted() * p_dM)) * AU_TO_RSOL;

Expand Down
Loading

0 comments on commit 5c962ac

Please sign in to comment.