diff --git a/posydon/binary_evol/DT/step_merged.py b/posydon/binary_evol/DT/step_merged.py index 9f48673688..0bc265a3fb 100644 --- a/posydon/binary_evol/DT/step_merged.py +++ b/posydon/binary_evol/DT/step_merged.py @@ -23,6 +23,7 @@ ) from posydon.config import PATH_TO_POSYDON_DATA from posydon.utils.common_functions import check_state_of_star +from posydon.utils.constants import Zsun, zams_table from posydon.utils.posydonerror import ModelError from posydon.utils.posydonwarning import Pwarn @@ -47,10 +48,15 @@ def __init__( grid_name_strippedHe=None, path=PATH_TO_POSYDON_DATA, merger_critical_rot = 0.4, - rel_mass_lost_HMS_HMS = 0.1, + rel_mass_lost_HMS_HMS = "Glebbeek+2013", # [0-1) or "Glebbeek+2013" (dependent on mass ratio q; Glebbeek E., Gaburov E., Portegies Zwart S., Pols O. R., 2013, MNRAS, 434, 3497) + HMS_HMS_merging_rejuvenation = True, # if True then new abundances based on total_mass_h1 or he4 + # ("Schneider+2016"; i.e. Schneider, F. R. N., Podsiadlowski, P., Langer, N., Castro, N., & Fossati, L. 2016, MNRAS,457, 2355 + # https://ui.adsabs.harvard.edu/abs/2016MNRAS.457.2355S/abstract) list_for_matching_HMS = [ - ["mass", "center_h1", "he_core_mass"], - [20.0, 1.0, 10.0], + #["mass", "center_h1", "he_core_mass"], #mostly suggested for HMS_HMS_merging_rejuvenation = False + #[20.0, 1.0, 10.0], + ["mass", "total_mass_h1", "he_core_mass"], #mostly suggested for HMS_HMS_merging_rejuvenation = True + [20.0, 20.0, 10.0], ["log_min_max", "min_max", "min_max"], #[m_min_H, m_max_H], [0, None] [None, None], [0, None] @@ -75,6 +81,7 @@ def __init__( self.merger_critical_rot = merger_critical_rot self.rel_mass_lost_HMS_HMS = rel_mass_lost_HMS_HMS + self.HMS_HMS_merging_rejuvenation = HMS_HMS_merging_rejuvenation super().__init__( grid_name_Hrich=grid_name_Hrich, @@ -93,12 +100,13 @@ def __call__(self,binary): if self.verbose: print("Before Merger", binary.star_1.state, binary.star_2.state, binary.state, binary.event) - print("M1 , M2, he_core_mass1, he_core_mass2: ", + print("M1 , M2, he_core_mass1, he_core_mass2, co_core_mass1, co_core_mass2: ", binary.star_1.mass, binary.star_2.mass, - binary.star_1.he_core_mass, binary.star_2.he_core_mass) - print("star_1.center_he4, star_2.center_he4, star_1.surface_he4, " - "star_2.surface_he4: ", binary.star_1.center_he4, - binary.star_2.center_he4, binary.star_1.surface_he4, + binary.star_1.he_core_mass, binary.star_2.he_core_mass, binary.star_1.co_core_mass, binary.star_2.co_core_mass) + print("star_1.center_h1, star_2.center_h1, star_1.center_he4, star_1.center_c12, star_2.center_c12, " + "star_2.center_he4, star_1.surface_he4, star_2.surface_he4: ", + binary.star_1.center_h1, binary.star_2.center_h1, binary.star_1.center_he4, + binary.star_2.center_he4, binary.star_1.center_c12, binary.star_2.center_c12, binary.star_1.surface_he4, binary.star_2.surface_he4) if binary.state == "merged": @@ -126,10 +134,16 @@ def __call__(self,binary): if self.verbose: print("After Merger", binary.star_1.state, binary.star_2.state, binary.state, binary.event) - print("M_merged , he_core_mass merged: ", binary.star_1.mass, - binary.star_1.he_core_mass) - print("star_1.center_he4, star_1.surface_he4: ", - binary.star_1.center_he4, binary.star_1.surface_he4) + if binary.event == 'oMerging1': + print("M_merged , he_core_mass merged, co_core_mass merged: ", binary.star_1.mass, + binary.star_1.he_core_mass, binary.star_1.co_core_mass) + print("star_1.center_h1, star_1.center_he4, star_1.center_c12, star_1.surface_he4: ", + binary.star_1.center_h1, binary.star_1.center_he4, binary.star_1.center_c12, binary.star_1.surface_he4) + elif binary.event=='oMerging2': + print("M_merged , he_core_mass merged, co_core_mass merged: ", binary.star_2.mass, + binary.star_2.he_core_mass, binary.star_2.co_core_mass) + print("star_2.center_h1, star_2.center_he4, star_2.center_c12, star_2.surface_he4: ", + binary.star_2.center_h1, binary.star_2.center_he4, binary.star_2.center_c12, binary.star_2.surface_he4) super().__call__(binary) @@ -150,6 +164,7 @@ def merged_star_properties(self, star_base, comp): s1 = star_base.state s2 = comp.state + def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", mass_weight1="mass", mass_weight2=None): A1 = getattr(star1, abundance_name) A2 = getattr(star2, abundance_name) @@ -168,27 +183,88 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma M2 = getattr(star2, "he_core_mass") - getattr(star2, "co_core_mass") else: M2 = getattr(star2, mass_weight2) - return (A1*M1 + A2*M2 ) / (M1+M2) + + try: + mass_weighted_avg_value=(A1*M1+A2*M2)/(M1+M2) + except TypeError: + mass_weighted_avg_value= np.nan + + if self.verbose: + print(abundance_name, mass_weight1, mass_weight2) + print("A_base, M_base_abundance, A_comp, M_comp_abundanc", A1, M1, A2, M2) + print("mass_weighted_avg= ", mass_weighted_avg_value) + + return mass_weighted_avg_value # MS + MS if ( s1 in LIST_ACCEPTABLE_STATES_FOR_HMS and s2 in LIST_ACCEPTABLE_STATES_FOR_HMS): - #these stellar attributes change value - merged_star.mass = (star_base.mass + comp.mass) * (1.-self.rel_mass_lost_HMS_HMS) - #TODO for key in ["center_h1", "center_he4", "center_c12", "center_n14","center_o16"]: - merged_star.center_h1 = mass_weighted_avg() - merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4") - merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12") - merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14") - merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16") - #TODO: should I check if the abundaces above end up in ~1 (?) - - # weigheted mixing on the surface abundances - merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1") - merged_star.surface_he4 = mass_weighted_avg(abundance_name = "surface_he4") - merged_star.surface_c12 = mass_weighted_avg(abundance_name = "surface_c12") - merged_star.surface_n14 = mass_weighted_avg(abundance_name = "surface_n14") - merged_star.surface_o16 = mass_weighted_avg(abundance_name = "surface_o16") + + q = comp.mass/star_base.mass + # first we calculate the mass loss percentage: + if self.rel_mass_lost_HMS_HMS == "Glebbeek+2013": + # Eq. 4 of Glebbeek+2013 + rel_mass_lost_HMS_HMS = (0.3 * q) / ((1.0 + q) ** 2.0) + elif isinstance(self.rel_mass_lost_HMS_HMS, (int, float)) and 0.0 <= self.rel_mass_lost_HMS_HMS < 1.0: + rel_mass_lost_HMS_HMS = self.rel_mass_lost_HMS_HMS + else: + raise ValueError( + 'rel_mass_lost_HMS_HMS must be either a float in [0, 1) ' + 'or the string "Glebbeek+2013".' + ) + + + if self.HMS_HMS_merging_rejuvenation : # according to Schneider+2016 + X_average1 = star_base.total_mass_h1 / star_base.mass + X_average2 = comp.total_mass_h1 / comp.mass + + Zinitial_div_Zsun = star_base.metallicity + if Zinitial_div_Zsun in zams_table.keys(): + Y_initial = zams_table[Zinitial_div_Zsun] + else: + raise KeyError(f"{Zinitial_div_Zsun} is a not defined metallicity") + Z_initial = Zinitial_div_Zsun*Zsun + X_initial = 1.0 - Y_initial - Z_initial + + merged_star.mass = (star_base.mass + comp.mass) * (1.-rel_mass_lost_HMS_HMS) + + X_average_merged = (star_base.total_mass_h1 + comp.total_mass_h1 - (1.-rel_mass_lost_HMS_HMS) * X_initial) / merged_star.mass + Y_average_merged = (star_base.total_mass_he4 + comp.total_mass_he4 - (1.-rel_mass_lost_HMS_HMS) * Y_initial) / merged_star.mass + + merged_star.total_mass_h1 = X_average_merged * merged_star.mass + merged_star.center_h1 = X_average_merged # this is only consistent with the merged_star.total_mass_h1 above if we assume full mixing during merger + merged_star.surface_h1 = X_average_merged # this is only consistent with the merged_star.total_mass_h1 above if we assume full mixing during merger + #in any case, the track_match after the step_merged will decide which parameter will be used for matching (center_h1, center_he4 or total_mass_h1) + + merged_star.total_mass_he4 = Y_average_merged * merged_star.mass + merged_star.center_he4 = Y_average_merged # this is only consistent with the merged_star.total_mass_he4 above if we assume full mixing during merger + merged_star.surface_he4 = Y_average_merged # this is only consistent with the merged_star.total_mass_he4 above if we assume full mixing during merger + #in any case, the track_match after the step_merged will decide which parameter will be used for matching (center_h1, center_he4 or total_mass_h1) + + for key in STARPROPERTIES: + # these stellar attributes become np.nan as + if key in [ "center_c12", "center_n14", "center_o16"]: + setattr(merged_star, key, np.nan) + + else: + merged_star.center_h1 = mass_weighted_avg() + merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4") + merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12") + merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14") + merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16") + #TODO: should I check if the abundaces above end up in ~1 (?) + + # weigheted mixing on the surface abundances + merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1") + merged_star.surface_he4 = mass_weighted_avg(abundance_name = "surface_he4") + merged_star.surface_c12 = mass_weighted_avg(abundance_name = "surface_c12") + merged_star.surface_n14 = mass_weighted_avg(abundance_name = "surface_n14") + merged_star.surface_o16 = mass_weighted_avg(abundance_name = "surface_o16") + + # The change of masses occurs after the calculation of weighted averages + merged_star.mass = (star_base.mass + comp.mass) * (1.-rel_mass_lost_HMS_HMS) + merged_star.total_mass_h1 = merged_star.center_h1 * merged_star.mass # very simplified assumption of X_average_merged = merged_star.center_h1 + #^ should work ok for both center_h1 and total_mass_h1 matching with single stars models for key in STARPROPERTIES: # these stellar attributes become np.nan @@ -204,8 +280,6 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_POSTMS and s2 in LIST_ACCEPTABLE_STATES_FOR_HMS): - merged_star.mass = star_base.mass + comp.mass #TODO: in step_CEE we need to eject part of the (common) envelope - # weigheted mixing on the surface abundances of the whole comp with the envelope of star_base merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1", mass_weight1="H-rich_envelope_mass", mass_weight2="mass") merged_star.surface_he4 = mass_weighted_avg(abundance_name = "surface_he4", mass_weight1="H-rich_envelope_mass", mass_weight2="mass") @@ -213,13 +287,15 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.surface_n14 = mass_weighted_avg(abundance_name = "surface_n14", mass_weight1="H-rich_envelope_mass", mass_weight2="mass") merged_star.surface_o16 = mass_weighted_avg(abundance_name = "surface_o16", mass_weight1="H-rich_envelope_mass", mass_weight2="mass") + # The change of masses occurs after the calculation of weighted averages + merged_star.mass = star_base.mass + comp.mass + for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.log_LHe = star_base.log_LHe @@ -232,7 +308,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_HMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - merged_star.mass = star_base.mass + comp.mass + merged_star = comp merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1", mass_weight2="H-rich_envelope_mass", mass_weight1="mass") merged_star.surface_he4 = mass_weighted_avg(abundance_name = "surface_he4", mass_weight2="H-rich_envelope_mass", mass_weight1="mass") @@ -240,51 +316,54 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.surface_n14 = mass_weighted_avg(abundance_name = "surface_n14", mass_weight2="H-rich_envelope_mass", mass_weight1="mass") merged_star.surface_o16 = mass_weighted_avg(abundance_name = "surface_o16", mass_weight2="H-rich_envelope_mass", mass_weight1="mass") + # The change of masses occurs after the calculation of weighted averages + merged_star.mass = star_base.mass + comp.mass + for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.log_LHe = comp.log_LHe merged_star.log_LZ = comp.log_LZ merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! - massless_remnant = convert_star_to_massless_remnant(comp) + massless_remnant = convert_star_to_massless_remnant(star_base) #postMS + postMS elif (s1 in LIST_ACCEPTABLE_STATES_FOR_POSTMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, key) - setattr(merged_star, key,current) - # weighted central abundances if merging cores. Else only from star_base if star_base.co_core_mass == 0 and comp.co_core_mass == 0: # two stars with Helium cores merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") #TODO : maybe he_core_mass makes more sense? merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") @@ -295,14 +374,17 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.surface_n14 = mass_weighted_avg(abundance_name = "surface_n14", mass_weight1="H-rich_envelope_mass", mass_weight2="H-rich_envelope_mass") merged_star.surface_o16 = mass_weighted_avg(abundance_name = "surface_o16", mass_weight1="H-rich_envelope_mass", mass_weight2="H-rich_envelope_mass") + # add total and core masses after calculations of weighter average + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, key) + setattr(merged_star, key,current) for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -311,30 +393,32 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_POSTMS and s2 in LIST_ACCEPTABLE_STATES_FOR_HeMS): - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, key) - setattr(merged_star, key,current) - # weighted central abundances if merging cores. Else only from star_base if star_base.co_core_mass == 0 and comp.co_core_mass == 0: # two stars with only Helium cores, not CO cores merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has just a He core (is a HeMS star) pass # the central abundances are kept as the ones of star_base else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") + # surface abundances from star_base + + # add total and core masses after abundance mass weighted calculations + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, key) + setattr(merged_star, key,current) for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -343,12 +427,9 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_HeMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, key) - setattr(merged_star, key,current) + merged_star = comp - # weighted central abundances if merging cores. Else only from star_base + # weighted central abundances if merging cores. Else only from comp if star_base.co_core_mass == 0 and comp.co_core_mass == 0: # two stars with only Helium cores, not CO cores merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") @@ -356,62 +437,72 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") elif (star_base.co_core_mass == 0 and comp.co_core_mass > 0): # star_base is the HeMS Star and comp has a CO core - pass # the central abundances are kept as the ones of star_base + pass # the central abundances are kept as the ones of comp else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") + # surface abundances from comp + + # add total and core masses after weighted averages above + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, key) + setattr(merged_star, key,current) for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! - massless_remnant = convert_star_to_massless_remnant(comp) + massless_remnant = convert_star_to_massless_remnant(star_base) #postMS + HeStar that is not in HeMS elif (s1 in LIST_ACCEPTABLE_STATES_FOR_POSTMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTHeMS): - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, key) - setattr(merged_star, key,current) - # weighted central abundances if merging cores. Else only from star_base if star_base.co_core_mass == 0 and comp.co_core_mass == 0: # two stars with Helium cores merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") + # add total and core masses + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, key) + setattr(merged_star, key,current) + for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -420,25 +511,25 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in LIST_ACCEPTABLE_STATES_FOR_POSTHeMS and s2 in LIST_ACCEPTABLE_STATES_FOR_POSTMS): - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, key) - setattr(merged_star, key,current) + merged_star = comp + # weighted central abundances if merging cores. Else only from star_base if star_base.co_core_mass == 0 and comp.co_core_mass == 0: # two stars with Helium cores merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core - pass # the central abundances are kept as the ones of star_base + merged_star.center_h1 = star_base.center_h1 + merged_star.center_he4 = star_base.center_he4 + merged_star.center_c12 = star_base.center_c12 + merged_star.center_n14 = star_base.center_n14 + merged_star.center_o16 = star_base.center_o16 elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core - merged_star.center_h1 = comp.center_h1 - merged_star.center_he4 = comp.center_he4 - merged_star.center_c12 = comp.center_c12 - merged_star.center_n14 = comp.center_n14 - merged_star.center_o16 = comp.center_o16 + pass # the central abundances are kept as the ones of comp elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") @@ -448,50 +539,60 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") + # add total and core masses + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, key) + setattr(merged_star, key,current) + for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! - massless_remnant = convert_star_to_massless_remnant(comp) + massless_remnant = convert_star_to_massless_remnant(star_base) # HeStar + HeStar elif (s1 in STAR_STATES_HE_RICH and s2 in STAR_STATES_HE_RICH): - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, key) - setattr(merged_star, key,current) - # weighted central abundances if merging cores. Else only from star_base if star_base.co_core_mass == 0 and comp.co_core_mass == 0: # two stars with Helium cores merged_star.center_h1 = mass_weighted_avg(mass_weight1="he_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="he_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="he_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="he_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="he_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="he_core_mass") + setattr(merged_star, "center_gamma", np.nan) elif (star_base.co_core_mass > 0 and comp.co_core_mass == 0): # star_base with CO core and the comp has a He core pass # the central abundances are kept as the ones of star_base elif (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has a He core merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 + merged_star.center_gamma = comp.center_gamma elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") + # add total and core masses + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, key) + setattr(merged_star, key,current) + # weigheted mixing on the surface abundances based on the He-rich envelopes of the two stars merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1", mass_weight1="He-rich_envelope_mass", mass_weight2="He-rich_envelope_mass") merged_star.surface_he4 = mass_weighted_avg(abundance_name = "surface_he4", mass_weight1="He-rich_envelope_mass", mass_weight2="He-rich_envelope_mass") @@ -504,8 +605,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) @@ -514,34 +614,40 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma elif (s1 in STAR_STATES_NOT_CO and s2 in ["WD"]): - #WD is considered a stripped CO core - # add total and core masses - for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: - current = getattr(star_base, key) + getattr(comp, "mass") - setattr(merged_star, key,current) + # WD is considered a stripped CO core + # WD would always be the comp, it cannot be the engulfing star (so no need to do the opposite stars case below) + # weighted central abundances if merging cores. Else only from star_base - if (comp.co_core_mass > 0 and star_base.co_core_mass == 0): # comp with CO core and the star_base has not + if (comp.co_core_mass is not None and star_base.co_core_mass == 0): # comp with CO core and the star_base has not merged_star.center_h1 = comp.center_h1 merged_star.center_he4 = comp.center_he4 merged_star.center_c12 = comp.center_c12 + merged_star.avg_c_in_c_core = comp.avg_c_in_c_core merged_star.center_n14 = comp.center_n14 merged_star.center_o16 = comp.center_o16 - elif (star_base.co_core_mass > 0 and comp.co_core_mass > 0): + merged_star.center_gamma = comp.center_gamma + elif (comp.co_core_mass is not None and star_base.co_core_mass > 0): merged_star.center_h1 = mass_weighted_avg(mass_weight1="co_core_mass") merged_star.center_he4 = mass_weighted_avg(abundance_name = "center_he4", mass_weight1="co_core_mass") merged_star.center_c12 = mass_weighted_avg(abundance_name = "center_c12", mass_weight1="co_core_mass") + merged_star.avg_c_in_c_core = mass_weighted_avg(abundance_name = "avg_c_in_c_core", mass_weight1="co_core_mass") merged_star.center_n14 = mass_weighted_avg(abundance_name = "center_n14", mass_weight1="co_core_mass") merged_star.center_o16 = mass_weighted_avg(abundance_name = "center_o16", mass_weight1="co_core_mass") + setattr(merged_star, "center_gamma", np.nan) else: Pwarn("weird compbination of CO core masses during merging", "EvolutionWarning") + # add total and core masses + for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: + current = getattr(star_base, key) + getattr(comp, "mass") + setattr(merged_star, key,current) + for key in STARPROPERTIES: # these stellar attributes become np.nan for substring in ["log_", "lg_", "surf_", "conv_", "lambda", "profile"]: if (substring in key) : setattr(merged_star, key, np.nan) - if key in [ "c12_c12", "center_gamma", - "avg_c_in_c_core", "total_moment_of_inertia", "spin"]: + if key in [ "c12_c12", "total_moment_of_inertia", "spin"]: setattr(merged_star, key, np.nan) merged_star.state = check_state_of_star(merged_star, star_CO=False) # TODO for sure this needs testing! massless_remnant = convert_star_to_massless_remnant(comp) diff --git a/posydon/binary_evol/DT/track_match.py b/posydon/binary_evol/DT/track_match.py index 07399cb708..ca54e85f3e 100644 --- a/posydon/binary_evol/DT/track_match.py +++ b/posydon/binary_evol/DT/track_match.py @@ -52,10 +52,10 @@ val_names = [" ", "mass", "log_R", "center_h1", "surface_h1", "he_core_mass", "center_he4", "surface_he4", - "center_c12"] + "center_c12", "co_core_mass", "total_mass_h1"] str_fmts = ["{:>14}", "{:>9}","{:>9}", "{:>9}", "{:>10}", "{:>12}", - "{:>10}", "{:>11}", "{:>10}"] + "{:>10}", "{:>11}", "{:>10}", "{:>12}", "{:>13}"] row_str = " ".join(str_fmts) DIVIDER_STR = "_"*len(row_str.format(*[""]*len(str_fmts))) # MAJOR.MINOR version of imported scipy package @@ -285,7 +285,9 @@ def __init__( "surface_he4", "surface_h1", "log_R", - "center_c12" + "center_c12", + "co_core_mass", + "total_mass_h1" ] ) @@ -743,7 +745,9 @@ def match_to_single_star(self, star): f'surface_he4 = {star.surface_he4:.4f}, ', f'surface_h1 = {star.surface_h1:.4f}, ', f'he_core_mass = {star.he_core_mass:.3f}, ', - f'center_c12 = {star.center_c12:.4f}' + f'center_c12 = {star.center_c12:.4f},', + f'co_core_mass = {star.co_core_mass:.3f}', + f'total_mass_h1 = {star.total_mass_h1:.3f}' ) # done with matching attempts diff --git a/posydon/popsyn/population_params_default.ini b/posydon/popsyn/population_params_default.ini index 2e06a9b3c4..f569d9c6ba 100644 --- a/posydon/popsyn/population_params_default.ini +++ b/posydon/popsyn/population_params_default.ini @@ -187,6 +187,10 @@ absolute_import = None # if given, use an absolute filepath to user defined step: # ['', ''] + rel_mass_lost_HMS_HMS = "Glebbeek+2013" + # [0-1) or "Glebbeek+2013" (dependent on mass ratio q) + HMS_HMS_merging_rejuvenation = True + # if True then new abundances based on total_mass_h1 or he4 (according to Schneider+2016) record_matching = False # True, False verbose = False