diff --git a/posydon/binary_evol/DT/step_detached.py b/posydon/binary_evol/DT/step_detached.py index f76728913c..03c10844db 100644 --- a/posydon/binary_evol/DT/step_detached.py +++ b/posydon/binary_evol/DT/step_detached.py @@ -387,7 +387,8 @@ def __init__( "surface_he4", "surface_h1", "log_R", - "center_c12" + "center_c12", + "co_core_mass" ] ) @@ -899,14 +900,15 @@ def sq_diff_function(x): f'{initials[0]:.3f} [Msun],', f'{initials[1]/1e6:.3f} [Myrs]', "\n", "with m(t0), log10(R(t0), center_he(t0), surface_he4(t0), " - "surface_h1(t0), he_core_mass(t0), center_c12(t0) = \n", + "surface_h1(t0), he_core_mass(t0), center_c12(t0), co_core_mass(t0) = \n", f'{self.get_track_val("mass", htrack, *sol.x):.3f}', f'{self.get_track_val("log_R", htrack, *sol.x):.3f}', f'{self.get_track_val("center_he4", htrack, *sol.x):.4f}', f'{self.get_track_val("surface_he4", htrack, *sol.x):.4f}', f'{self.get_track_val("surface_h1", htrack, *sol.x):.4f}', f'{self.get_track_val("he_core_mass", htrack, *sol.x):.3f}', - f'{self.get_track_val("center_c12", htrack, *sol.x):.4f}\n', + f'{self.get_track_val("center_c12", htrack, *sol.x):.4f}', + f'{self.get_track_val("co_core_mass", htrack, *sol.x):.3f}\n', "The same values of the secondary at the end of the previous " "step was = \n", f'{star.mass:.3f}', @@ -915,7 +917,8 @@ def sq_diff_function(x): f'{star.surface_he4:.4f}', f'{star.surface_h1:.4f}', f'{star.he_core_mass:.3f}', - f'{star.center_c12:.4f}' + f'{star.center_c12:.4f}', + f'{star.co_core_mass:.4f}' ) return initials[0], initials[1], htrack @@ -1082,15 +1085,15 @@ def get_star_data(binary, star1, star2, htrack, if self.verbose or self.verbose == 1: print("Matching duration: " f"{t_after_matching-t_before_matching:.6g}") - + if pd.isna(m0) or pd.isna(t0): return None, None, None - + if htrack: self.grid = self.grid_Hrich else: self.grid = self.grid_strippedHe - + # check if m0 is in the grid if m0 < self.grid.grid_mass.min() or m0 > self.grid.grid_mass.max(): set_binary_to_failed(binary) diff --git a/posydon/binary_evol/DT/step_merged.py b/posydon/binary_evol/DT/step_merged.py index 70ca96e425..7fafdcb8ab 100644 --- a/posydon/binary_evol/DT/step_merged.py +++ b/posydon/binary_evol/DT/step_merged.py @@ -15,6 +15,7 @@ from posydon.utils.common_functions import check_state_of_star from posydon.binary_evol.DT.step_isolated import IsolatedStep from posydon.utils.posydonerror import FlowError +import copy from posydon.utils.posydonwarning import Pwarn @@ -88,8 +89,8 @@ def __call__(self,binary): merged_star_properties = self.merged_star_properties 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: ", 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_2.surface_he4) + 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, 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_2.center_he4, star_1.center_c12, star_2.center_c12, 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": if binary.event == 'oMerging1': binary.star_1,binary.star_2 = merged_star_properties(binary.star_1,binary.star_2) @@ -100,12 +101,16 @@ def __call__(self,binary): else: raise FlowError("step_merging initiated but binary.state != 'merged'") - binary.event = None 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, co_core_mass merged: ", binary.star_1.mass, binary.star_1.he_core_mass, binary.star_1.co_core_mass) + print("center_h1, center_he4, center_c12, surface_h1, surface_he4: ", binary.star_1.center_h1, binary.star_1.center_he4, binary.star_1.center_c12, binary.star_1.surface_h1, binary.star_1.surface_he4) + elif binary.event == 'oMerging2': + print("M_merged , he_core_mass, co_core_mass merged: ", binary.star_2.mass, binary.star_2.he_core_mass, binary.star_2.co_core_mass) + print("center_h1, center_he4, center_c12, surface_h1, surface_he4: ", binary.star_2.center_h1, binary.star_2.center_he4, binary.star_2.center_c12, binary.star_2.surface_h1, binary.star_2.surface_he4) + binary.event = None super().__call__(binary) def merged_star_properties(self,star_base,comp): @@ -120,8 +125,7 @@ def merged_star_properties(self,star_base,comp): is the star that is engulfed """ #by default the stellar attributes that keep the same value from the - #merged_star = copy.copy(star_base) - merged_star = star_base + merged_star = copy.deepcopy(star_base) s1 = star_base.state s2 = comp.state @@ -143,7 +147,19 @@ 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_abundance", 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 @@ -194,8 +210,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.log_LHe = star_base.log_LHe @@ -208,7 +223,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 = comp + merged_star = copy.deepcopy(comp) merged_star.mass = star_base.mass + comp.mass merged_star.surface_h1 = mass_weighted_avg(abundance_name = "surface_h1", mass_weight2="H-rich_envelope_mass", mass_weight1="mass") @@ -222,8 +237,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.log_LHe = comp.log_LHe @@ -246,22 +260,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma 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") @@ -278,8 +298,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) @@ -298,8 +317,10 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma 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: @@ -310,8 +331,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) @@ -320,7 +340,7 @@ 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): - merged_star = comp + merged_star = copy.deepcopy(comp) # add total and core masses for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: @@ -332,8 +352,10 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma 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 is the HeMS Star and comp has a CO core pass # the central abundances are kept as the ones of star_base else: @@ -344,8 +366,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(star_base) @@ -364,22 +385,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma 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") @@ -388,8 +415,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) @@ -398,7 +424,7 @@ 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): - merged_star = comp + merged_star = copy.deepcopy(comp) # add total and core masses for key in ["mass", "he_core_mass", "c_core_mass", "o_core_mass", "co_core_mass"]: current = getattr(merged_star, key) + getattr(star_base, key) @@ -408,22 +434,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma 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") @@ -432,8 +464,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(star_base) @@ -452,22 +483,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma 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") @@ -483,8 +520,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) @@ -496,21 +532,28 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma #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(merged_star, key) + getattr(comp, "mass") + if key not in ["c_core_mass", "o_core_mass"]: + current = getattr(merged_star, key) + getattr(comp, key) + else: + current = getattr(merged_star, key) + getattr(comp, "co_core_mass") setattr(merged_star, key,current) # 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") @@ -519,8 +562,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) @@ -528,7 +570,7 @@ def mass_weighted_avg(star1=star_base,star2=comp, abundance_name="center_h1", ma # Star + NS/BH elif (s1 in STAR_STATES_NOT_CO and s2 in ["NS", "BH"]): - merged_star = comp + merged_star = copy.deepcopy(comp) # TODO: potentially flag a Thorne-Zytkov object massless_remnant = convert_star_to_massless_remnant(star_base) diff --git a/posydon/binary_evol/MESA/step_mesa.py b/posydon/binary_evol/MESA/step_mesa.py index c8ec84fc29..caf5fe296d 100644 --- a/posydon/binary_evol/MESA/step_mesa.py +++ b/posydon/binary_evol/MESA/step_mesa.py @@ -640,7 +640,7 @@ def update_properties_NN(self, star_1_CO=False, star_2_CO=False, getattr(star, key_h).extend(cb_bh[key_p][:-1]) elif key in ['state', 'lg_mdot']: continue - elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16']: + elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16', 'avg_c_in_c_core']: continue else: setattr(star, key, None) @@ -886,7 +886,7 @@ def initial_final_interpolation(self, star_1_CO=False, star_2_CO=False): setattr(star, key, fv[key_p]) elif key == 'state': continue - elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16']: + elif star.state == 'WD' and key in ['co_core_mass','he_core_mass','center_h1','center_he4','center_c12','center_n14','center_o16', 'avg_c_in_c_core']: continue else: setattr(star, key, None)