"""Detached evolution step."""
__authors__ = [
"Devina Misra <devina.misra@unige.ch>",
"Zepei Xing <Zepei.Xing@unige.ch>",
"Emmanouil Zapartas <ezapartas@gmail.com>",
"Nam Tran <tranhn03@gmail.com>",
"Simone Bavera <Simone.Bavera@unige.ch>",
"Konstantinos Kovlakas <Konstantinos.Kovlakas@unige.ch>",
"Kyle Akira Rocha <kylerocha2024@u.northwestern.edu>",
"Jeffrey Andrews <jeffrey.andrews@northwestern.edu>",
]
import os
import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import PchipInterpolator
from scipy.optimize import minimize
from scipy.optimize import root
from posydon.utils.data_download import PATH_TO_POSYDON_DATA
from posydon.binary_evol.binarystar import BINARYPROPERTIES
from posydon.binary_evol.singlestar import STARPROPERTIES
from posydon.interpolation import GRIDInterpolator
from posydon.interpolation.data_scaling import DataScaler
from posydon.utils.common_functions import (
bondi_hoyle,
orbital_period_from_separation,
roche_lobe_radius,
check_state_of_star,
PchipInterpolator2
)
from posydon.binary_evol.flow_chart import (STAR_STATES_CC)
import posydon.utils.constants as const
LIST_ACCEPTABLE_STATES_FOR_HMS = ["H-rich_Core_H_burning"]
LIST_ACCEPTABLE_STATES_FOR_postMS = [
"H-rich_Shell_H_burning",
"H-rich_Core_He_burning",
"H-rich_Central_He_depleted",
"H-rich_Core_C_burning",
"H-rich_Central_C_depletion",
"H-rich_non_burning"]
LIST_ACCEPTABLE_STATES_FOR_HeStar = [
'stripped_He_Core_He_burning',
'stripped_He_Shell_He_burning', # includes stars burning C in core
'stripped_He_Central_He_depleted', # includes stars burning C in core
'stripped_He_Central_C_depletion',
'stripped_He_non_burning' # includes stars burning C in core
]
STAR_STATES_H_RICH = [
'H-rich_Core_H_burning',
'H-rich_Core_He_burning',
'H-rich_Shell_H_burning',
'H-rich_Central_He_depleted',
'H-rich_Shell_He_burning',
'H-rich_Core_C_burning',
'H-rich_Central_C_depletion',
'H-rich_non_burning'
]
[docs]class detached_step:
"""Evolve a detached binary.
The binary will be evolved until Roche-lobe overflow, core-collapse or
maximum simulation time, using the standard equations that govern the
orbital evolution.
Parameters
----------
path : str
Path to the directory that contains a HDF5 grid.
dt : float
The timestep size, in years, to be appended to the history of the
binary. None means only the final step.
Note: do not select very small timesteps cause it may mess with the
solving of the ODE.
n_o_steps_history: int
Alternatively, we can define the number of timesteps to be appended to
the history of the binary. None means only the final step. If both `dt`
and `n_o_steps_history` are different than None, `dt` has priority.
matching_method: str
Method to find the best match between a star from a previous step and a
point in a single MIST-like stellar track. Options "root" (which tries
to find a root of two matching quantities, and it is possible to not
achieve it) or "minimize" (minimizes the sum of squares of differences
of various quantities between the previous step and the track).
verbose : Boolean
True if we want to print stuff.
do_wind_loss: Boolean
If True, take into account change of separation due to mass loss from
the star.
do_tides: Booleans
If True, take into account change of separation, eccentricity and star
spin due to tidal forces.
do_gravitational_radiation: Boolean
If True, take into account change of separation and eccentricity due to
gravitational wave radiation.
do_magnetic_braking: Boolean
If True, take into account change of star spin due to magnetic braking.
do_stellar_evolution_and_spin_from_winds: Boolean
If True, take into account change of star spin due to change of its
moment of inertia during its evolution and due to spin angular momentum
loss due to winds.
Attributes
----------
KEYS : list of str
Contains valid keywords which is used
to extract quantities from the grid.
grid : GRIDInterpolator
Object to interpolate between the time-series in
the h5 grid.
initial_mass : list of float
Contains the initial masses of the stars in the grid.
Note
----
A matching between the properties of the star, and the h5 tracks are
required. In the "root" solver matching_method, if the root solver fails
then the evolution will immediately end, and the binary state will be
tagged with "Root solver failed". In the "minimize" matching_method, we
minimize the sum of squares of differences of various quantities between
the previous step and the h5 track.
Warns
-----
UserWarning
If the call cannot determine the primary or secondary in the binary.
Raises
------
Exception
If the ode-solver fails to solve the differential equation
that governs the orbital evolution.
"""
def __init__(
self,
grid=None,
path=PATH_TO_POSYDON_DATA,
dt=None,
n_o_steps_history=None,
matching_method="minimize",
initial_mass=None,
rootm=None,
verbose=False,
do_wind_loss=True,
do_tides=True,
do_gravitational_radiation=True,
do_magnetic_braking=True,
do_stellar_evolution_and_spin_from_winds=True,
RLO_orbit_at_orbit_with_same_am=False
):
"""Initialize the step. See class documentation for details."""
self.dt = dt
self.n_o_steps_history = n_o_steps_history
self.matching_method = matching_method
self.do_wind_loss = do_wind_loss
self.do_tides = do_tides
self.do_gravitational_radiation = do_gravitational_radiation
self.do_magnetic_braking = do_magnetic_braking
self.do_stellar_evolution_and_spin_from_winds = (
do_stellar_evolution_and_spin_from_winds
)
self.RLO_orbit_at_orbit_with_same_am = RLO_orbit_at_orbit_with_same_am
self.grid = grid
self.initial_mass = initial_mass
self.rootm = rootm
self.verbose = verbose
if verbose:
print(
dt,
n_o_steps_history,
matching_method,
do_wind_loss,
do_tides,
do_gravitational_radiation,
do_magnetic_braking,
do_stellar_evolution_and_spin_from_winds)
self.translate = {
"time": "time",
"orbital_period": "porb",
"eccentricity": "ecc",
"separation": "sep",
"state": None,
"event": None,
"rl_relative_overflow_1": "rl_relative_overflow_1",
"rl_relative_overflow_2": "rl_relative_overflow_2",
"lg_mtransfer_rate": "lg_mtransfer_rate",
"V_sys": None,
"mass": "mass",
"log_R": "log_R",
"R": "R",
"lg_mdot": "mdot",
"log_L": "log_L",
"lg_wind_mdot": "mdot",
"lg_system_mdot": "lg_mdot",
"he_core_mass": "he_core_mass",
"he_core_radius": "he_core_radius",
"c_core_mass": "c_core_mass",
"c_core_radius": "c_core_radius",
"o_core_mass": "o_core_mass",
"o_core_radius": "o_core_radius",
"center_h1": "center_h1",
"center_he4": "center_he4",
"center_c12": "center_c12",
"center_o16": "center_o16",
"center_n14": "center_n14",
"surface_h1": "surface_h1",
"surface_he4": "surface_he4",
"surface_c12": "surface_c12",
"surface_n14": "surface_n14",
"surface_o16": "surface_o16",
"center_gamma": "center_gamma",
"log_LH": "log_LH",
"log_LHe": "log_LHe",
"log_LZ": "log_LZ",
"log_Lnuc": "log_Lnuc",
"c12_c12": "c12_c12",
"avg_c_in_c_core": "avg_c_in_c_core",
"surf_avg_omega_div_omega_crit": "surf_avg_omega_div_omega_crit",
"surf_avg_omega": "omega",
"total_moment_of_inertia": "inertia",
"log_total_angular_momentum": "log_total_angular_momentum",
"profile": None,
"metallicity": None,
"spin": "spin_parameter",
"log_total_angular_momentum": "log_total_angular_momentum",
"conv_env_top_mass": "conv_env_top_mass",
"conv_env_bot_mass": "conv_env_bot_mass",
"conv_env_top_radius": "conv_env_top_radius",
"conv_env_bot_radius": "conv_env_bot_radius",
"conv_env_turnover_time_g": "conv_env_turnover_time_g",
"conv_env_turnover_time_l_b": "conv_env_turnover_time_l_b",
"conv_env_turnover_time_l_t": "conv_env_turnover_time_l_t",
"envelope_binding_energy": "envelope_binding_energy",
"mass_conv_reg_fortides": "mass_conv_reg_fortides",
"thickness_conv_reg_fortides": "thickness_conv_reg_fortides",
"radius_conv_reg_fortides": "radius_conv_reg_fortides",
"lambda_CE_1cent": "lambda_CE_1cent",
"lambda_CE_10cent": "lambda_CE_10cent",
"lambda_CE_30cent": "lambda_CE_30cent",
"co_core_mass": "co_core_mass",
"co_core_radius": "co_core_radius",
"lambda_CE_pure_He_star_10cent": "lambda_CE_pure_He_star_10cent",
"trap_radius": "trap_radius",
"acc_radius": "acc_radius",
"t_sync_rad_1": "t_sync_rad_1",
"t_sync_conv_1": "t_sync_conv_1",
"t_sync_rad_2": "t_sync_rad_2",
"t_sync_conv_2": "t_sync_conv_2",
"mass_transfer_case": None,
"nearest_neighbour_distance": None,
}
# these are the KEYS read from POSYDON h5 grid files (after translating
# them to the appropriate columns)
self.KEYS = (
'age',
'mass',
'mdot',
'inertia',
'conv_mx1_top_r',
'conv_mx1_bot_r',
'surface_h1',
'center_h1',
'mass_conv_reg_fortides',
'thickness_conv_reg_fortides',
'radius_conv_reg_fortides',
'log_Teff',
'surface_he3',
'surface_he4',
'center_he4',
'avg_c_in_c_core',
'log_LH',
'log_LHe',
'log_LZ',
'log_Lnuc',
'c12_c12',
'center_c12',
'he_core_mass',
'log_L',
'log_R',
'c_core_mass',
'o_core_mass',
'co_core_mass',
'c_core_radius',
'o_core_radius',
'co_core_radius',
'spin_parameter',
'log_total_angular_momentum',
'center_n14',
'center_o16',
'surface_n14',
'surface_o16',
'conv_env_top_mass',
'conv_env_bot_mass',
'conv_env_top_radius',
'conv_env_bot_radius',
'conv_env_turnover_time_g',
'conv_env_turnover_time_l_b',
'conv_env_turnover_time_l_t',
'envelope_binding_energy',
'lambda_CE_1cent',
'lambda_CE_10cent',
'lambda_CE_30cent',
'lambda_CE_pure_He_star_10cent',
'center_gamma'
)
self.KEYS_POSITIVE = (
'mass_conv_reg_fortides',
'thickness_conv_reg_fortides',
'radius_conv_reg_fortides'
)
self.root_keys = np.array( # for the matching
[
"age",
"mass",
"he_core_mass",
"center_h1",
"center_he4",
"surface_he4",
"surface_h1",
"log_R",
"center_c12"
]
)
# keys for the final value interpolation
self.final_keys = (
'avg_c_in_c_core_at_He_depletion',
'co_core_mass_at_He_depletion',
'm_core_CE_1cent',
'm_core_CE_10cent',
'm_core_CE_30cent',
'm_core_CE_pure_He_star_10cent',
'r_core_CE_1cent',
'r_core_CE_10cent',
'r_core_CE_30cent',
'r_core_CE_pure_He_star_10cent'
)
# keys for the star profile interpolation
self.profile_keys = (
'radius',
'mass',
'logRho',
'energy',
'x_mass_fraction_H',
'y_mass_fraction_He',
'z_mass_fraction_metals',
'neutral_fraction_H',
'neutral_fraction_He',
'avg_charge_He'
)
grid_name1 = os.path.join('single_HMS', 'grid_0.0142.h5')
grid_name2 = os.path.join('single_HeMS', 'grid_0.0142.h5')
self.grid1 = GRIDInterpolator(os.path.join(path, grid_name1))
self.grid2 = GRIDInterpolator(os.path.join(path, grid_name2))
[docs] def get_track_val(self, key, htrack, m0, t):
"""Return a single value from the interpolated time-series.
Parameters
----------
key : str
Keyword of the required quantity.
m0 : float
The associated initial mass of the required quantity.
t : float
The required time in the time-series.
Returns
-------
float
The value of the quantity `key` from a MIST-like track of
initial mass `m0` at the time `t0`.
"""
# htrack as a boolean determines whether H or He grid is used
if htrack:
self.grid = self.grid1
else:
self.grid = self.grid2
try:
x = self.grid.get("age", m0)
y = self.grid.get(key, m0)
except ValueError:
return np.array(t) * np.nan
try:
val = np.interp(t, x, y, left=1e99, right=1e99)
except ValueError:
i_bad = [None]
while len(i_bad):
i_bad = np.where(np.diff(x) <= 0)[0]
x = np.delete(x, i_bad)
y = np.delete(y, i_bad)
val = np.interp(t, x, y)
return val
[docs] def scale(self, key, htrack, method):
"""Nomarlize quantities in the single star grids to (0,1).
Parameters
----------
key : str
Keyword of the required quantity.
method : str
Scalling method in the data normalization class
Returns
-------
class
Data normalization class
"""
if htrack:
self.grid = self.grid1
else:
self.grid = self.grid2
self.initial_mass = self.grid.grid_mass
all_attribute = []
for mass in self.initial_mass:
for i in self.grid.get(key, mass):
all_attribute.append(i)
all_value = np.array(all_attribute)
sc = DataScaler()
xt = sc.fit_and_transform(
all_value, method=method, lower=0.0, upper=1.0)
# xtnew = sc.transform(x)
return sc
[docs] def get_root0(self, keys, x, htrack, rs=None):
"""Determine the closest associated initial mass and time in the grid
which has a specific value as tentative solutions for matching.
Parameters
----------
keys : list of str
Contains the keys of the required specific quantities that will be
matched in the MIST-like track.
x : list of floats, of same length as "keys"
Contains the latest values (from a previous POSYDON step) of the
quantities of "keys" in the POSYDON SingleStar object.
rs : list of floats, same length as "keys"
Contains normalization factors to be divided for rescaling
x values.
Returns
-------
list of 2 float values
Contains the associated initial mass (in solar units) and the time
(in years) such that the time-series of the `keys` at that time has
the closest values to `x`. These will become m0, t0 for the later
integration during the detached binary evolution.
If there is no match then NaNs will be returned instead.
"""
if htrack:
self.grid = self.grid1
else:
self.grid = self.grid2
self.initial_mass = self.grid.grid_mass
n = 0
for mass in self.grid.grid_mass:
n = max(n, len(self.grid.get("age", mass)))
self.rootm = np.inf * np.ones((len(self.grid.grid_mass),
n, len(self.root_keys)))
for i, mass in enumerate(self.grid.grid_mass):
for j, key in enumerate(self.root_keys):
track = self.grid.get(key, mass)
self.rootm[i, : len(track), j] = track
if rs is None:
rs = np.ones_like(keys)
else:
rs = np.asanyarray(rs)
x = np.asanyarray(x)
idx = np.argmax(np.asanyarray(keys)[:, None] == self.root_keys, axis=1)
X = self.rootm[:, :, idx]
d = np.linalg.norm((X - x[None, None, :]) / rs[None, None, :], axis=-1)
idx = np.unravel_index(d.argmin(), X.shape[:-1])
t = self.rootm[idx][np.argmax("age" == self.root_keys)]
m0 = self.grid.grid_mass[idx[0]]
return m0, t
[docs] def get_mist0(self, star, htrack):
"""Determine the associated initial mass and time in the grid that
matches the properties of the binary.
For "root" matching_method, the properties that are matched is always
the mass of the secondary star.
If the secondary has the state `MS` then
the center hydrogen abundance will also be matched
otherwise the mass of helium-core will be matched.
Parameters
----------
star : SingleStar
The star which properties are required
to be matched with the single MIST-like grid.
Returns
-------
list of 2 float values
Contains the associated (in solar units) and the time (in years)
such that the time-series in the grid matches
the properties of the secondary.
"""
if htrack:
self.grid = self.grid1
else:
self.grid = self.grid2
m_min_H = np.min(self.grid1.grid_mass)
m_max_H = np.max(self.grid1.grid_mass)
m_min_He = np.min(self.grid2.grid_mass)
m_max_He = np.max(self.grid2.grid_mass)
get_root0 = self.get_root0
get_track_val = self.get_track_val
matching_method = self.matching_method
tran = self.transform()
sc_mass_H = tran[0]
sc_mass_He = tran[1]
sc_log_R_H = tran[2]
sc_log_R_He = tran[3]
sc_he_core_mass_H = tran[4]
sc_he_core_mass_He = tran[5]
sc_center_h1 = tran[6]
sc_center_he4_H = tran[7]
sc_center_he4_He = tran[8]
sc_center_c12 = tran[9]
initials = None
# tolerance 1e-8
tolerance_mist_integration = 1e-2
if self.verbose:
print(matching_method)
if matching_method == "root":
if star.state in LIST_ACCEPTABLE_STATES_FOR_HMS:
x0 = get_root0(["center_h1", "mass"],
[star.center_h1, star.mass],
htrack, rs=[0.7, 300])
sol = root(
lambda x: [
get_track_val("center_h1", htrack, *x)
- star.center_h1,
get_track_val("mass", htrack, *x) - star.mass,
],
x0,
method="hybr",
)
else:
x0 = get_root0(
["he_core_mass", "mass"],
[star.he_core_mass, star.mass],
htrack,
rs=[11, 300],
)
sol = root(
lambda x: [
get_track_val("he_core_mass", htrack, *x)
- star.he_core_mass,
get_track_val("mass", htrack, *x) - star.mass,
],
x0,
method="hybr",
)
if not sol.success or sol.x[1] < 0:
initials = (np.nan, np.nan)
else:
initials = sol.x
elif matching_method == "minimize":
if star.state in LIST_ACCEPTABLE_STATES_FOR_HMS:
# Initialezed ZAMS systems have no information of radius
if (star.log_R is None) or np.isnan(star.log_R):
initials = (star.mass, 0)
else:
MESA_label = ["mass", "center_h1", "log_R", "he_core_mass"]
posydon_attribute = [star.mass, star.center_h1,
star.log_R, star.he_core_mass]
rs = [20.0, 1.0, 2.0, 10.0]
if self.verbose:
print("Matching attributes and their normalizations :",
MESA_label, rs)
for i in MESA_label:
if i not in self.root_keys:
raise Exception("Expected matching parameter not "
"added in the MIST model options.")
x0 = get_root0(MESA_label, posydon_attribute,
htrack, rs=rs)
bnds = ([m_min_H, m_max_H], [0, None])
sol = minimize(
lambda x: (
sc_mass_H.transform(
get_track_val(MESA_label[0], htrack, *x))
- sc_mass_H.transform(posydon_attribute[0])) ** 2
+ (sc_center_h1.transform(
get_track_val(MESA_label[1], htrack, *x))
- sc_center_h1.transform(posydon_attribute[1]))**2
+ (sc_log_R_H.transform(
get_track_val(MESA_label[2], htrack, *x))
- sc_log_R_H.transform(posydon_attribute[2])) ** 2
+ (sc_he_core_mass_H.transform(
get_track_val(MESA_label[3], htrack, *x))
- sc_he_core_mass_H.transform(
posydon_attribute[3])) ** 2,
x0,
method="TNC",
bounds=bnds
)
# stars after mass transfer could swell up so that log_R
# is not appropriate for matching
if np.abs(sol.fun) > tolerance_mist_integration:
MESA_label = ["mass", "center_h1", "he_core_mass"]
posydon_attribute = [star.mass, star.center_h1,
star.he_core_mass]
rs = [20.0, 1.0, 10.0]
x0 = get_root0(
MESA_label, posydon_attribute, htrack, rs=rs)
bnds = ([m_min_H, m_max_H], [0, None])
sol = minimize(
lambda x: (
sc_mass_H.transform(
get_track_val(MESA_label[0], htrack, *x))
- sc_mass_H.transform(
posydon_attribute[0])) ** 2
+ (sc_center_h1.transform(
get_track_val(MESA_label[1], htrack, *x))
- sc_center_h1.transform(
posydon_attribute[1])) ** 2
+ (sc_he_core_mass_H.transform(
get_track_val(MESA_label[2], htrack, *x))
- sc_he_core_mass_H.transform(
posydon_attribute[2])) ** 2,
x0, method="TNC", bounds=bnds
)
elif star.state in LIST_ACCEPTABLE_STATES_FOR_postMS:
MESA_label = ["he_core_mass", "mass", "center_he4", "log_R"]
posydon_attribute = [star.he_core_mass, star.mass,
star.center_he4, star.log_R]
rs = [10.0, 20.0, 1.0, 2.0]
if self.verbose:
print("Matching attributes and their normalizations : ",
MESA_label, rs)
for i in MESA_label:
if i not in self.root_keys:
raise Exception("Expected matching parameter not added"
" in the MIST model options.")
x0 = get_root0(MESA_label, posydon_attribute, htrack, rs=rs)
bnds = ([m_min_H, m_max_H], [0, None])
sol = minimize(
lambda x: (
(sc_he_core_mass_H.transform(
get_track_val(MESA_label[0], htrack, *x))
- sc_he_core_mass_H.transform(posydon_attribute[0]))
/ sc_he_core_mass_H.transform(posydon_attribute[0]))
** 2 + ((sc_mass_H.transform(
get_track_val(MESA_label[1], htrack, *x))
- sc_mass_H.transform(posydon_attribute[1]))
/ sc_mass_H.transform(posydon_attribute[1])) ** 2
+ (sc_center_he4_H.transform(
get_track_val(MESA_label[2], htrack, *x))
- sc_center_he4_H.transform(posydon_attribute[2])) ** 2
+ (sc_log_R_H.transform(
get_track_val(MESA_label[3], htrack, *x))
- sc_log_R_H.transform(posydon_attribute[3])) ** 2,
x0, method="TNC", bounds=bnds,
)
elif star.state in LIST_ACCEPTABLE_STATES_FOR_HeStar:
MESA_label = [
"he_core_mass",
# "surface_he4",
"center_he4",
# "mass",
# "center_c12",
"log_R"
]
posydon_attribute = [
star.he_core_mass,
# star.surface_he4,
star.center_he4,
# star.mass,
# star.center_c12,
star.log_R
]
# rs = [10, 2.0, 2.0, 20, 2.0]
rs = [10.0, 1.0, 2.0]
if self.verbose:
print("Matching attributes and their normalizations : ",
MESA_label, rs)
for i in MESA_label:
if i not in self.root_keys:
raise Exception("Expected matching parameter not added"
" in the MIST model options.")
x0 = get_root0(MESA_label, posydon_attribute, False, rs=rs)
bnds = ([m_min_He, m_max_He], [0, None])
# match he4 abundance using definite value
sol = minimize(
lambda x: (
(sc_he_core_mass_He.transform(
get_track_val(MESA_label[0], htrack, *x))
- sc_he_core_mass_He.transform(posydon_attribute[0]))
/ sc_he_core_mass_He.transform(posydon_attribute[0]))
** 2 + (sc_center_he4_He.transform(
get_track_val(MESA_label[1], htrack, *x))
- sc_center_he4_He.transform(posydon_attribute[1]))
** 2 + (sc_log_R_He.transform(
get_track_val(MESA_label[2], htrack, *x))
- sc_log_R_He.transform(posydon_attribute[2]))**2,
x0, method="TNC", bounds=bnds,
)
if (np.abs(sol.fun) > tolerance_mist_integration
or not sol.success):
sol = minimize(
lambda x: (
(sc_he_core_mass_He.transform(
get_track_val(MESA_label[0], htrack, *x))
- sc_he_core_mass_He.transform(
posydon_attribute[0]))
/ sc_he_core_mass_He.transform(
posydon_attribute[0])) ** 2
+ (sc_center_he4_He.transform(
get_track_val(MESA_label[1], htrack, *x))
- sc_center_he4_He.transform(posydon_attribute[1]))
** 2 + (sc_log_R_He.transform(
get_track_val(MESA_label[2], htrack, *x))
- sc_log_R_He.transform(posydon_attribute[2]))
** 2,
x0, method="Powell",
)
if (np.abs(sol.fun) > tolerance_mist_integration
or not sol.success):
star.htrack = True
x0 = get_root0(
MESA_label, posydon_attribute, star.htrack, rs=rs)
bnds = ([m_min_H, m_max_H], [0, None])
sol = minimize(
lambda x: (
(sc_he_core_mass_He.transform(
get_track_val(MESA_label[0], htrack, *x))
- sc_he_core_mass_He.transform(
posydon_attribute[0]))
/ sc_he_core_mass_He.transform(
posydon_attribute[0])) ** 2
+ (sc_center_he4_H.transform(
get_track_val(MESA_label[1], htrack, *x))
- sc_center_he4_H.transform(
posydon_attribute[1]))
** 2 + (sc_log_R_H.transform(
get_track_val(MESA_label[2], htrack, *x))
- sc_log_R_H.transform(
posydon_attribute[2])) ** 2,
x0, method="TNC", bounds=bnds,
)
if ((self.get_track_val("mass", star.htrack, *sol.x)
- self.get_track_val(
"he_core_mass", star.htrack, *sol.x))
/ self.get_track_val(
"mass", star.htrack, *sol.x) >= 0.05):
initials = (np.nan, np.nan)
if initials is None:
if self.verbose:
print("sol.fun = ", np.abs(sol.fun))
if np.abs(sol.fun) < tolerance_mist_integration:
if self.verbose:
print("minimization in matching considered acceptable,"
" with", np.abs(sol.fun), "<",
tolerance_mist_integration, "tolerance")
initials = sol.x
star.fun = sol.fun
star.stiching_rel_mass_difference = (
self.get_track_val("mass", htrack, *sol.x) - star.mass
) / star.mass
if star.log_R in MESA_label:
star.stiching_rel_logRadius_difference = (
self.get_track_val("log_R", htrack, *sol.x)
- star.log_R) / star.log_R
else:
star.stiching_rel_logRadius_difference = np.nan
if (star.total_moment_of_inertia is not None
and not np.isnan(star.total_moment_of_inertia)):
star.stiching_rel_inertia_difference = (
self.get_track_val("inertia", htrack, *sol.x)
- star.total_moment_of_inertia
) / star.total_moment_of_inertia
else:
if self.verbose:
print("minimization in matching not successful, with",
np.abs(sol.fun), ">", tolerance_mist_integration,
"tolerance")
initials = (np.nan, np.nan)
star.fun = np.nan
star.stiching_rel_mass_difference = np.nan
star.stiching_rel_radius_difference = np.nan
star.stiching_rel_inertia_difference = np.nan
if self.verbose:
print(
"matching with track of intial mass m0, at time t0 ",
initials[0],
initials[1],
"with m(t0), log10(R(t0), center_he(t0), surface_he4, "
"surface_h1, he_core_mass, center_c12) = ",
self.get_track_val("mass", htrack, *sol.x),
self.get_track_val("log_R", htrack, *sol.x),
self.get_track_val("center_he4", htrack, *sol.x),
self.get_track_val("surface_he4", htrack, *sol.x),
self.get_track_val("surface_h1", htrack, *sol.x),
self.get_track_val("he_core_mass", htrack, *sol.x),
self.get_track_val("center_c12", htrack, *sol.x),
"Mass / Radius of the secondary at the end of the previous "
"step was = ",
star.mass,
star.log_R,
star.center_he4,
star.surface_he4,
star.surface_h1,
star.he_core_mass,
star.center_c12
)
return initials
def __repr__(self):
"""Return the type of evolution type."""
return "Detached Step."
def __call__(self, binary):
"""Evolve the binary until RLO or compact object formation."""
get_mist0 = self.get_mist0
KEYS = self.KEYS
KEYS_POSITIVE = self.KEYS_POSITIVE
# the primary is the more evolved star
if (binary.star_1.state in ("BH", "NS", "WD")
and binary.star_2.state in STAR_STATES_H_RICH):
primary = binary.star_1
secondary = binary.star_2
secondary.htrack = True
primary.htrack = secondary.htrack
primary.co = True
elif (binary.star_1.state in ("BH", "NS", "WD")
and binary.star_2.state in LIST_ACCEPTABLE_STATES_FOR_HeStar):
primary = binary.star_1
secondary = binary.star_2
secondary.htrack = False
primary.htrack = secondary.htrack
primary.co = True
elif (binary.star_2.state in ("BH", "NS", "WD")
and binary.star_1.state in STAR_STATES_H_RICH):
primary = binary.star_2
secondary = binary.star_1
secondary.htrack = True
primary.htrack = secondary.htrack
primary.co = True
elif (binary.star_2.state in ("BH", "NS", "WD")
and binary.star_1.state in LIST_ACCEPTABLE_STATES_FOR_HeStar):
primary = binary.star_2
secondary = binary.star_1
secondary.htrack = False
primary.htrack = secondary.htrack
primary.co = True
elif (binary.star_1.state in STAR_STATES_H_RICH
and binary.star_2.state in STAR_STATES_H_RICH):
primary = binary.star_1
secondary = binary.star_2
secondary.htrack = True
primary.htrack = True
primary.co = False
elif (binary.star_1.state in LIST_ACCEPTABLE_STATES_FOR_HeStar
and binary.star_2.state in STAR_STATES_H_RICH):
primary = binary.star_1
secondary = binary.star_2
secondary.htrack = True
primary.htrack = False
primary.co = False
elif (binary.star_2.state in LIST_ACCEPTABLE_STATES_FOR_HeStar
and binary.star_1.state in STAR_STATES_H_RICH):
primary = binary.star_2
secondary = binary.star_1
secondary.htrack = True
primary.htrack = False
primary.co = False
elif (binary.star_1.state in LIST_ACCEPTABLE_STATES_FOR_HeStar
and binary.star_2.state in LIST_ACCEPTABLE_STATES_FOR_HeStar):
primary = binary.star_1
secondary = binary.star_2
secondary.htrack = False
primary.htrack = False
primary.co = False
else:
raise Exception("States not recognized!")
if not primary.co:
m1, t1 = get_mist0(primary, primary.htrack)
m2, t2 = get_mist0(secondary, secondary.htrack)
def get_star_data(binary, star1, star2, htrack, co):
"""Get and interpolate the properties of stars.
The data of a compact object can be stored as a copy of its
companion for convenience except its mass, radius, mdot, Idot,
and radius, mdot, Idot are set to be zero.
Parameters
----------
htrack : bool
htrack of star1
co: bool
co of star2
Return
-------
interp1d
Contains the properties of star1 if co is false,
if co is true, star2 is a compact object,
return the properties of star2
"""
if htrack:
self.grid = self.grid1
elif not htrack:
self.grid = self.grid2
get_track = self.grid.get
with np.errstate(all="ignore"):
# get the initial m0, t0 track
m0, t0 = get_mist0(star1, htrack)
if np.any(np.isnan([m0, t0])):
# binary.event = "END"
# binary.state += " (GridMatchingFailed)"
# if self.verbose:
# print("Failed matching")
return None
max_time = binary.properties.max_simulation_time
assert max_time > 0.0, "max_time is non-positive"
age = get_track("age", m0)
t_max = age.max() # max timelength of the track
interp1d = dict()
kvalue = dict()
for key in KEYS[1:]:
kvalue[key] = get_track(key, m0)
try:
for key in KEYS[1:]:
if key in KEYS_POSITIVE:
positive = True
interp1d[key] = PchipInterpolator2(age, kvalue[key],
positive=positive)
else:
interp1d[key] = PchipInterpolator2(age, kvalue[key])
except ValueError:
i_bad = [None]
while len(i_bad) != 0:
i_bad = np.where(np.diff(age) <= 0)[0]
age = np.delete(age, i_bad)
for key in KEYS[1:]:
kvalue[key] = np.delete(kvalue[key], i_bad)
for key in KEYS[1:]:
if key in KEYS_POSITIVE:
positive = True
interp1d[key] = PchipInterpolator2(age, kvalue[key],
positive=positive)
else:
interp1d[key] = PchipInterpolator2(age, kvalue[key])
interp1d["inertia"] = PchipInterpolator(
age, kvalue["inertia"] / (const.msol * const.rsol**2))
interp1d["Idot"] = interp1d["inertia"].derivative()
interp1d["L"] = PchipInterpolator(age, 10 ** kvalue["log_L"])
interp1d["R"] = PchipInterpolator(age, 10 ** kvalue["log_R"])
interp1d["t_max"] = t_max
interp1d["max_time"] = max_time
interp1d["t0"] = t0
interp1d["m0"] = m0
if co:
kvalue["mass"] = np.zeros_like(kvalue["mass"]) + star2.mass
kvalue["R"] = np.zeros_like(kvalue["log_R"])
kvalue["mdot"] = np.zeros_like(kvalue["mdot"])
interp1d["mass"] = PchipInterpolator(age, kvalue["mass"])
interp1d["R"] = PchipInterpolator(age, kvalue["R"])
interp1d["mdot"] = PchipInterpolator(age, kvalue["mdot"])
interp1d["Idot"] = PchipInterpolator(age, kvalue["mdot"])
return interp1d
# get the matched data of two stars, respectively
interp1d_sec = get_star_data(
binary, secondary, primary, secondary.htrack, False)
if primary.co:
interp1d_pri = get_star_data(
binary, secondary, primary, secondary.htrack, True)
elif not primary.co:
interp1d_pri = get_star_data(
binary, primary, secondary, primary.htrack, False)
if interp1d_sec is None or interp1d_pri is None:
# binary.event = "END"
binary.state += " (GridMatchingFailed)"
if self.verbose:
print("Failed matching")
return
t0_sec = interp1d_sec["t0"]
t0_pri = interp1d_pri["t0"]
m01 = interp1d_sec["m0"]
m02 = interp1d_pri["m0"]
t_max_sec = interp1d_sec["t_max"]
t_max_pri = interp1d_pri["t_max"]
t_offset_sec = binary.time - t0_sec
t_offset_pri = binary.time - t0_pri
max_time = interp1d_sec["max_time"]
@event(True, 1)
def ev_rlo1(t, y):
"""Difference between radius and Roche lobe at a given time.
Used to check if there is RLOF mass transfer during the detached
binary evolution interpolation.
Parameters
----------
t : float
Time of the evolution, in years.
y : tuple of floats
[separation, eccentricity] at that time. Separation should be
in solar radii.
Returns
-------
float
Difference between stellar radius and Roche lobe radius in
solar radii.
"""
sep = y[0]
ecc = y[1]
RL = roche_lobe_radius(interp1d_sec["mass"](t - t_offset_sec)
/ interp1d_pri["mass"](t - t_offset_pri),
(1 - ecc) * sep)
# 95% filling of the RL is enough to assume beginning of RLO,
# as we do in CO-HMS_RLO grid
return interp1d_sec["R"](t - t_offset_sec) - 0.95*RL
@event(True, 1)
def ev_rlo2(t, y):
"""Difference between radius and Roche lobe at a given time.
Used to check if there is RLOF mass transfer during the detached
binary evolution interpolation.
Parameters
----------
t : float
Time of the evolution, in years
y : tuple of floats
[separation, eccentricity] at that time. Separation should be
in solar radii.
Returns
-------
float
Difference between stellar radius and Roche lobe radius in
solar radii.
"""
sep = y[0]
ecc = y[1]
RL = roche_lobe_radius(interp1d_pri["mass"](t - t_offset_pri)
/ interp1d_sec["mass"](t - t_offset_sec),
(1 - ecc) * sep)
return interp1d_pri["R"](t - t_offset_pri) - 0.95*RL
@event(True, 1)
def ev_rel_rlo1(t, y):
"""Relative difference between radius and Roche lobe.
Used to check if there is RLOF mass transfer during the detached
binary evolution interpolation.
Parameters
----------
t : float
Time of the evolution, in years.
y : tuple of floats
[separation, eccentricity] at that time. Separation should be
in solar radii.
Returns
-------
float
Relative difference between stellar radius and Roche lobe
radius.
"""
sep = y[0]
ecc = y[1]
RL = roche_lobe_radius(interp1d_sec["mass"](t - t_offset_sec)
/ interp1d_pri["mass"](t - t_offset_pri),
(1 - ecc) * sep)
return (interp1d_sec["R"](t - t_offset_sec) - RL) / RL
@event(True, 1)
def ev_rel_rlo2(t, y):
"""Relative difference between radius and Roche lobe.
Used to check if there is RLOF mass transfer during the detached
binary evolution interpolation.
Parameters
----------
t : float
Time of the evolution, in years.
y : tuple of floats
[separation, eccentricity] at that time. Separation should be
in solar radii.
Returns
-------
float
Relative difference between stellar radius and Roche lobe
radius.
"""
sep = y[0]
ecc = y[1]
RL = roche_lobe_radius(interp1d_pri["mass"](t - t_offset_pri)
/ interp1d_sec["mass"](t - t_offset_sec),
(1 - ecc) * sep)
return (interp1d_pri["R"](t - t_offset_pri) - RL) / RL
@event(True, -1)
def ev_max_time1(t, y):
return t_max_sec + t_offset_sec - t
@event(True, -1)
def ev_max_time2(t, y):
return t_max_pri + t_offset_pri - t
# make a function to get the spin of two stars
def get_omega(star):
if (star.log_total_angular_momentum is not None
and star.total_moment_of_inertia is not None
and not np.isnan(star.log_total_angular_momentum)
and not np.isnan(star.total_moment_of_inertia)):
omega_in_rad_per_year = (
10.0 ** star.log_total_angular_momentum
/ star.total_moment_of_inertia * const.secyer
) # the last factor transforms it from rad/s to rad/yr
if self.verbose:
print("calculating initial omega from angular momentum and"
" moment of inertia", omega_in_rad_per_year)
# except TypeError:
else:
# we equate secondary's initial omega to surf_avg_omega
# (although the critical rotation should be improved to
# take into accoun radiation pressure)
if (star.surf_avg_omega is not None
and not np.isnan(star.surf_avg_omega)):
omega_in_rad_per_year = star.surf_avg_omega * const.secyer
# the last factor transforms it from rad/s to rad/yr.
if self.verbose:
print("calculating initial omega from surf_avg_omega",
omega_in_rad_per_year)
elif (star.surf_avg_omega_div_omega_crit is not None
and not np.isnan(star.surf_avg_omega_div_omega_crit)):
omega_in_rad_per_year = (
star.surf_avg_omega_div_omega_crit * np.sqrt(
const.standard_cgrav * star.mass * const.msol
/ ((10.0 ** (star.log_R) * const.rsol) ** 3))
* const.secyer)
# the last factor transforms it from rad/s to rad/yr
# EDIT: We assume POSYDON surf_avg_omega is provided in
# rad/yr already.
if self.verbose:
print("calculating initial omega from "
"surf_avg_omega_div_omega_crit",
omega_in_rad_per_year)
else:
omega_in_rad_per_year = 0.0
if self.verbose:
print("could calculate initial omega",
omega_in_rad_per_year)
if self.verbose:
print("initial omega_in_rad_per_year", omega_in_rad_per_year)
return omega_in_rad_per_year
if (ev_rlo1(binary.time, [binary.separation, binary.eccentricity]) >= 0
or ev_rlo2(binary.time,
[binary.separation, binary.eccentricity])
>= 0):
binary.state = "initial_RLOF"
return
# binary.event = "END"
# TODO: put it out of its misery here!
else:
if not (max_time - binary.time > 0.0):
raise Exception("max_time is lower than the current time. "
"Evolution of the detached binary will go to "
"lower times.")
with np.errstate(all="ignore"):
omega_in_rad_per_year_sec = get_omega(secondary)
if primary.co:
# omega of compact objects won't be used for intergration
omega_in_rad_per_year_pri = omega_in_rad_per_year_sec
elif not primary.co:
omega_in_rad_per_year_pri = get_omega(primary)
try:
s = solve_ivp(
lambda t, y: diffeq(
t,
y,
interp1d_sec["R"](t - t_offset_sec),
interp1d_sec["L"](t - t_offset_sec),
*[
interp1d_sec[key](t - t_offset_sec)
for key in KEYS[1:11]
],
interp1d_sec["Idot"](t - t_offset_sec),
interp1d_pri["R"](t - t_offset_pri),
interp1d_pri["L"](t - t_offset_pri),
*[
interp1d_pri[key](t - t_offset_pri)
for key in KEYS[1:11]
],
interp1d_pri["Idot"](t - t_offset_pri),
self.do_wind_loss,
self.do_tides,
self.do_gravitational_radiation,
self.do_magnetic_braking,
self.do_stellar_evolution_and_spin_from_winds
# ,self.verbose
),
events=[ev_rlo1, ev_rlo2, ev_max_time1, ev_max_time2],
method="Radau",
t_span=(binary.time, max_time),
y0=[
binary.separation,
binary.eccentricity,
omega_in_rad_per_year_sec,
omega_in_rad_per_year_pri,
],
dense_output=True,
# vectorized=True
)
except Exception:
s = solve_ivp(
lambda t, y: diffeq(
t,
y,
interp1d_sec["R"](t - t_offset_sec),
interp1d_sec["L"](t - t_offset_sec),
*[
interp1d_sec[key](t - t_offset_sec)
for key in KEYS[1:11]
],
interp1d_sec["Idot"](t - t_offset_sec),
interp1d_pri["R"](t - t_offset_pri),
interp1d_pri["L"](t - t_offset_pri),
*[
interp1d_pri[key](t - t_offset_pri)
for key in KEYS[1:11]
],
interp1d_pri["Idot"](t - t_offset_pri),
self.do_wind_loss,
self.do_tides,
self.do_gravitational_radiation,
self.do_magnetic_braking,
self.do_stellar_evolution_and_spin_from_winds
# ,self.verbose
),
events=[ev_rlo1, ev_rlo2, ev_max_time1, ev_max_time2],
method="RK45",
t_span=(binary.time, max_time),
y0=[
binary.separation,
binary.eccentricity,
omega_in_rad_per_year_sec,
omega_in_rad_per_year_pri,
],
dense_output=True,
# vectorized=True
)
if self.verbose:
print("solution of ODE", s)
if s.status == -1:
print("Integration failed", s.message)
binary.state += ' (Integration failure)'
# binary.event = "END"
return
# raise RuntimeError("Integration failed", s.message)
if self.dt is not None and self.dt > 0:
t = np.arange(binary.time, s.t[-1] + self.dt/2.0, self.dt)[1:]
if t[-1] < s.t[-1]:
t = np.hstack([t, s.t[-1]])
elif (self.n_o_steps_history is not None
and self.n_o_steps_history > 0):
t_step = (s.t[-1] - binary.time) / self.n_o_steps_history
t = np.arange(binary.time, s.t[-1] + t_step / 2.0, t_step)[1:]
if t[-1] < s.t[-1]:
t = np.hstack([t, s.t[-1]])
else: # self.dt is None and self.n_o_steps_history is None
t = np.array([s.t[-1]])
orb_params = s.sol(t)
sep_interp, ecc_interp, omega_interp_sec, omega_interp_pri = s.sol(
t)
mass_interp_sec = interp1d_sec[self.translate["mass"]]
mass_interp_pri = interp1d_pri[self.translate["mass"]]
# s.sol(t)[0]
# interp1d = dict()
interp1d_sec["sep"] = s.sol(t)[0]
# lambda x: orb_params[0][np.argmax(x == t[:, None], 0)]
interp1d_sec["ecc"] = s.sol(t)[1]
# lambda x: orb_params[1][np.argmax(x == t[:, None], 0)]
interp1d_sec["omega"] = s.sol(t)[2]
# lambda x: orb_params[2][np.argmax(x == t[:, None], 0)]
interp1d_pri["omega"] = s.sol(t)[3]
interp1d_sec["porb"] = orbital_period_from_separation(
sep_interp, mass_interp_sec(t - t_offset_sec),
mass_interp_pri(t - t_offset_pri))
interp1d_pri["porb"] = orbital_period_from_separation(
sep_interp, mass_interp_pri(t - t_offset_pri),
mass_interp_sec(t - t_offset_sec))
interp1d_sec["time"] = t # binary.time + x - t0
# time_interp = binary.time + t - t0
for obj, prop in zip(
[secondary, primary, binary],
[STARPROPERTIES, STARPROPERTIES, BINARYPROPERTIES],
):
for key in prop:
if key in ["event",
"mass_transfer_case",
"nearest_neighbour_distance",
"state", "metallicity", "V_sys"]:
current = getattr(obj, key)
# For star objects, the state is calculated
# further below
history = [current] * len(t[:-1])
# elif key in ("state") and obj == secondary:
# current = check_state_of_star(obj, star_CO=False)
# history = [current] * len(t[:-1])
elif (key in ["surf_avg_omega_div_omega_crit"]
and obj == secondary):
# In fact I replace the actual surf_avg_w with the ef-
# fective omega which takes into account the whole star
# key = 'effective_omega' # in rad/sec
# current = s.y[2][-1] / 3.1558149984e7
# history_of_attribute = s.y[2][:-1] / 3.1558149984e7
omega_crit_current_sec = np.sqrt(
const.standard_cgrav
* interp1d_sec[self.translate["mass"]](
t[-1] - t_offset_sec).item() * const.msol
/ (interp1d_sec[self.translate["R"]](
t[-1] - t_offset_sec).item() * const.rsol) ** 3
)
omega_crit_hist_sec = np.sqrt(
const.standard_cgrav
* interp1d_sec[self.translate["mass"]](
t[:-1] - t_offset_sec) * const.msol
/ (interp1d_sec[self.translate["R"]](
t[:-1] - t_offset_sec) * const.rsol) ** 3
)
current = (interp1d_sec["omega"][-1] / const.secyer
/ omega_crit_current_sec)
history = (interp1d_sec["omega"][:-1] / const.secyer
/ omega_crit_hist_sec)
elif (key in ["surf_avg_omega_div_omega_crit"]
and obj == primary):
if primary.co:
current = None
history = [current] * len(t[:-1])
elif not primary.co:
# TODO: change `item()` to 0
omega_crit_current_pri = np.sqrt(
const.standard_cgrav
* interp1d_pri[self.translate["mass"]](
t[-1] - t_offset_pri).item() * const.msol
/ (interp1d_pri[self.translate["R"]](
t[-1] - t_offset_pri).item() * const.rsol)
** 3)
omega_crit_hist_pri = np.sqrt(
const.standard_cgrav
* interp1d_pri[self.translate["mass"]](
t[:-1] - t_offset_pri) * const.msol
/ (interp1d_pri[self.translate["R"]](
t[:-1] - t_offset_pri) * const.rsol) ** 3)
current = (interp1d_pri["omega"][-1]
/ const.secyer / omega_crit_current_pri)
history = (interp1d_pri["omega"][:-1]
/ const.secyer / omega_crit_hist_pri)
elif key in ["surf_avg_omega"] and obj == secondary:
# current = interp1d["omega"](t[-1]) / const.secyer
current = interp1d_sec["omega"][-1] / const.secyer
history = interp1d_sec["omega"][:-1] / const.secyer
elif key in ["surf_avg_omega"] and obj == primary:
if primary.co:
current = None
history = [current] * len(t[:-1])
else:
current = interp1d_pri["omega"][-1] / const.secyer
history = interp1d_pri["omega"][:-1] / const.secyer
elif key in ["rl_relative_overflow_1"] and obj == binary:
if binary.star_1.state in ("BH", "NS", "WD"):
current = None
history = [current] * len(t[:-1])
elif secondary == binary.star_1:
current = ev_rel_rlo1(t[-1],
[interp1d_sec["sep"][-1],
interp1d_sec["ecc"][-1]])
history = ev_rel_rlo1(t[:-1],
[interp1d_sec["sep"][:-1],
interp1d_sec["ecc"][:-1]])
elif secondary == binary.star_2:
current = ev_rel_rlo2(t[-1],
[interp1d_sec["sep"][-1],
interp1d_sec["ecc"][-1]])
history = ev_rel_rlo2(t[:-1],
[interp1d_sec["sep"][:-1],
interp1d_sec["ecc"][:-1]])
elif key in ["rl_relative_overflow_2"] and obj == binary:
if binary.star_2.state in ("BH", "NS", "WD"):
current = None
history = [current] * len(t[:-1])
elif secondary == binary.star_2:
current = ev_rel_rlo1(t[-1],
[interp1d_sec["sep"][-1],
interp1d_sec["ecc"][-1]])
history = ev_rel_rlo1(t[:-1],
[interp1d_sec["sep"][:-1],
interp1d_sec["ecc"][:-1]])
elif secondary == binary.star_1:
current = ev_rel_rlo2(t[-1],
[interp1d_sec["sep"][-1],
interp1d_sec["ecc"][-1]])
history = ev_rel_rlo2(t[:-1],
[interp1d_sec["sep"][:-1],
interp1d_sec["ecc"][:-1]])
elif key in ["separation", "orbital_period",
"eccentricity", "time"]:
current = interp1d_sec[self.translate[key]][-1].item()
history = interp1d_sec[self.translate[key]][:-1]
elif (key in ["total_moment_of_inertia"]
and obj == secondary):
current = interp1d_sec[self.translate[key]](
t[-1] - t_offset_sec).item() * (
const.msol * const.rsol ** 2)
history = interp1d_sec[self.translate[key]](
t[:-1] - t_offset_sec) * (
const.msol * const.rsol ** 2)
elif key in ["total_moment_of_inertia"] and obj == primary:
if primary.co:
current = getattr(obj, key)
history = [current] * len(t[:-1])
else:
current = interp1d_pri[self.translate[key]](
t[-1] - t_offset_pri).item() * (
const.msol * const.rsol ** 2)
history = interp1d_pri[self.translate[key]](
t[:-1] - t_offset_pri) * (
const.msol * const.rsol ** 2)
elif (key in ["log_total_angular_momentum"]
and obj == secondary):
current = np.log10(
(interp1d_sec["omega"][-1] / const.secyer)
* (interp1d_sec[
self.translate["total_moment_of_inertia"]](
t[-1] - t_offset_sec).item() * (
const.msol * const.rsol ** 2)))
history = np.log10(
(interp1d_sec["omega"][:-1] / const.secyer)
* (interp1d_sec[
self.translate["total_moment_of_inertia"]](
t[:-1] - t_offset_sec) * (
const.msol * const.rsol ** 2)))
elif (key in ["log_total_angular_momentum"]
and obj == primary):
if primary.co:
current = getattr(obj, key)
history = [current] * len(t[:-1])
else:
current = np.log10(
(interp1d_pri["omega"][-1] / const.secyer)
* (interp1d_pri[
self.translate["total_moment_of_inertia"]](
t[-1] - t_offset_pri).item() * (
const.msol * const.rsol ** 2)))
history = np.log10(
(interp1d_pri["omega"][:-1] / const.secyer)
* (interp1d_pri[
self.translate["total_moment_of_inertia"]](
t[:-1] - t_offset_pri) * (
const.msol * const.rsol ** 2)))
elif key in ["spin"] and obj == secondary:
current = (
const.clight
* (interp1d_sec["omega"][-1] / const.secyer)
* interp1d_sec[
self.translate["total_moment_of_inertia"]](
t[-1] - t_offset_sec).item()
* (const.msol * const.rsol ** 2)
/ (const.standard_cgrav * (
interp1d_sec[self.translate["mass"]](
t[-1] - t_offset_sec).item()
* const.msol) ** 2))
history = (
const.clight
* (interp1d_sec["omega"][:-1] / const.secyer)
* interp1d_sec[
self.translate["total_moment_of_inertia"]](
t[:-1] - t_offset_sec)
* (const.msol * const.rsol ** 2)
/ (const.standard_cgrav * (
interp1d_sec[self.translate["mass"]](
t[:-1] - t_offset_sec) * const.msol)**2))
elif key in ["spin"] and obj == primary:
if primary.co:
current = getattr(obj, key)
history = [current] * len(t[:-1])
else:
current = (
const.clight
* (interp1d_pri["omega"][-1] / const.secyer)
* interp1d_pri[
self.translate["total_moment_of_inertia"]](
t[-1] - t_offset_pri).item()
* (const.msol * const.rsol ** 2)
/ (const.standard_cgrav * (
interp1d_pri[self.translate["mass"]](
t[-1] - t_offset_pri).item()
* const.msol)**2))
history = (
const.clight * (interp1d_pri["omega"][:-1]
/ const.secyer)
* interp1d_pri[
self.translate["total_moment_of_inertia"]](
t[:-1] - t_offset_pri)
* (const.msol * const.rsol ** 2)
/ (const.standard_cgrav * (interp1d_pri[
self.translate["mass"]](
t[:-1] - t_offset_pri)
* const.msol)**2))
elif (key in ["lg_mdot", "lg_wind_mdot"]
and obj == secondary):
# in detached step, lg_mdot = lg_wind_mdot
if interp1d_sec[self.translate[key]](
t[-1] - t_offset_sec) == 0:
current = -98.99
else:
current = np.log10(
np.abs(interp1d_sec[self.translate[key]](
t[-1] - t_offset_sec))).item()
history = np.ones_like(t[:-1])
for i in range(len(t)-1):
if interp1d_sec[self.translate[key]](
t[i] - t_offset_sec) == 0:
history[i] = -98.99
else:
history[i] = np.log10(
np.abs(interp1d_sec[self.translate[key]](
t[i] - t_offset_sec)))
elif key in ["lg_mdot", "lg_wind_mdot"] and obj == primary:
if primary.co:
current = None
history = [current] * len(t[:-1])
else:
if interp1d_sec[self.translate[key]](
t[-1] - t_offset_sec) == 0:
current = -98.99
else:
current = np.log10(np.abs(
interp1d_sec[self.translate[key]](
t[-1] - t_offset_sec))).item()
history = np.ones_like(t[:-1])
for i in range(len(t)-1):
if (interp1d_sec[self.translate[key]](
t[i] - t_offset_sec) == 0):
history[i] = -98.99
else:
history[i] = np.log10(np.abs(
interp1d_sec[self.translate[key]](
t[i] - t_offset_sec)))
elif (self.translate[key] in interp1d_sec
and obj == secondary):
current = interp1d_sec[self.translate[key]](
t[-1] - t_offset_sec).item()
history = interp1d_sec[self.translate[key]](
t[:-1] - t_offset_sec)
elif (self.translate[key] in interp1d_pri
and obj == primary):
if primary.co:
current = getattr(obj, key)
history = [current] * len(t[:-1])
else:
current = interp1d_pri[self.translate[key]](
t[-1] - t_offset_pri).item()
history = interp1d_pri[self.translate[key]](
t[:-1] - t_offset_pri)
elif key in ["profile"]:
current = None
history = [current] * len(t[:-1])
else:
current = np.nan
history = np.ones_like(t[:-1]) * current
setattr(obj, key, current)
getattr(obj, key + "_history").extend(history)
secondary.state = check_state_of_star(secondary, star_CO=False)
for timestep in range(-len(t[:-1]), 0):
secondary.state_history[timestep] = check_state_of_star(
secondary, i=timestep, star_CO=False)
if primary.co:
mdot_acc = np.atleast_1d(bondi_hoyle(
binary, primary, secondary, slice(-len(t), None),
wind_disk_criteria=True, scheme='Kudritzki+2000'))
primary.lg_mdot = np.log10(mdot_acc.item(-1))
primary.lg_mdot_history[len(primary.lg_mdot_history) - len(t)
+ 1:] = np.log10(mdot_acc[:-1])
elif not primary.co:
primary.state = check_state_of_star(primary, star_CO=False)
for timestep in range(-len(t[:-1]), 0):
primary.state_history[timestep] = check_state_of_star(
primary, i=timestep, star_CO=False)
def get_star_final_values(star, htrack, m0):
if htrack:
self.grid = self.grid1
elif not htrack:
self.grid = self.grid2
get_final_values = self.grid.get_final_values
get_final_state = self.grid.get_final_state
for key in self.final_keys:
setattr(star, key, get_final_values('S1_%s' % (key), m0))
def get_star_profile(star, htrack, m0):
if htrack:
self.grid = self.grid1
elif not htrack:
self.grid = self.grid2
get_profile = self.grid.get_profile
profile_new = np.array(get_profile('mass', m0)[1])
for i in self.profile_keys:
profile_new[i] = get_profile(i, m0)[0]
profile_new['omega'] = star.surf_avg_omega
star.profile = profile_new
if s.t_events[0] or s.t_events[1]: # reached RLOF
if self.RLO_orbit_at_orbit_with_same_am:
# final circular orbit conserves angular momentum
# compared to the eccentric orbit
binary.separation *= (1 - s.y[1][-1]**2)
binary.orbital_period *= (1 - s.y[1][-1]**2) ** 1.5
else:
# final circular orbit is at periastron of the ecc. orbit
binary.separation *= (1 - s.y[1][-1])
binary.orbital_period *= (1 - s.y[1][-1]) ** 1.5
assert np.abs(
binary.orbital_period
- orbital_period_from_separation(
binary.separation, secondary.mass, primary.mass
)
) / binary.orbital_period < 10 ** (-2)
binary.eccentricity = 0
if s.t_events[0]:
if secondary == binary.star_1:
binary.state = "RLO1"
binary.event = "oRLO1"
else:
binary.state = "RLO2"
binary.event = "oRLO2"
elif s.t_events[1]:
if secondary == binary.star_1:
binary.state = "RLO2"
binary.event = "oRLO2"
else:
binary.state = "RLO1"
binary.event = "oRLO1"
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"
else:
binary.event = "CC2"
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:
# simultaneous core-collapse of the other star as well
primary_time = t_max_pri + t_offset_pri - t[-1]
secondary_time = t_max_sec + t_offset_sec - t[-1]
if primary_time == secondary_time:
# we manually check if s.t_events[3] should also
# be happening simultaneously
get_star_final_values(primary, primary.htrack, m02)
get_star_profile(primary, primary.htrack, m02)
if primary.mass != secondary.mass:
raise ValueError(
"Both stars are found to be ready for collapse "
"(i.e. end of their life) during the detached "
"step, but do not have the same mass")
elif s.t_events[3]:
# reached t_max of track. End of life (possible collapse) of
# primary
if secondary == binary.star_1:
binary.event = "CC2"
else:
binary.event = "CC1"
get_star_final_values(primary, primary.htrack, m02)
get_star_profile(primary, primary.htrack, m02)
else: # Reached max_time asked.
if binary.properties.max_simulation_time - binary.time < 0.0:
binary.event = "MaxTime_exceeded"
else:
binary.event = "maxtime"
# binary.event = "MaxTime_exceeded"
[docs]def event(terminal, direction=0):
"""Return a helper function to set attributes for solve_ivp events."""
def dec(f):
f.terminal = True
f.direction = direction
return f
return dec
[docs]def diffeq(
t,
y,
R_sec,
L_sec,
M_sec,
Mdot_sec,
I_sec,
# he_core_mass,
# mass_conv_core,
# conv_mx1_top,
# conv_mx1_bot,
conv_mx1_top_r_sec,
conv_mx1_bot_r_sec,
surface_h1_sec,
center_h1_sec,
M_env_sec,
DR_env_sec,
Renv_middle_sec,
Idot_sec,
R_pri,
L_pri,
M_pri,
Mdot_pri,
I_pri,
# he_core_mass,
# mass_conv_core,
# conv_mx1_top,
# conv_mx1_bot,
conv_mx1_top_r_pri,
conv_mx1_bot_r_pri,
surface_h1_pri,
center_h1_pri,
M_env_pri,
DR_env_pri,
Renv_middle_pri,
Idot_pri,
do_wind_loss=True,
do_tides=True,
do_gravitational_radiation=True,
do_magnetic_braking=True,
do_stellar_evolution_and_spin_from_winds=True,
verbose=False,
):
"""Diff. equation describing the orbital evolution of a detached binary.
The equation handles wind mass-loss [1_], tidal [2_], gravational [3_]
effects and magnetic braking [4_]. It also handles the change of the
secondary's stellar spin due to its change of moment of intertia and due to
mass-loss from its spinning surface. It is assumed that the mass loss is
fully non-conservative. Magnetic braking is fully applied to secondary
stars with mass less than 1.3 Msun and fully off for stars with mass larger
then 1.5 Msun. The effect of magnetic braking falls linearly for stars with
mass between 1.3 Msun and 1.5 Msun.
TODO: exaplin new features (e.g., double COs)
Parameters
----------
t : float
The age of the system in years
y : list of float
Contains the separation, eccentricity and angular velocity, in Rsolar,
dimensionless and rad/year units, respectively.
M_pri : float
Mass of the primary in Msolar units.
M_sec : float
Mass of the secondary in Msolar units.
Mdot : float
Rate of change of mass of the star in Msolar/year units.
(Negative for wind mass loss.)
R : float
Radius of the star in Rsolar units.
I : float
Moment of inertia of the star in Msolar*Rsolar^2.
L : float
Luminosity of the star in solar units.
#mass_conv_core : float
# Convective core mass of the secondary in Msolar units.
conv_mx1_top_r : float
Coordinate of top convective mixing zone coordinate in Rsolar.
conv_mx1_bot_r : float
Coordinate of bottom convective mixing zone coordinate in Rsolar.
surface_h1 : float
surface mass Hydrogen abundance
center_h1 : float
center mass Hydrogen abundance
M_env : float
mass of the dominant convective region for tides above the core,
in Msolar.
DR_env : float
thickness of the dominant convective region for tides above the core,
in Rsolar.
Renv_middle : float
position of the dominant convective region for tides above the core,
in Rsolar.
Idot : float
Rate of change of the moment of inertia of the star in
Msolar*Rsolar^2 per year.
do_wind_loss: Boolean
If True, take into account change of separation due to mass loss from
the secondary. Default: True.
do_tides: Booleans
If True, take into account change of separation, eccentricity and
secondary spin due to tidal forces. Default: True.
do_gravitational_radiation: Boolean
If True, take into account change of separation and eccentricity due to
gravitational wave radiation. Default: True
do_magnetic_braking: Boolean
If True, take into account change of star spin due to magnetic braking.
Default: True.
do_stellar_evolution_and_spin_from_winds: Boolean
If True, take into account change of star spin due to change of its
moment of inertia during its evolution and due to spin angular momentum
loss due to winds. Default: True.
verbose : Boolean
If we want to print stuff. Default: False.
Returns
-------
list of float
Contains the change of
the separation, eccentricity and angular velocity, in Rsolar,
dimensionless and rad/year units, respectively.
References
----------
.. [1] Tauris, T. M., & van den Heuvel, E. 2006,
Compact stellar X-ray sources, 1, 623
.. [2] Hut, P. 1981, A&A, 99, 126
.. [3] Junker, W., & Schafer, G. 1992, MNRAS, 254, 146
.. [4] Rappaport, S., Joss, P. C., & Verbunt, F. 1983, ApJ, 275, 713
"""
y[0] = np.max([y[0], 0]) # We limit separation to non-negative values
a = y[0]
y[1] = np.max([y[1], 0]) # We limit eccentricity to non-negative values
e = y[1]
if e > 0 and e < 10.0 ** (-3):
# we force a negligible eccentricity to become 0
# for computational stability
e = 0.0
if verbose:
print("negligible eccentricity became 0 for "
"computational stability")
y[2] = np.max([y[2], 0]) # We limit omega spin to non-negative values
Omega_sec = y[2] # in rad/yr
y[3] = np.max([y[3], 0])
Omega_pri = y[3]
da = 0.0
de = 0.0
dOmega_sec = 0.0
dOmega_pri = 0.0
# Mass Loss
if do_wind_loss:
q1 = M_sec / M_pri
k11 = (1 / (1 + q1)) * (Mdot_sec / M_sec)
k21 = Mdot_sec / M_sec
k31 = Mdot_sec / (M_pri + M_sec)
# This is simplified to da_mt = -a * Mdot/(M+Macc), for only (negative)
# wind Mdot from star M.
da_mt_sec = a * (2 * k11 - 2 * k21 + k31)
q2 = M_pri / M_sec
k12 = (1 / (1 + q2)) * (Mdot_pri / M_pri)
k22 = Mdot_pri / M_pri
k32 = Mdot_pri / (M_pri + M_sec)
da_mt_pri = a * (
2 * k12 - 2 * k22 + k32
)
if verbose:
print("da_mt = ", da_mt_sec, da_mt_pri)
da = da + da_mt_sec + da_mt_pri
# Tidal forces
if do_tides:
q1 = M_pri / M_sec
q2 = M_sec / M_pri
# P_orb in years. From 3rd Kepler's law, transforming separation from
# Rsolar to AU, to avoid using constants
# TODO: we aleady have this function!
P_orb = np.sqrt((a / const.aursun) ** 3 / (M_pri + M_sec))
n = 2.0 * const.pi / P_orb # mean orbital ang. vel. in rad/year
f1 = (
1
+ (31 / 2) * e ** 2
+ (255 / 8) * e ** 4
+ (185 / 16) * e ** 6
+ (25 / 64) * e ** 8
)
f2 = 1 + (15 / 2) * e ** 2 + (45 / 8) * e ** 4 + (5 / 16) * e ** 6
f3 = 1 + (15 / 4) * e ** 2 + (15 / 8) * e ** 4 + (5 / 64) * e ** 6
f4 = 1 + (3 / 2) * e ** 2 + (1 / 8) * e ** 4
f5 = 1 + 3 * e ** 2 + (3 / 8) * e ** 4
# equilibrium timecale
if ((M_env_sec != 0.0 and not np.isnan(M_env_sec))
and (DR_env_sec != 0.0 and not np.isnan(DR_env_sec)) and (
Renv_middle_sec != 0.0 and not np.isnan(Renv_middle_sec))):
# eq. (31) of Hurley et al. 2002, generalized for convective layers
# not on surface too
tau_conv_sec = 0.431 * ((M_env_sec * DR_env_sec * Renv_middle_sec
/ (3 * L_sec)) ** (1.0 / 3.0))
else:
if verbose:
print("something wrong with M_env/DR_env/Renv_middle",
M_env_sec, DR_env_sec, Renv_middle_sec)
tau_conv_sec = 1.0e99
if ((M_env_pri != 0.0 and not np.isnan(M_env_pri))
and (DR_env_pri != 0.0 and not np.isnan(DR_env_pri)) and (
Renv_middle_pri != 0.0 and not np.isnan(Renv_middle_pri))):
# eq. (31) of Hurley et al. 2002, generalized for convective layers
# not on surface too
tau_conv_pri = 0.431 * ((M_env_pri * DR_env_pri * Renv_middle_pri
/ (3 * L_pri)) ** (1.0/3.0))
else:
if verbose:
print("something wrong with M_env/DR_env/Renv_middle",
M_env_pri, DR_env_pri, Renv_middle_pri)
tau_conv_pri = 1.0e99
P_spin_sec = 2 * np.pi / Omega_sec
P_spin_pri = 2 * np.pi / Omega_pri
P_tid_sec = np.abs(1 / (1 / P_orb - 1 / P_spin_sec))
P_tid_pri = np.abs(1 / (1 / P_orb - 1 / P_spin_pri))
f_conv_sec = np.min([1, (P_tid_sec / (2 * tau_conv_sec)) ** 2])
f_conv_pri = np.min([1, (P_tid_pri / (2 * tau_conv_pri)) ** 2])
F_tid = 1. # not 50 as before
kT_conv_sec = (
(2. / 21) * (f_conv_sec / tau_conv_sec) * (M_env_sec / M_sec)
) # eq. (30) of Hurley et al. 2002
kT_conv_pri = (
(2. / 21) * (f_conv_pri / tau_conv_pri) * (M_env_pri / M_pri)
)
if kT_conv_sec is None or not np.isfinite(kT_conv_sec):
kT_conv_sec = 0.0
if kT_conv_pri is None or not np.isfinite(kT_conv_pri):
kT_conv_pri = 0.0
if verbose:
print("kT_conv_sec is", kT_conv_sec, ", set to 0.")
print("kT_conv_pri is", kT_conv_pri, ", set to 0.")
# this is the 1/timescale of all d/dt calculted below in yr^-1
if verbose:
print(
"Equilibrium tides in deep convective envelope",
M_env_sec,
DR_env_sec,
Renv_middle_sec,
R_sec,
M_sec,
M_env_pri,
DR_env_pri,
Renv_middle_pri,
R_pri,
M_pri
)
print("convective tiimescales and efficiencies",
tau_conv_sec, P_orb, P_spin_sec, P_tid_sec,
f_conv_sec,
F_tid,
tau_conv_pri, P_orb, P_spin_pri, P_tid_pri,
f_conv_pri,
F_tid,
)
# dynamical timecale
F_tid = 1
# E2 = 1.592e-9*M**(2.84) # eq. (43) of Hurley et al. 2002. Deprecated
R_conv_sec = conv_mx1_top_r_sec - conv_mx1_bot_r_sec
R_conv_pri = conv_mx1_top_r_pri - conv_mx1_bot_r_pri # convective core
# R_conv = conv_mx1_top_r # convective core
if (R_conv_sec > R_sec or R_conv_sec <= 0.0
or conv_mx1_bot_r_sec / R_sec > 0.1):
# R_conv = 0.5*R
# if verbose:
# print(
# "R_conv of the convective core is not behaving well or "
# "we are not calculating the convective core, we make it "
# "equal to half of Rstar",
# R_conv,
# R,
# conv_mx1_top_r,
# conv_mx1_bot_r,
# )
# we switch to Zahn+1975 calculation of E2
E21 = 1.592e-9 * M_sec ** (2.84)
else:
if R_sec <= 0:
E21 = 0
elif surface_h1_sec > 0.4:
E21 = 10.0 ** (-0.42) * (R_conv_sec / R_sec) ** (
7.5
) # eq. (9) of Qin et al. 2018, 616, A28
elif surface_h1_sec <= 0.4: # "HeStar":
E21 = 10.0 ** (-0.93) * (R_conv_sec / R_sec) ** (
6.7
) # eq. (9) of Qin et al. 2018, 616, A28
else: # in principle we should not go here
E21 = 1.592e-9 * M_sec ** (
2.84
) # eq. (43) of Hurley et al. 2002 from Zahn+1975 Depreciated
# kT = 1.9782e4 * np.sqrt(M * R**2 / a**5) * (1 + q)**(5. / 6) * E2
# eq. (42) of Hurley et al. 2002. Depreciated
if (R_conv_pri > R_pri or R_conv_pri <= 0.0
or conv_mx1_bot_r_pri / R_pri > 0.1):
E22 = 1.592e-9 * M_pri ** (2.84)
if verbose:
print(
"R_conv of the convective core is not behaving well or we "
"are not calculating the convective core, we switch to "
"Zahn+1975 calculation of E2",
R_conv_sec,
R_sec,
conv_mx1_top_r_sec,
conv_mx1_bot_r_sec,
E21,
R_conv_pri,
R_pri,
conv_mx1_top_r_pri,
conv_mx1_bot_r_pri,
E22
)
else:
if R_pri <= 0:
E22 = 0
elif surface_h1_pri > 0.4:
E22 = 10.0 ** (-0.42) * (R_conv_pri / R_pri) ** (
7.5
) # eq. (9) of Qin et al. 2018, 616, A28
elif surface_h1_pri <= 0.4: # "HeStar":
E22 = 10.0 ** (-0.93) * (R_conv_pri / R_pri) ** (
6.7
) # eq. (9) of Qin et al. 2018, 616, A28
else: # in principle we should not go here
E22 = 1.592e-9 * M_pri ** (
2.84
)
kT_rad_sec = (
np.sqrt(const.standard_cgrav * (M_sec * const.msol)
* (R_sec * const.rsol)**2 / (a * const.rsol)**5)
* (1 + q1) ** (5.0 / 6)
* E21
* const.secyer)
kT_rad_pri = (
np.sqrt(const.standard_cgrav * (M_pri * const.msol)
* (R_pri * const.rsol)**2 / (a * const.rsol)**5)
* (1 + q2) ** (5.0 / 6)
* E22
* const.secyer)
# this is the 1/timescale of all d/dt calculted below in yr^-1
if verbose:
print(
"Dynamical tides in radiative envelope",
conv_mx1_top_r_sec,
conv_mx1_bot_r_sec,
R_conv_sec,
E21,
conv_mx1_top_r_pri,
conv_mx1_bot_r_pri,
R_conv_pri,
E22,
F_tid
)
kT_sec = max(kT_conv_sec, kT_rad_sec)
kT_pri = max(kT_conv_pri, kT_rad_pri)
if verbose:
print("kT_conv/rad of tides is ", kT_conv_sec, kT_rad_sec,
kT_conv_pri, kT_rad_pri, "in 1/yr, and we picked the ",
kT_sec, kT_pri)
da_tides_sec = (
-6
* F_tid
* kT_sec
* q1
* (1 + q1)
* (R_sec / a) ** 8
* (a / (1 - e ** 2) ** (15 / 2))
* (f1 - (1 - e ** 2) ** (3 / 2) * f2 * Omega_sec / n)
) # eq. (9) of Hut 1981, 99, 126
da_tides_pri = (
-6
* F_tid
* kT_pri
* q2
* (1 + q2)
* (R_pri / a) ** 8
* (a / (1 - e ** 2) ** (15 / 2))
* (f1 - (1 - e ** 2) ** (3 / 2) * f2 * Omega_pri / n)
)
de_tides_sec = (
-27
* F_tid
* kT_sec
* q1
* (1 + q1)
* (R_sec / a) ** 8
* (e / (1 - e ** 2) ** (13 / 2))
* (f3 - (11 / 18) * (1 - e ** 2) ** (3 / 2) * f4 * Omega_sec/n)
) # eq. (10) of Hut 1981, 99, 126
de_tides_pri = (
-27
* F_tid
* kT_pri
* q2
* (1 + q2)
* (R_pri / a) ** 8
* (e / (1 - e ** 2) ** (13 / 2))
* (f3 - (11 / 18) * (1 - e ** 2) ** (3 / 2) * f4 * Omega_pri/n)
)
dOmega_tides_sec = (
(3 * F_tid * kT_sec * q1 ** 2 * M_sec * R_sec ** 2 / I_sec)
* (R_sec / a) ** 6
* n
/ (1 - e ** 2) ** 6
* (f2 - (1 - e ** 2) ** (3 / 2) * f5 * Omega_sec / n)
) # eq. (11) of Hut 1981, 99, 126
dOmega_tides_pri = (
(3 * F_tid * kT_pri * q2 ** 2 * M_pri * R_pri ** 2 / I_pri)
* (R_pri / a) ** 6
* n
/ (1 - e ** 2) ** 6
* (f2 - (1 - e ** 2) ** (3 / 2) * f5 * Omega_pri / n)
)
if verbose:
print("da,de,dOmega_tides = ",
da_tides_sec, de_tides_sec, dOmega_tides_sec,
da_tides_pri, de_tides_pri, dOmega_tides_pri)
da = da + da_tides_sec + da_tides_pri
de = de + de_tides_sec + de_tides_pri
dOmega_sec = dOmega_sec + dOmega_tides_sec
dOmega_pri = dOmega_pri + dOmega_tides_pri
# Gravitional radiation
if do_gravitational_radiation:
v = (M_pri * M_sec / (M_pri + M_sec) ** 2)
da_gr = (
(-2 * const.clight / 15) * (v / ((1 - e ** 2) ** (9 / 2)))
* (const.standard_cgrav * (M_pri + M_sec) * const.msol
/ (a * const.rsol * const.clight ** 2)) ** 3
* ((96 + 292 * e ** 2 + 37 * e ** 4) * (1 - e ** 2)
- (1 / 28 * const.standard_cgrav * (M_pri + M_sec) * const.msol
/ (a * const.rsol * const.clight ** 2))
* ((14008 + 4707 * v)+(80124 + 21560 * v) * e ** 2
+ (17325 + 10458) * e ** 4 - 0.5 * (5501 - 1036 * v) * e ** 6))
) * const.secyer / const.rsol
# eq. (35) of Junker et al. 1992, 254, 146
de_gr = (
(-1 / 15) * ((v * const.clight ** 3) / (
const.standard_cgrav * (M_pri + M_sec) * const.msol))
* (const.standard_cgrav * (M_pri + M_sec) * const.msol / (
a * const.rsol * const.clight ** 2)) ** 4
* (e / (1 - e ** 2) ** (7 / 2))
* ((304 + 121 * e ** 2) * (1 - e ** 2)
- (1 / 56 * const.standard_cgrav * (M_pri + M_sec) * const.msol
/ (a * const.rsol * const.clight ** 2))
* (8 * (16705 + 4676 * v) + 12 * (9082 + 2807 * v) * e ** 2
- (25211 - 3388 * v) * e ** 4))
) * const.secyer
# eq. (36) of Junker et al. 1992, 254, 146
if verbose:
print("da,de_gr = ", da_gr, de_gr)
da = da + da_gr
de = de + de_gr
# Magnetic braking
if do_magnetic_braking:
# domega_mb / dt = torque_mb / I,
# where torque_mb is in eq.36 of Rapport+1983, with γ = 4
dOmega_mb_sec = (
-6.82e34
* M_sec
* R_sec ** 4
* (Omega_sec * 365.25 / (2 * np.pi)) ** 3
* 1.48763e-49
/ I_sec
* np.clip((1.5 - M_sec) / (1.5 - 1.3), 0, 1)
)
dOmega_mb_pri = (
-6.82e34
* M_pri
* R_pri ** 4
* (Omega_pri * 365.25 / (2 * np.pi)) ** 3
* 1.48763e-49
/ I_pri
* np.clip((1.5 - M_pri) / (1.5 - 1.3), 0, 1)
)
# -3.8*10e-30
# * M
# * R ** 4
# * Omega ** 3
# / I
# * (const.rsol)**2 * (const.secyer)**4
# (rsol^2 because I ~ MR^2 in the denominator and secyer^4
# from omega^3 and to make domega_mb in rad/yr)
if verbose:
print("dOmega_mb = ", dOmega_mb_sec, dOmega_mb_pri)
dOmega_sec = dOmega_sec + dOmega_mb_sec
dOmega_pri = dOmega_pri + dOmega_mb_pri
if do_stellar_evolution_and_spin_from_winds:
# Due to the secondary's own evolution, we have:
# domega_spin/dt = d(Jspin/I)/dt = dJspin/dt * 1/I + Jspin*d(1/I)/dt.
# These are the two terms calculated below.
dOmega_spin_wind_sec = (
2.0 / 3.0 * R_sec ** 2 * Omega_sec * Mdot_sec / I_sec
)
# jshell*Mdot/I : specific angular momentum of a thin sperical shell
# * mdot / moment of inertia
# Omega is in rad/yr here, and R, M in solar (Mdot solar/yr).
dOmega_deformation_sec = np.min(
[-Omega_sec * Idot_sec / I_sec, 100]
)
# This is term of Jspin*d(1/I)/dt term of the domega_spin/dt. #
# We limit its increase due to contraction to 100 [(rad/yr)/yr]
# increase, otherwise the integrator will fail
# (usually when we have WD formation).
dOmega_spin_wind_pri = (
2.0 / 3.0 * R_pri ** 2 * Omega_pri * Mdot_pri / I_pri
)
# jshell*Mdot/I : specific angular momentum of a thin sperical shell
# * mdot / moment of inertia
# Omega is in rad/yr here, and R, M in solar (Mdot solar/yr).
dOmega_deformation_pri = np.min(
[-Omega_pri * Idot_pri / I_pri, 100]
)
if verbose:
print(
"dOmega_spin_wind , dOmega_deformation = ",
dOmega_spin_wind_sec,
dOmega_deformation_sec,
dOmega_spin_wind_pri,
dOmega_deformation_pri,
)
dOmega_sec = dOmega_sec + dOmega_spin_wind_sec
dOmega_pri = dOmega_pri + dOmega_spin_wind_pri
dOmega_sec = dOmega_sec + dOmega_deformation_sec
dOmega_pri = dOmega_pri + dOmega_deformation_pri
if verbose:
print("a[Ro],e,Omega[rad/yr] have been =", a, e, Omega_sec, Omega_pri)
print("da,de,dOmega (all) in 1yr is = ",
da, de, dOmega_sec, dOmega_pri)
return [
da,
de,
dOmega_sec,
dOmega_pri,
]