Skip to content

Fix output NaNs and event names #93

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
astroJeff wants to merge 3 commits into
base: development-deprecated
Choose a base branch
from
Draft

Fix output NaNs and event names #93

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
d8bc6c5
Updating strings for _start and _end for CE and CC, removing NaN's fr…
astroJeff Apr 19, 2023
06b54e7
Completing merge
astroJeff Jun 8, 2023
d9f5982
Adding more _start and _end consistency across states
astroJeff Jun 8, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
112 changes: 90 additions & 22 deletions posydon/binary_evol/CE/step_CEE.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -186,18 +186,18 @@ def __init__(
def __call__(self, binary):
"""Perform the CEE step for a BinaryStar object."""
# Determine which star is the donor and which is the companion
if binary.event in ["oCE1", "oDoubleCE1"]:
if binary.event in ["CE1_start", "DCE_start"]:
donor_star = binary.star_1
comp_star = binary.star_2
elif binary.event in ["oCE2", "oDoubleCE2"]:
elif binary.event in ["CE2_start", "DCE_start"]:
donor_star = binary.star_2
comp_star = binary.star_1
else:
raise ValueError("CEE does not apply if `event` is not "
"`oCE1`, 'oDoubleCE1' or `oCE2`, 'oDoubleCE1'")
"'CE1_start', 'CE2_start', or 'DCE_start'")

# Check for double CE
if binary.event in ["oDoubleCE1", "oDoubleCE2"]:
if binary.event is "DCE_start":
double_CE = True
else:
double_CE = False
Expand All @@ -214,23 +214,61 @@ def __call__(self, binary):
'stripped_He_Core_He_burning']:
# system merges
binary.state = 'merged'
if binary.event in ["oCE1", "oDoubleCE1"]:
binary.event = "oMerging1"
if binary.event in ["oCE2", "oDoubleCE2"]:
binary.event = "oMerging2"
if double_CE:
binary.event = 'DCE_merger'
Copy link
Copy Markdown
Contributor

@ezapartas ezapartas Jun 15, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe make it obvious it is at the onset/start of the process of merging (not just in general a merger product), by making it "DCE_merger_start", consistent with your implementation of "*_start" in this PR

else:
binary.event = 'CE_merger'
Copy link
Copy Markdown
Contributor

@ezapartas ezapartas Jun 15, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Same here, CE_merger_start (or CE_merging_start)


for key in BINARYPROPERTIES:
# the binary attributes that are changed in the CE step
if key not in ["time", "state", "event", 'V_sys',
'mass_transfer_case',
'nearest_neighbour_distance']:
setattr(binary, key, np.nan) # the rest become np.nan
if key == 'mass_transfer_case':
setattr(binary, key, 'None')
stars = [donor_star, comp_star]
# for now we just add all initial mass to the (merger) star_1
masses = [donor_star.mass + comp_star.mass, np.nan]
star_states = [donor_star.state, comp_star.state]
for star, star_state, mass in zip(stars, star_states, masses):
star.mass = mass
star.state = star_state
for key in STARPROPERTIES:
# the binary attributes that are changed in the CE step
if key not in ["mass", "state"]:
setattr(star, key, np.nan)
return

if self.common_envelope_option_for_HG_star == "pessimistic":
# Merging if HG donor
if donor_star.state in ['H-rich_Shell_H_burning']:
# system merges
binary.state = 'merged'
if binary.event in ["oCE1", "oDoubleCE1"]:
binary.event = "oMerging1"
if binary.event in ["oCE2", "oDoubleCE2"]:
binary.event = "oMerging2"

if double_CE:
binary.event = 'DCE_merger'
else:
binary.event = 'CE_merger'

for key in BINARYPROPERTIES:
# the binary attributes that are changed in the CE step
if key not in ["time", "state", "event", 'V_sys',
'mass_transfer_case',
'nearest_neighbour_distance']:
setattr(binary, key, np.nan) # the rest become np.nan
if key == 'mass_transfer_case':
setattr(binary, key, 'None')
stars = [donor_star, comp_star]
# for now we just add all initial mass to the (merger) star_1
masses = [donor_star.mass + comp_star.mass, np.nan]
star_states = [donor_star.state, comp_star.state]
for star, star_state, mass in zip(stars, star_states, masses):
star.mass = mass
star.state = star_state
for key in STARPROPERTIES:
# the binary attributes that are changed in the CE step
if key not in ["mass", "state"]:
setattr(star, key, np.nan)
return

# Calculate binary's evolution
Expand Down Expand Up @@ -668,26 +706,38 @@ def CEE_simple_alpha_prescription(
print("during the assumed windloss phase after postCE")
print("with 'common_envelope_option_after_succ_CEE' :",
common_envelope_option_after_succ_CEE)
print("the orbit cahnged from postCEE : ",
print("the orbit changed from postCEE : ",
orbital_period_postCEE)
print("to : ", orbital_period_f)
separation_f = cf.orbital_separation_from_period(orbital_period_f,
mc1_f, mc2_f)

if state1_i == 'stripped_He_Central_He_depleted':
if donor == binary.star_1:
binary.event = 'CC1'
binary.event = 'CC1_start'
elif donor == binary.star_2:
binary.event = 'CC2'
binary.event = 'CC2_start'
else:
binary.event = None
if binary.event == 'CE1_start':
binary.event = 'CE1_end'
elif binary.event == 'CE2_start':
binary.event = 'CE2_end'
elif double_CE:
binary.event = 'DCE_end'
else:
raise ValueError(
"Impossible binary.event condition. binary.event = {}"
"dont know how to proceed",
format(binary.event))




# Adjust binary properties
binary.separation = separation_f
binary.orbital_period = orbital_period_f
binary.eccentricity = 0.0
binary.state = 'detached'
# binary.event = None

for key in BINARYPROPERTIES:
# the binary attributes that are changed in the CE step
Expand Down Expand Up @@ -803,10 +853,28 @@ def CEE_simple_alpha_prescription(
else:
# system merges
binary.state = 'merged'
if binary.event in ["oCE1", "oDoubleCE1"]:
binary.event = "oMerging1"
if binary.event in ["oCE2", "oDoubleCE2"]:
binary.event = "oMerging2"
binary.event = 'CE_merger'

for key in BINARYPROPERTIES:
# the binary attributes that are changed in the CE step
if key not in ["time", "state", "event", 'V_sys',
'mass_transfer_case',
'nearest_neighbour_distance']:
setattr(binary, key, np.nan) # the rest become np.nan
if key == 'mass_transfer_case':
setattr(binary, key, 'None')
# binary.event = 'END'
stars = [donor, comp_star]
# for now we just add all initial mass to the (merger) star_1
masses = [m1_i+m2_i, np.nan]
star_states = [donor.state, comp_star.state]
for star, star_state, mass in zip(stars, star_states, masses):
star.mass = mass
star.state = star_state
for key in STARPROPERTIES:
# the binary attributes that are changed in the CE step
if key not in ["mass", "state"]:
setattr(star, key, np.nan)

if verbose:
print("system merges due to one of the two star's core filling"
Expand Down
46 changes: 26 additions & 20 deletions posydon/binary_evol/DT/step_detached.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1017,7 +1017,7 @@ def __call__(self, binary):
raise Exception("State not recognized!")
else:
raise Exception("Non existent companion has not a recognized value!")

def get_star_data(binary, star1, star2, htrack,
co, copy_prev_m0=None, copy_prev_t0=None):
"""Get and interpolate the properties of stars.
Expand Down Expand Up @@ -1128,10 +1128,10 @@ def get_star_data(binary, star1, star2, htrack,
# get the matched data of two stars, respectively
interp1d_sec, m0, t0 = get_star_data(
binary, secondary, primary, secondary.htrack, co=False)

primary_not_normal = (primary.co) or (self.non_existent_companion in [1,2])
primary_normal = (not primary.co) and self.non_existent_companion == 0
primary_normal = (not primary.co) and self.non_existent_companion == 0

if primary_not_normal:
# copy the secondary star except mass which is of the primary,
# and radius, mdot, Idot = 0
Expand All @@ -1143,8 +1143,8 @@ def get_star_data(binary, star1, star2, htrack,
binary, primary, secondary, primary.htrack, False)[0]
else:
raise Exception("During matching primary is either should be either normal or not normal. `non_existent_companion` should be zero.")


if interp1d_sec is None or interp1d_pri is None:
# binary.event = "END"
binary.state += " (GridMatchingFailed)"
Expand Down Expand Up @@ -1844,25 +1844,25 @@ def get_star_profile(star, htrack, m0):
if s.t_events[0]:
if secondary == binary.star_1:
binary.state = "RLO1"
binary.event = "oRLO1"
binary.event = "RLO1_start"
else:
binary.state = "RLO2"
binary.event = "oRLO2"
binary.event = "RLO2_start"
elif s.t_events[1]:
if secondary == binary.star_1:
binary.state = "RLO2"
binary.event = "oRLO2"
binary.event = "RLO2_start"
else:
binary.state = "RLO1"
binary.event = "oRLO1"
binary.event = "RLO1_start"

elif s.t_events[2]:
# reached t_max of track. End of life (possible collapse) of
# secondary
if secondary == binary.star_1:
binary.event = "CC1"
binary.event = "CC1_start"
else:
binary.event = "CC2"
binary.event = "CC2_start"
get_star_final_values(secondary, secondary.htrack, m01)
get_star_profile(secondary, secondary.htrack, m01)
if not primary.co and primary.state in STAR_STATES_CC:
Expand All @@ -1883,9 +1883,9 @@ def get_star_profile(star, htrack, m0):
# reached t_max of track. End of life (possible collapse) of
# primary
if secondary == binary.star_1:
binary.event = "CC2"
binary.event = "CC2_start"
else:
binary.event = "CC1"
binary.event = "CC1_start"
get_star_final_values(primary, primary.htrack, m02)
get_star_profile(primary, primary.htrack, m02)
else: # Reached max_time asked.
Expand Down Expand Up @@ -1987,6 +1987,9 @@ def diffeq(
Convective turnover time of the star, calculated @
0.5*pressure_scale_height above the bottom of the outer convection
zone in yr.
wdivwc: float
The ratio of angular rotation rate (omega) to critical angular
rotation rate (omega_crit). This is dimensionless.
L : float
Luminosity of the star in solar units.
#mass_conv_core : float
Expand Down Expand Up @@ -2401,10 +2404,11 @@ def diffeq(

if magnetic_braking_mode == "RVJ83":
# Torque from Rappaport, Verbunt, and Joss 1983, ApJ, 275, 713
# The torque is eq.36 of Rapport+1983, with γ = 4. Torque units
# converted from cgs units to [Msol], [Rsol], [yr] as all stellar
# parameters are given in units of [Msol], [Rsol], [yr] and so that
# dOmega_mb/dt is in units of [yr^-2].
# The torque is eq.36 of Rapport+1983, with γ = 4
# Torque units converted from cgs units to [Msol], [Rsol], [yr]
# as all stellar parameters are given in units of [Msol], [Rsol], [yr]
# and so that dOmega_mb/dt is in units of [yr^-2].

dOmega_mb_sec = (
-3.8e30 * (const.rsol**2 / const.secyer)
* M_sec
Expand Down Expand Up @@ -2515,8 +2519,10 @@ def diffeq(
Prot_sec = 2 * np.pi / Omega_sec # [yr]
Rossby_number_sec = Prot_sec / tau_conv_sec

n_pri = (0.03 / Rossby_number_pri) + 0.5 * Rossby_number_pri + 1.0
n_sec = (0.03 / Rossby_number_sec) + 0.5 * Rossby_number_sec + 1.0
n_pri = (0.03 / Rossby_number_pri) \
+ 0.5 * Rossby_number_pri + 1.0
n_sec = (0.03 / Rossby_number_sec) \
+ 0.5 * Rossby_number_sec + 1.0

Qn_pri = 4.05 * np.exp(-1.4 * n_pri)
Qn_sec = 4.05 * np.exp(-1.4 * n_sec)
Expand Down
18 changes: 9 additions & 9 deletions posydon/binary_evol/MESA/step_mesa.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1282,10 +1282,10 @@ def __call__(self, binary):
self.binary.event = 'redirect'
return
elif (state_1 == 'H-rich_Central_C_depletion'): # redirect if CC1
self.binary.event = 'CC1'
self.binary.event = 'CC1_start'
return
elif (state_2 == 'H-rich_Central_C_depletion'): # redirect if CC2
self.binary.event = 'CC2'
self.binary.event = 'CC2_start'
return
else:
raise ValueError('The star_1.state = %s, star_2.state = %s, '
Expand Down Expand Up @@ -1333,29 +1333,29 @@ def __call__(self, binary):

# check the star states
if (state_2 in ["WD", "NS", "BH"]
and (state_1 in FOR_RLO_STATES) and event == "oRLO1"):
and (state_1 in FOR_RLO_STATES) and event == "RLO1_start"):
self.flip_stars_before_step = False
m1 = self.binary.star_1.mass
m2 = self.binary.star_2.mass
# catch and redirect double core collapse, this happens if q=1:
if self.binary.star_1.state == 'H-rich_Central_C_depletion':
self.binary.event = 'CC1'
self.binary.event = 'CC1_start'
return
# super().__call__(binary)
elif (state_1 in ["WD", "NS", "BH"] and (state_2 in FOR_RLO_STATES)
and event == "oRLO2"):
and event == "RLO2_start"):
self.flip_stars_before_step = True
m1 = self.binary.star_2.mass
m2 = self.binary.star_1.mass
# catch and redirect double core collapse, this happens if q=1:
if self.binary.star_2.state == 'H-rich_Central_C_depletion':
self.binary.event = 'CC2'
self.binary.event = 'CC2_start'
return
# super().__call__(binary)
else:
raise ValueError(
'The star_1.state = %s, star_2.state = %s, binary.state = %s, '
'binary.event = %s and not CO - HMS - oRLO1/oRLO2!'
'binary.event = %s and not CO - HMS - RLO1_start/RLO2_start!'
% (state_1, state_2, state, event))
# redirect if outside grids
if 0.466 <= m1 <= 128.735 and 0.092 <= m2 <= 39.25 and p <= 3780.83:
Expand Down Expand Up @@ -1411,7 +1411,7 @@ def __call__(self, binary):
m2 = self.binary.star_1.mass
# catch and redirect double core collapse, this happens if q=1:
if self.binary.star_1.state == 'stripped_He_Central_C_depletion':
self.binary.event = 'CC1'
self.binary.event = 'CC1_start'
# REMOVED assume circularisation after first CC
# new_separation = self.binary.separation*(
# 1.-self.binary.eccentricity**2)
Expand All @@ -1427,7 +1427,7 @@ def __call__(self, binary):
m2 = self.binary.star_2.mass
# catch and redirect double core collapse, this happens if q=1:
if self.binary.star_2.state == 'stripped_He_Central_C_depletion':
self.binary.event = 'CC2'
self.binary.event = 'CC2_start'
return
else:
raise ValueError(
Expand Down
Loading