diff --git a/chemcomp/accretion_models/__init__.py b/chemcomp/accretion_models/__init__.py
index e594695..ed1cf95 100755
--- a/chemcomp/accretion_models/__init__.py
+++ b/chemcomp/accretion_models/__init__.py
@@ -1,3 +1,3 @@
from ._accretion_class import Accretion
from .gas import GasAccretion
-from .pebbles import PebbleAccretion
\ No newline at end of file
+from .pebbles import PebbleAccretion
diff --git a/chemcomp/accretion_models/_accretion_class.py b/chemcomp/accretion_models/_accretion_class.py
index af17ae6..b6cd5b1 100755
--- a/chemcomp/accretion_models/_accretion_class.py
+++ b/chemcomp/accretion_models/_accretion_class.py
@@ -13,11 +13,11 @@ def __init__(self, config=None):
self.regime = ""
def calc_m_dot(self):
- """ function that calculates and updates the accretion model """
+ """function that calculates and updates the accretion model"""
pass
def accretion_domain(self):
- """ function that calculates the regimes of the accretion model """
+ """function that calculates the regimes of the accretion model"""
pass
def init_accretion(self, planet):
@@ -32,11 +32,11 @@ def init_accretion(self, planet):
self._init_params()
def _init_params(self):
- """ function that inits parameters needed """
+ """function that inits parameters needed"""
pass
def update_z(self):
- """update the total mass of heavy elements that have been accreted. """
+ """update the total mass of heavy elements that have been accreted."""
self.M_z[0] += np.sum(self.m_dot_c_chem[0, :-7] * self.planet.h)
self.M_z[1] += np.sum(self.m_dot_a_chem[0, :-7] * self.planet.h)
diff --git a/chemcomp/accretion_models/gas.py b/chemcomp/accretion_models/gas.py
index 27e1fbc..9e9039f 100755
--- a/chemcomp/accretion_models/gas.py
+++ b/chemcomp/accretion_models/gas.py
@@ -10,14 +10,14 @@
class GasAccretion(Accretion):
- """ Model for GasAccretion.
+ """Model for GasAccretion.
Minimum of Ikoma2000, Machida2010 and disk gas accretion rate (including horseshoe region).
"""
def __init__(self, config):
super().__init__()
self.name = "GasAccretion"
- self.kappa_env = config.get("kappa_env", 0.05 * u.cm ** 2 / u.g)
+ self.kappa_env = config.get("kappa_env", 0.05 * u.cm**2 / u.g)
self.kappa_env = self.kappa_env.cgs.value
self.f = config.get("f", 0.2)
self.f_machida = config.get("f_machida", 1)
@@ -25,10 +25,7 @@ def __init__(self, config):
def _init_params(self):
self.const_terms = (
- 1e3
- * (30*earth_mass_cgs_value)**2.5
- * self.kappa_env/0.05
- * year
+ 1e3 * (30 * earth_mass_cgs_value) ** 2.5 * self.kappa_env / 0.05 * year
)
def accretion_domain(self):
@@ -41,23 +38,23 @@ def remove_mass_from_flux(self):
return
def get_min(self):
- """" Get the minimal accretion rate """
+ """ " Get the minimal accretion rate"""
low = (
- 0.83
- * self.f_machida
- * self.planet.omega_k
- * self.planet.H_g ** 2
- * self.planet.sigma_g
- * (self.planet.r_h / self.planet.H_g) ** (9 / 2)
+ 0.83
+ * self.f_machida
+ * self.planet.omega_k
+ * self.planet.H_g**2
+ * self.planet.sigma_g
+ * (self.planet.r_h / self.planet.H_g) ** (9 / 2)
)
low = np.maximum(low, 1e-60)
high = (
- 0.14
- * self.f_machida
- * self.planet.omega_k
- * self.planet.H_g ** 2
- * self.planet.sigma_g
+ 0.14
+ * self.f_machida
+ * self.planet.omega_k
+ * self.planet.H_g**2
+ * self.planet.sigma_g
)
high = np.maximum(high, 1e-60)
@@ -69,7 +66,7 @@ def get_min(self):
mdot_HS = np.maximum(mdot_HS, 1e-60)
disk_max = disk_max + mdot_HS
- ikoma = self.planet.M/(self.planet.M_c**(-2.5)*self.const_terms)
+ ikoma = self.planet.M / (self.planet.M_c ** (-2.5) * self.const_terms)
ikoma = np.maximum(ikoma, 1e-60)
valid_accrates = [low, high, disk_max, ikoma]
@@ -80,7 +77,7 @@ def get_min(self):
return str(minimum), mdot
def calc_m_dot(self):
- """ main routine that calculates the accretion rate."""
+ """main routine that calculates the accretion rate."""
self.update_parameters()
self.accretion_domain()
self.m_dot_c = 0
@@ -96,6 +93,5 @@ def calc_m_dot(self):
self.m_dot = self.m_dot_a + self.m_dot_c
self.m_dot_a_chem = self.chem_mask_array * self.m_dot_a
-
assert np.all(self.m_dot_a_chem >= 0)
- return
\ No newline at end of file
+ return
diff --git a/chemcomp/accretion_models/pebbles.py b/chemcomp/accretion_models/pebbles.py
index 29b12f1..8646859 100755
--- a/chemcomp/accretion_models/pebbles.py
+++ b/chemcomp/accretion_models/pebbles.py
@@ -8,7 +8,7 @@
class PebbleAccretion(Accretion):
- """ Model for PebbleAccretion. """
+ """Model for PebbleAccretion."""
def __init__(self, config):
super().__init__()
@@ -18,7 +18,7 @@ def __init__(self, config):
self.u_frag = config.get("u_frag", None)
self.t_f = config.get("STOKES", None)
- self.factor_peb_accretion = config.get('factor_peb_accretion', 1.0)
+ self.factor_peb_accretion = config.get("factor_peb_accretion", 1.0)
if self.rho_solid is None:
self._calculate_rho_solid = True
@@ -28,7 +28,9 @@ def __init__(self, config):
if self.u_frag is None:
self._calculate_u_frag = True
- print("Warning: No fixed fragmentation velocity has been set! The code will calculate the fragmentation velocity according to the Drazkowska 2017 approach: 1m/s if silicate, 10m/s if icy.")
+ print(
+ "Warning: No fixed fragmentation velocity has been set! The code will calculate the fragmentation velocity according to the Drazkowska 2017 approach: 1m/s if silicate, 10m/s if icy."
+ )
else:
self.u_frag = self.u_frag.cgs.value
self._calculate_u_frag = False
@@ -37,34 +39,42 @@ def __init__(self, config):
self.twoD = False
def _init_params(self):
- """ Pass some pebble parameters to the disk. """
- self.planet.disk.init_pebble_parameters(self.epsilon_p, self.rho_solid, self.t_f, self.u_frag,
- self._calculate_rho_solid, self._calculate_u_frag)
+ """Pass some pebble parameters to the disk."""
+ self.planet.disk.init_pebble_parameters(
+ self.epsilon_p,
+ self.rho_solid,
+ self.t_f,
+ self.u_frag,
+ self._calculate_rho_solid,
+ self._calculate_u_frag,
+ )
def calc_stokes(self):
- """ Compute the stokes number by interpolating the stokes number of the disk to the planets location. """
- self.t_f = np.interp(self.planet.a_p, self.planet.disk.r, self.planet.disk.stokes_number_pebbles)
+ """Compute the stokes number by interpolating the stokes number of the disk to the planets location."""
+ self.t_f = np.interp(
+ self.planet.a_p, self.planet.disk.r, self.planet.disk.stokes_number_pebbles
+ )
self.tau_f = self.t_f / self.planet.omega_k
def calc_densities(self):
- """ calculate the pebble scalehight and the volume density of pebbles. """
+ """calculate the pebble scalehight and the volume density of pebbles."""
self.H_p = np.sqrt(self.planet.disk.alpha_height / self.t_f) * self.planet.H_g
- self.rho_peb_components = self.planet.sigma_peb_components / (np.sqrt(2 * np.pi) * self.H_p)
+ self.rho_peb_components = self.planet.sigma_peb_components / (
+ np.sqrt(2 * np.pi) * self.H_p
+ )
return
def accretion_domain(self): # checked
- """ decide on the regimes of accretion: bondi vs hill """
+ """decide on the regimes of accretion: bondi vs hill"""
if self.regime == "Bondi":
transition_mass = (
- np.sqrt(1 / 3)
- * self.planet.del_v ** 3
- / (G * self.planet.omega_k)
+ np.sqrt(1 / 3) * self.planet.del_v**3 / (G * self.planet.omega_k)
)
if self.planet.M > transition_mass:
self.regime = "Hill"
def twoD_transition(self): # checked
- """ Transition between 2D and 3D pebble accretion. """
+ """Transition between 2D and 3D pebble accretion."""
transition = np.pi * self.R_acc / (2 * np.sqrt(2 * np.pi))
if transition > self.H_p:
self.twoD = True
@@ -72,7 +82,7 @@ def twoD_transition(self): # checked
self.twoD = False
def calc_R_acc(self): # checked
- """Calculate the capture radius of pebble accretion. """
+ """Calculate the capture radius of pebble accretion."""
if self.regime == "Bondi":
t_b = self.planet.r_b / self.planet.del_v # checked
self.R_acc = np.sqrt(4 * self.tau_f / t_b) * self.planet.r_b # checked
@@ -83,17 +93,22 @@ def calc_R_acc(self): # checked
if self.regime in ["Hill", "Bondi"]:
t_p = (
- G
- * self.planet.M
- / ((np.abs(self.planet.del_v) + self.planet.omega_k * self.planet.r_h) ** 3)
+ G
+ * self.planet.M
+ / (
+ (np.abs(self.planet.del_v) + self.planet.omega_k * self.planet.r_h)
+ ** 3
+ )
+ ) # checked
+ self.R_acc = self.R_acc * np.exp(
+ -0.4 * (self.tau_f / t_p) ** 0.65
) # checked
- self.R_acc = self.R_acc * np.exp(-0.4 * (self.tau_f / t_p) ** 0.65) # checked
self.del_v = self.planet.del_v + self.planet.omega_k * self.R_acc # checked
return
def remove_mass_from_flux(self):
- """ Ask the disk to remove the accreted pebbles."""
+ """Ask the disk to remove the accreted pebbles."""
self.planet.disk.remove_from_peb(self.m_dot_chem, self.planet.idx)
return
@@ -107,13 +122,19 @@ def calc_m_dot(self):
self.twoD_transition()
if self.twoD:
- self.m_dot_chem = 2 * self.R_acc * self.planet.sigma_peb_components * self.del_v # checked
+ self.m_dot_chem = (
+ 2 * self.R_acc * self.planet.sigma_peb_components * self.del_v
+ ) # checked
else:
- self.m_dot_chem = np.pi * self.R_acc ** 2 * self.rho_peb_components * self.del_v # checked
+ self.m_dot_chem = (
+ np.pi * self.R_acc**2 * self.rho_peb_components * self.del_v
+ ) # checked
self.m_dot_chem = self.factor_peb_accretion * self.m_dot_chem
- self.m_dot_c_chem = (1 - self.planet.solid_frac) * self.m_dot_chem # checked
+ self.m_dot_c_chem = (
+ 1 - self.planet.solid_frac
+ ) * self.m_dot_chem # checked
self.m_dot_a_chem = self.planet.solid_frac * self.m_dot_chem # checked
self.m_dot = np.sum(self.m_dot_chem[0]) # sum over all elements
@@ -121,7 +142,6 @@ def calc_m_dot(self):
self.m_dot_c = (1 - self.planet.solid_frac) * self.m_dot # checked
self.m_dot_a = self.planet.solid_frac * self.m_dot # checked
-
else:
self.m_dot = 0
self.m_dot_c = 1 * self.m_dot
diff --git a/chemcomp/disks/_0pacity/opacity_models.py b/chemcomp/disks/_0pacity/opacity_models.py
index 337026b..d4adc9b 100755
--- a/chemcomp/disks/_0pacity/opacity_models.py
+++ b/chemcomp/disks/_0pacity/opacity_models.py
@@ -2,10 +2,15 @@
from chemcomp.helpers.units.berts_units import *
-opacity_array = np.genfromtxt(os.path.join(os.path.dirname(os.path.realpath(__file__)), "meanopac5.dat"),
- usecols=(0, 1), skip_header=1).T
+opacity_array = np.genfromtxt(
+ os.path.join(os.path.dirname(os.path.realpath(__file__)), "meanopac5.dat"),
+ usecols=(0, 1),
+ skip_header=1,
+).T
opacity_array_derivative = np.gradient(opacity_array[1], opacity_array[0])
-opacity_array_derivative_derivative = np.gradient(opacity_array_derivative, opacity_array[0])
+opacity_array_derivative_derivative = np.gradient(
+ opacity_array_derivative, opacity_array[0]
+)
def belllin(rho, temp, Z=0.01, onlygas=False):
@@ -19,26 +24,30 @@ def belllin(rho, temp, Z=0.01, onlygas=False):
onlygas If set, then do not include the dust part of the _0pacity
"""
if onlygas:
- ki = np.array([1e-8,1e-36,1.5e20,0.348])
- a = np.array([0.6667,0.3333,1.,0.])
- b = np.array([3.,10.,-2.5,0.])
+ ki = np.array([1e-8, 1e-36, 1.5e20, 0.348])
+ a = np.array([0.6667, 0.3333, 1.0, 0.0])
+ b = np.array([3.0, 10.0, -2.5, 0.0])
else:
- ki = np.array([2e-4,2e16,0.1,2e81,1e-8,1e-36,1.5e20,0.348])
- a = np.array([0.,0.,0.,1.,0.6667,0.3333,1.,0.])
- b = np.array([2.,-7.,0.5,-24.,3.,10.,-2.5,0.])
- nn = len(ki)
+ ki = np.array([2e-4, 2e16, 0.1, 2e81, 1e-8, 1e-36, 1.5e20, 0.348])
+ a = np.array([0.0, 0.0, 0.0, 1.0, 0.6667, 0.3333, 1.0, 0.0])
+ b = np.array([2.0, -7.0, 0.5, -24.0, 3.0, 10.0, -2.5, 0.0])
+ nn = len(ki)
- n = len(rho)
- dm = np.zeros((nn,n))
+ n = len(rho)
+ dm = np.zeros((nn, n))
for ii in range(nn):
- dm[ii,:] = ki[ii]*rho[:]**a[ii]
- tc = np.zeros((nn+1,n))
- for ii in range(nn-1):
- tc[ii+1,:] = (dm[ii,:]/dm[ii+1,:])**(1./(b[ii+1]-b[ii]))
- tc[nn,:] = 1e99
+ dm[ii, :] = ki[ii] * rho[:] ** a[ii]
+ tc = np.zeros((nn + 1, n))
+ for ii in range(nn - 1):
+ tc[ii + 1, :] = (dm[ii, :] / dm[ii + 1, :]) ** (1.0 / (b[ii + 1] - b[ii]))
+ tc[nn, :] = 1e99
kappa = np.zeros_like(rho)
for ii in range(nn):
- kappa = np.where(np.logical_and(temp>tc[ii,:],temp<tc[ii+1,:]), dm[ii]*temp**b[ii],kappa)
+ kappa = np.where(
+ np.logical_and(temp > tc[ii, :], temp < tc[ii + 1, :]),
+ dm[ii] * temp ** b[ii],
+ kappa,
+ )
return kappa * Z / 0.01
@@ -50,7 +59,7 @@ def oplin(_Rho, _Temp, Z=0.01):
"""
ts4 = 1.0e-4 * _Temp
- rho13 = _Rho ** 0.333333333
+ rho13 = _Rho**0.333333333
rho23 = rho13 * rho13
ts42 = ts4 * ts4
ts44 = ts42 * ts42
@@ -59,9 +68,9 @@ def oplin(_Rho, _Temp, Z=0.01):
# init with no scattering
opacity = bk7 * _Rho / (ts42 * np.sqrt(ts4))
- m_0 = _Temp <= t456 * _Rho ** power2
- m_0_1 = np.logical_and(m_0, _Temp > t234 * _Rho ** power1)
- m_1 = np.logical_or(_Temp < t678 * _Rho ** power3, _Rho < 1.0e-10)
+ m_0 = _Temp <= t456 * _Rho**power2
+ m_0_1 = np.logical_and(m_0, _Temp > t234 * _Rho**power1)
+ m_1 = np.logical_or(_Temp < t678 * _Rho**power3, _Rho < 1.0e-10)
# disjoint _0pacity laws for 5, 6, and 7.
o5 = bk5 * rho23[m_1] * ts42[m_1] * ts4[m_1]
@@ -82,8 +91,8 @@ def oplin(_Rho, _Temp, Z=0.01):
o4 = bk4 * rho23[m_0_1] / (ts48[m_0_1] * ts4[m_0_1])
o5 = bk5 * rho23[m_0_1] * ts42[m_0_1] * ts4[m_0_1]
# parameters used for smoothing
- o4an = o4 ** 4
- o3an = o3 ** 4
+ o4an = o4**4
+ o3an = o3**4
# smoothed and continuous _0pacity law for regions 3, 4, and 5.
opacity[m_0_1] = (
@@ -107,7 +116,7 @@ def oplin(_Rho, _Temp, Z=0.01):
opacity[m_0] = (
(o1an * o2an / (o1an + o2an)) ** 2 + (o3 / (1 + 1.0e22 / t10)) ** 4
- ) ** 0.25
+ ) ** 0.25
return opacity * Z / 0.01
diff --git a/chemcomp/disks/_chemistry.py b/chemcomp/disks/_chemistry.py
index e9304a7..33162d7 100755
--- a/chemcomp/disks/_chemistry.py
+++ b/chemcomp/disks/_chemistry.py
@@ -74,87 +74,105 @@ def split_molecules(abundances_inp):
"""
(
- rest_mol,
- TotMassCO,
- TotMassN2,
- TotMassCH4,
- TotMassCO2,
- TotMassNH3,
- TotMassH2S,
- TotMassH2O,
- TotMassFe3O4,
- TotMassC,
- TotMassFeS,
- TotMassNaAlSi3O8,
- TotMassKAlSi3O8,
- TotmassMg2SiO4,
- TotMassFe2O3,
- TotMassVO,
- TotMassMgSiO3,
- TotMassAl2O3,
- TotMassTiO,
- ) = abundances_inp
+ rest_mol,
+ TotMassCO,
+ TotMassN2,
+ TotMassCH4,
+ TotMassCO2,
+ TotMassNH3,
+ TotMassH2S,
+ TotMassH2O,
+ TotMassFe3O4,
+ TotMassC,
+ TotMassFeS,
+ TotMassNaAlSi3O8,
+ TotMassKAlSi3O8,
+ TotmassMg2SiO4,
+ TotMassFe2O3,
+ TotMassVO,
+ TotMassMgSiO3,
+ TotMassAl2O3,
+ TotMassTiO,
+ ) = abundances_inp
TotC = (
- TotMassCO * MassC / MassCO
- + TotMassCH4 * MassC / MassCH4
- + TotMassCO2 * MassC / MassCO2
- + TotMassC * MassC / MassC
+ TotMassCO * MassC / MassCO
+ + TotMassCH4 * MassC / MassCH4
+ + TotMassCO2 * MassC / MassCO2
+ + TotMassC * MassC / MassC
)
TotO = (
- TotMassCO * MassO / MassCO
- + TotMassCO2 * MassO * 2.0 / MassCO2
- + TotMassH2O * MassO / MassH2O
- + TotMassFe3O4 * 4.0 * MassO / MassFe3O4
- + TotmassMg2SiO4 * 4.0 * MassO / MassMg2SiO4
- + TotMassFe2O3 * 3 * MassO / MassFe2O3
- + TotMassMgSiO3 * MassO * 3.0 / MassMgSiO3
- + TotMassTiO * MassO / MassTiO
- + TotMassAl2O3 * 3 * MassO / MassAl2O3
- + TotMassKAlSi3O8 * 8 * MassO / MassKAlSi3O8
- + TotMassNaAlSi3O8 * 8 * MassO / MassNaAlSi3O8
- + TotMassVO * MassO / MassVO
+ TotMassCO * MassO / MassCO
+ + TotMassCO2 * MassO * 2.0 / MassCO2
+ + TotMassH2O * MassO / MassH2O
+ + TotMassFe3O4 * 4.0 * MassO / MassFe3O4
+ + TotmassMg2SiO4 * 4.0 * MassO / MassMg2SiO4
+ + TotMassFe2O3 * 3 * MassO / MassFe2O3
+ + TotMassMgSiO3 * MassO * 3.0 / MassMgSiO3
+ + TotMassTiO * MassO / MassTiO
+ + TotMassAl2O3 * 3 * MassO / MassAl2O3
+ + TotMassKAlSi3O8 * 8 * MassO / MassKAlSi3O8
+ + TotMassNaAlSi3O8 * 8 * MassO / MassNaAlSi3O8
+ + TotMassVO * MassO / MassVO
)
TotFe = (
- TotMassFe3O4 * 3.0 * MassFe / MassFe3O4
- + TotMassFe2O3 * 2.0 * MassFe / MassFe2O3
- + TotMassFeS * MassFe / MassFeS
+ TotMassFe3O4 * 3.0 * MassFe / MassFe3O4
+ + TotMassFe2O3 * 2.0 * MassFe / MassFe2O3
+ + TotMassFeS * MassFe / MassFeS
)
TotS = TotMassFeS * MassS / MassFeS + TotMassH2S * MassS / MassH2S
TotMg = (
- TotmassMg2SiO4 * MassMg * 2.0 / MassMg2SiO4
- + TotMassMgSiO3 * MassMg / MassMgSiO3
+ TotmassMg2SiO4 * MassMg * 2.0 / MassMg2SiO4
+ + TotMassMgSiO3 * MassMg / MassMgSiO3
)
TotSi = (
- TotmassMg2SiO4 * MassSi / MassMg2SiO4
- + TotMassMgSiO3 * MassSi / MassMgSiO3
- + TotMassKAlSi3O8 * 3 * MassSi / MassKAlSi3O8
- + TotMassNaAlSi3O8 * 3 * MassSi / MassNaAlSi3O8
-
+ TotmassMg2SiO4 * MassSi / MassMg2SiO4
+ + TotMassMgSiO3 * MassSi / MassMgSiO3
+ + TotMassKAlSi3O8 * 3 * MassSi / MassKAlSi3O8
+ + TotMassNaAlSi3O8 * 3 * MassSi / MassNaAlSi3O8
)
TotK = TotMassKAlSi3O8 * MassK / MassKAlSi3O8
TotNa = TotMassNaAlSi3O8 * MassNa / MassNaAlSi3O8
TotAl = (
- TotMassKAlSi3O8 * MassAl / MassKAlSi3O8
- + TotMassNaAlSi3O8 * MassAl / MassNaAlSi3O8
- + TotMassAl2O3 * 2 * MassAl / MassAl2O3
+ TotMassKAlSi3O8 * MassAl / MassKAlSi3O8
+ + TotMassNaAlSi3O8 * MassAl / MassNaAlSi3O8
+ + TotMassAl2O3 * 2 * MassAl / MassAl2O3
)
TotTi = TotMassTiO * MassTi / MassTiO
TotV = TotMassVO * MassV / MassVO
TotN = TotMassNH3 * MassN / MassNH3 + TotMassN2 * 2 * MassN / MassN2
TotH_mol = (
- TotMassCH4 * 4.0 * MassH / MassCH4
- + TotMassH2O * 2.0 * MassH / MassH2O
- + TotMassNH3 * 3 * MassH / MassNH3
- + TotMassH2S * 2 * MassH / MassH2S
+ TotMassCH4 * 4.0 * MassH / MassCH4
+ + TotMassH2O * 2.0 * MassH / MassH2O
+ + TotMassNH3 * 3 * MassH / MassNH3
+ + TotMassH2S * 2 * MassH / MassH2S
)
# Note how the last entry of the array is TotH_mol
- return np.array([TotC, TotO, TotFe, TotS, TotMg, TotSi, TotNa, TotK, TotN, TotAl, TotTi, TotV, TotH_mol])
+ return np.array(
+ [
+ TotC,
+ TotO,
+ TotFe,
+ TotS,
+ TotMg,
+ TotSi,
+ TotNa,
+ TotK,
+ TotN,
+ TotAl,
+ TotTi,
+ TotV,
+ TotH_mol,
+ ]
+ )
+
-def calc_composition(abundances_inp: np.array, solid: bool, He_Mbg: float) -> Tuple[np.array, np.array]:
+def calc_composition(
+ abundances_inp: np.array, solid: bool, He_Mbg: float
+) -> Tuple[np.array, np.array]:
"""
Parameters
@@ -169,21 +187,23 @@ def calc_composition(abundances_inp: np.array, solid: bool, He_Mbg: float) -> Tu
Nd numpy array: the total mass of Hydrogen in molecules per grid cell
"""
-
+
molecule_elements = split_molecules(abundances_inp)
rest_mol = abundances_inp[0]
rest = np.sum(molecule_elements[:-1], axis=0)
TotH_mol = molecule_elements[-1]
- if solid: # These are the solid components, getting the H, He abundance is trivial
+ if solid: # These are the solid components, getting the H, He abundance is trivial
TotH = TotH_mol
TotHe = np.zeros_like(rest)
- else: # The gas components, where information about the Helium content is required and therefore the supplied He_Mbg argument is used
+ else: # The gas components, where information about the Helium content is required and therefore the supplied He_Mbg argument is used
TotHe = He_Mbg * rest_mol
- TotH = (1 - rest - TotHe)
-
+ TotH = 1 - rest - TotHe
+
# Test about the He+H split
- assert np.isclose(TotH + TotHe - TotH_mol, rest_mol).all(), "Background gas after split does not equal the total background gas from the molecular abundances"
+ assert np.isclose(
+ TotH + TotHe - TotH_mol, rest_mol
+ ).all(), "Background gas after split does not equal the total background gas from the molecular abundances"
elements = np.array(
[
@@ -198,15 +218,18 @@ def calc_composition(abundances_inp: np.array, solid: bool, He_Mbg: float) -> Tu
],
)
- molecules = abundances_inp.copy() # Avoid unwanted side-effects when returning the supplied molecular abundances
+ molecules = (
+ abundances_inp.copy()
+ ) # Avoid unwanted side-effects when returning the supplied molecular abundances
return np.transpose(np.stack([elements, molecules]), (2, 0, 1)), TotH_mol
class Chemistry:
- """ Chemistry class. Please be careful when changing or adding molecular species,
+ """Chemistry class. Please be careful when changing or adding molecular species,
some occurancies of the chemistry rely (hardcoded) on the order of the species!
"""
+
def __init__(self, config):
# Adundances in terms of H (Asplund et al. 2009) + variation from Chiara plot
if not config.get("use_FeH", True):
@@ -263,7 +286,9 @@ def __init__(self, config):
self.abuFe2O3 = 0.25 * (self.FeHabu - 0.9 * self.SHabu)
self.abuFe3O4 = (1.0 / 6.0) * (self.FeHabu - 0.9 * self.SHabu)
self.abuFeS = 0.9 * self.SHabu
- self.abuMgSiO3 = self.MgHabu - 2 * (self.MgHabu - (self.SiHabu - 3 * self.KHabu - 3 * self.NaHabu))
+ self.abuMgSiO3 = self.MgHabu - 2 * (
+ self.MgHabu - (self.SiHabu - 3 * self.KHabu - 3 * self.NaHabu)
+ )
self.abuMg2SiO4 = self.MgHabu - (self.SiHabu - 3 * self.KHabu - 3 * self.NaHabu)
self.abuCO = float(config.get("CO_frac", 0.45)) * self.CHabu
@@ -284,21 +309,22 @@ def __init__(self, config):
self.abuAl2O3 = 0.5 * (self.AlHabu - (self.KHabu + self.NaHabu))
self.abuH2O = self.OHabu - (
- 3 * self.abuMgSiO3
- + 4 * self.abuMg2SiO4
- + self.abuCO
- + 2 * self.abuCO2
- + 3 * self.abuFe2O3
- + 4 * self.abuFe3O4
- + self.abuTiO
- + 3 * self.abuAl2O3
- + 8 * self.abuNaAlSi3O8
- + 8 * self.abuKAlSi3O8
- + self.abuVO
+ 3 * self.abuMgSiO3
+ + 4 * self.abuMg2SiO4
+ + self.abuCO
+ + 2 * self.abuCO2
+ + 3 * self.abuFe2O3
+ + 4 * self.abuFe3O4
+ + self.abuTiO
+ + 3 * self.abuAl2O3
+ + 8 * self.abuNaAlSi3O8
+ + 8 * self.abuKAlSi3O8
+ + self.abuVO
)
- assert np.isclose(self.abuC + self.abuCO + self.abuCO2 + self.abuCH4, self.CHabu), \
- "Broken chem model, we loose or gain C"
+ assert np.isclose(
+ self.abuC + self.abuCO + self.abuCO2 + self.abuCH4, self.CHabu
+ ), "Broken chem model, we loose or gain C"
self.abu_array = np.array(
[
@@ -383,22 +409,28 @@ def get_position_of_ice(self, T: np.array) -> np.array:
idx = np.minimum(idx, T.size - 2)
return idx
-
- def calc_He_Mbg(self, abundances_gas: np.array, abundances_solid: np.array, nHe_nH: np.array, T: np.array) -> float:
+
+ def calc_He_Mbg(
+ self,
+ abundances_gas: np.array,
+ abundances_solid: np.array,
+ nHe_nH: np.array,
+ T: np.array,
+ ) -> float:
"""
Calculates the mass ratio of Helium with respect to the background gas.
-
+
Parameters
----------
abundances_gas: 9xN numpy array containing gas composition. Can be generated with Chemistry.abundances
abundances_solid: 9xN numpy array containing solid composition. Can be generated with Chemistry.abundances
nHe_nH: nHe/nH number density ratio of Helium and Hydrogen. This should match with your initial condition
T: Temperature array
-
+
Returns
-------
A float describing the mass ratio of Helium w.r.t. the background gas
-
+
"""
TotH_solid = split_molecules(abundances_solid)[-1]
molecule_elements = split_molecules(abundances_gas)
@@ -406,36 +438,58 @@ def calc_He_Mbg(self, abundances_gas: np.array, abundances_solid: np.array, nHe_
rest = np.sum(molecule_elements[:-1], axis=0)
rest_mol = abundances_gas[0]
Z = self.get_solid_heavies(T)
- Mdust_Mgas = Z*rest_mol # = Mmol_solid / Mbg * Mbg / Mgas = Mmol_solid/Mgas = Mdust/Mgas
- TotH_solid_Mgas = TotH_solid*Mdust_Mgas # Total mass of Hydrogen in solids over the total gas mass
-
+ Mdust_Mgas = (
+ Z * rest_mol
+ ) # = Mmol_solid / Mbg * Mbg / Mgas = Mmol_solid/Mgas = Mdust/Mgas
+ TotH_solid_Mgas = (
+ TotH_solid * Mdust_Mgas
+ ) # Total mass of Hydrogen in solids over the total gas mass
+
# This gives the total hydrogen mass in both gaseous and solid phase, normed to the gas mass, (Hgas+Hsolid)/gas
TotH = (1 - rest + TotH_solid_Mgas) / (1 + nHe_nH * MassHe)
# This gives Hgas/gas
TotH_gas = TotH - TotH_solid_Mgas
# This gives H/gas * He/H = He/gas
TotHe = TotH * nHe_nH * MassHe
-
+
# Some tests to check for some errors for the calculation of TotH and TotHe
- assert np.isclose(1-rest-molecule_elements[-1], rest_mol).all(), "Background gas from the molecular abundance array does not match background gas after splitting the molecules"
- assert np.isclose(TotH_gas + TotHe - TotH_mol_gas, rest_mol).all(), "Splitting the background gas into hydrogen and helium did not work"
- assert np.isclose(TotHe / (TotH*MassHe), nHe_nH).all(), "After separating solid and gaseous H, [He/H] is no longer correct"
- assert np.isclose(TotH_gas+TotHe+rest, 1).all(), "Adding up all gaseous components does not result in the total gas mass"
+ assert np.isclose(
+ 1 - rest - molecule_elements[-1], rest_mol
+ ).all(), "Background gas from the molecular abundance array does not match background gas after splitting the molecules"
+ assert np.isclose(
+ TotH_gas + TotHe - TotH_mol_gas, rest_mol
+ ).all(), "Splitting the background gas into hydrogen and helium did not work"
+ assert np.isclose(
+ TotHe / (TotH * MassHe), nHe_nH
+ ).all(), "After separating solid and gaseous H, [He/H] is no longer correct"
+ assert np.isclose(
+ TotH_gas + TotHe + rest, 1
+ ).all(), (
+ "Adding up all gaseous components does not result in the total gas mass"
+ )
-
# It is expected that He_Mbg is constant in space, so it's just saved as a float
assert np.isclose((TotHe / rest_mol).std(), 0), "He_Mbg is not constant."
- He_Mbg = (TotHe / rest_mol)[0]
+ He_Mbg = (TotHe / rest_mol)[0]
return He_Mbg
def get_solid_composition(self, sigma_dust_mol: np.array, T: np.array) -> np.array:
sigma_dust = np.sum(sigma_dust_mol, axis=1)
abundances_solid = np.transpose(
- np.divide(sigma_dust_mol, sigma_dust[:, np.newaxis], out=np.zeros_like(sigma_dust_mol),
- where=np.logical_and(sigma_dust[:, np.newaxis] > 1.0e-60, sigma_dust_mol > 1.0e-60)))
- chem_solid, _ = calc_composition(abundances_solid, solid=True, He_Mbg=self.He_Mbg)
+ np.divide(
+ sigma_dust_mol,
+ sigma_dust[:, np.newaxis],
+ out=np.zeros_like(sigma_dust_mol),
+ where=np.logical_and(
+ sigma_dust[:, np.newaxis] > 1.0e-60, sigma_dust_mol > 1.0e-60
+ ),
+ )
+ )
+ chem_solid, _ = calc_composition(
+ abundances_solid, solid=True, He_Mbg=self.He_Mbg
+ )
if np.max(T) > np.max(T_array):
- chem_solid[:np.argmax(T) + 1] = 0
+ chem_solid[: np.argmax(T) + 1] = 0
return chem_solid
@@ -460,29 +514,36 @@ def get_gas_composition(self, sigma_g_components_mol: np.array = None) -> np.arr
sigma = np.sum(sigma_g_components_mol, axis=1)
abundances_gas = np.transpose(
- np.divide(sigma_g_components_mol, sigma[:, np.newaxis], out=np.zeros_like(sigma_g_components_mol),
- where=sigma[:, np.newaxis] != 0))
+ np.divide(
+ sigma_g_components_mol,
+ sigma[:, np.newaxis],
+ out=np.zeros_like(sigma_g_components_mol),
+ where=sigma[:, np.newaxis] != 0,
+ )
+ )
chem_gas, _ = calc_composition(abundances_gas, solid=False, He_Mbg=self.He_Mbg)
self.mu = sigma / np.sum(sigma_g_components_mol.T / self.M_array, axis=0)
return chem_gas
def set_z(self, Z0):
- """ set the fraction of heavy elements in the disk """
+ """set the fraction of heavy elements in the disk"""
# The fraction of heavy !molecules! in either gas or solids with respect to the backgroundgas (H+He):
self.dtg0 = Z0
# ratio of heavies from the sum over all molecular species:
self._ZH = np.sum(self.abu_array[1:] * np.squeeze(self.M_array)[1:])
- rezHabu = MassH + self.HeHabu*MassHe + (self.elabu_init*elmasses).sum()
- self._dtg0_of_chem_model = self._ZH/rezHabu
+ rezHabu = MassH + self.HeHabu * MassHe + (self.elabu_init * elmasses).sum()
+ self._dtg0_of_chem_model = self._ZH / rezHabu
return
def get_solid_heavies(self, T):
- """ returns a modified dust to gas ratio according to the temperature of the disk """
+ """returns a modified dust to gas ratio according to the temperature of the disk"""
abu_array = np.ones_like(T) * self.abu_array[:, np.newaxis]
- M_masked_solid = np.where(T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array))
+ M_masked_solid = np.where(
+ T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array)
+ )
M_masked_solid[0, :] = 0 # exclude restmol
Tot_M_solid = np.sum(M_masked_solid * abu_array, axis=0) + 1e-300
@@ -508,60 +569,81 @@ def get_composition(self, T: np.array) -> np.array:
"""
abundances_solid, abundances_gas = self.abundances(T)
- self.mu = np.sum(abundances_gas, axis=0) / np.sum(abundances_gas / self.M_array, axis=0)
+ self.mu = np.sum(abundances_gas, axis=0) / np.sum(
+ abundances_gas / self.M_array, axis=0
+ )
# massratio tests for the elemental composition:
- assert np.all(np.isclose(np.sum(abundances_gas, axis=0),
- np.ones_like(T))), "error in chemistry, total sum of gas abundancies is not 1"
- assert np.all(np.isclose(np.sum(abundances_solid, axis=0),
- np.sum(abundances_solid, axis=0) != 0, np.ones_like(T),
- np.zeros_like(T))), "error in chemistry, total sum of solid abundancies is not 1"
+ assert np.all(
+ np.isclose(np.sum(abundances_gas, axis=0), np.ones_like(T))
+ ), "error in chemistry, total sum of gas abundancies is not 1"
+ assert np.all(
+ np.isclose(
+ np.sum(abundances_solid, axis=0),
+ np.sum(abundances_solid, axis=0) != 0,
+ np.ones_like(T),
+ np.zeros_like(T),
+ )
+ ), "error in chemistry, total sum of solid abundancies is not 1"
self.He_Mbg = self.calc_He_Mbg(abundances_gas, abundances_solid, self.HeHabu, T)
- chem_solid, _ = calc_composition(abundances_solid, solid=True, He_Mbg=self.He_Mbg)
+ chem_solid, _ = calc_composition(
+ abundances_solid, solid=True, He_Mbg=self.He_Mbg
+ )
chem_gas, _ = calc_composition(abundances_gas, solid=False, He_Mbg=self.He_Mbg)
if np.max(T) > np.max(T_array):
- chem_solid[:np.argmax(T) + 1] = 0
+ chem_solid[: np.argmax(T) + 1] = 0
else:
# massratio tests for the molecular composition of the solids:
# only test, when there are solids!
- assert np.all(np.isclose(np.sum(chem_solid[:, 0, :], axis=1),
- np.ones_like(T))), "error in chemistry, total sum of solid abundancies is not 1"
+ assert np.all(
+ np.isclose(np.sum(chem_solid[:, 0, :], axis=1), np.ones_like(T))
+ ), "error in chemistry, total sum of solid abundancies is not 1"
# massratio tests for the molecular composition of the gas:
- assert np.all(np.isclose(np.sum(chem_gas[:, 0, :], axis=1),
- np.ones_like(T))), "error in chemistry, total sum of gas abundancies is not 1"
+ assert np.all(
+ np.isclose(np.sum(chem_gas[:, 0, :], axis=1), np.ones_like(T))
+ ), "error in chemistry, total sum of gas abundancies is not 1"
###############
# Sanity checks
###############
TotH_mol_test = (
- self.abuCH4 * 4.0
- + self.abuH2O * 2.0
- + self.abuNH3 * 3
- + self.abuH2S * 2
+ self.abuCH4 * 4.0 + self.abuH2O * 2.0 + self.abuNH3 * 3 + self.abuH2S * 2
)
abu_mol_test = (self.abu_array * self.M_array.squeeze()).sum() - TotH_mol_test
abu_asp_test = (self.elabu_init * elmasses).sum()
- assert np.isclose(abu_asp_test, abu_mol_test), "error in chemistry, total sum of molecular abundancies is not equal to the total sum of the X/H stellar vaules. Partioning is wrong."
+ assert np.isclose(
+ abu_asp_test, abu_mol_test
+ ), "error in chemistry, total sum of molecular abundancies is not equal to the total sum of the X/H stellar vaules. Partioning is wrong."
# Test if the heavy elements are identical to the heavy molecules (gas):
- _, abu_gas_test = self.abundances([np.max(T_array)+100])
- chem_gas_test, restH = calc_composition(abu_gas_test, solid=False, He_Mbg=self.He_Mbg)
- elabus_heavy_chem_test = np.sum(chem_gas_test[0,0,:-7])
- assert np.isclose(chem_gas_test[0, 1, 1:].sum(), elabus_heavy_chem_test+restH), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies."
+ _, abu_gas_test = self.abundances([np.max(T_array) + 100])
+ chem_gas_test, restH = calc_composition(
+ abu_gas_test, solid=False, He_Mbg=self.He_Mbg
+ )
+ elabus_heavy_chem_test = np.sum(chem_gas_test[0, 0, :-7])
+ assert np.isclose(
+ chem_gas_test[0, 1, 1:].sum(), elabus_heavy_chem_test + restH
+ ), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies."
# Test validity of dust to gas ratio:
- assert np.isclose(chem_gas_test[0,1,1:].sum()*(1+self.dtg0),self.dtg0), "Error in dust to gas ratio"
+ assert np.isclose(
+ chem_gas_test[0, 1, 1:].sum() * (1 + self.dtg0), self.dtg0
+ ), "Error in dust to gas ratio"
# Test if the heavy elements are identical to the heavy molecules (solids):
if np.max(T) < np.max(T_array):
# only test, when there are solids!
abu_solid_test, _ = self.abundances([0])
- chem_solid_test, restH = calc_composition(abu_solid_test, solid=True, He_Mbg=self.He_Mbg)
- elabus_heavy_chem_test = np.sum(chem_solid_test[0,0,:-7])
- assert np.isclose(chem_solid_test[0, 1, 1:].sum(), elabus_heavy_chem_test+restH), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies."
+ chem_solid_test, restH = calc_composition(
+ abu_solid_test, solid=True, He_Mbg=self.He_Mbg
+ )
+ elabus_heavy_chem_test = np.sum(chem_solid_test[0, 0, :-7])
+ assert np.isclose(
+ chem_solid_test[0, 1, 1:].sum(), elabus_heavy_chem_test + restH
+ ), "error in chemistry, total sum of elemental mass abundancies is not equal to the total sum of molecular mass abundancies."
return chem_solid, chem_gas
@@ -582,19 +664,27 @@ def abundances(self, T: np.array) -> Tuple[np.array, np.array]:
"""
abu_array = np.ones_like(T) * self.abu_array[:, np.newaxis]
- M_masked_solid = np.where(T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array))
+ M_masked_solid = np.where(
+ T < T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array)
+ )
M_masked_solid[0, :] = 0 # exclude restmol
Tot_M_solid = np.sum(M_masked_solid * abu_array, axis=0) + 1e-300
- Masses_solid = np.where(T < T_array, abu_array * M_masked_solid / Tot_M_solid,
- np.zeros_like(abu_array))
+ Masses_solid = np.where(
+ T < T_array,
+ abu_array * M_masked_solid / Tot_M_solid,
+ np.zeros_like(abu_array),
+ )
Z = self.get_solid_heavies(T)
- M_masked_gas = np.where(T >= T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array))
+ M_masked_gas = np.where(
+ T >= T_array, np.ones_like(T) * self.M_array, np.zeros_like(abu_array)
+ )
M_masked_gas[0, :] = 0 # exclude restmol
Tot_M_gas = np.sum(M_masked_gas * abu_array, axis=0) + 1e-300
- Masses_gas = np.where(T >= T_array, abu_array * M_masked_gas / Tot_M_gas,
- np.zeros_like(abu_array))
+ Masses_gas = np.where(
+ T >= T_array, abu_array * M_masked_gas / Tot_M_gas, np.zeros_like(abu_array)
+ )
Masses_gas[0, :] = 1 / (1 + (self.dtg0 - Z))
Masses_gas[1:, :] = Masses_gas[1:, :] * (self.dtg0 - Z) / (1 + (self.dtg0 - Z))
diff --git a/chemcomp/disks/_disk_class.py b/chemcomp/disks/_disk_class.py
index dc45bab..75cd1d9 100755
--- a/chemcomp/disks/_disk_class.py
+++ b/chemcomp/disks/_disk_class.py
@@ -17,17 +17,19 @@
AU = const.au.cgs.value
ME = const.M_earth.cgs.value
-FWHM_TO_SIGMA = 2*(2*np.log(2))**0.5
+FWHM_TO_SIGMA = 2 * (2 * np.log(2)) ** 0.5
+
def gaussian(x, max, mu, FWHM):
- """Gaussian helper function. """
+ """Gaussian helper function."""
# prevent underflow by maximizing exponent to -1e2 (np.exp(-1e2)=3e-44)
- exp = np.minimum(0.5*((x-mu)/FWHM*FWHM_TO_SIGMA)**2.0, 1e2)
- return max*np.exp(-exp)
+ exp = np.minimum(0.5 * ((x - mu) / FWHM * FWHM_TO_SIGMA) ** 2.0, 1e2)
+ return max * np.exp(-exp)
+
def tvisc_fct(T, sigma_g, tirr, omega_k, alpha, Z, mu):
"""Calculate viscous heating and return the difference between our current guess of the temperature
- to the last guess of the temperature """
+ to the last guess of the temperature"""
c_s = (k_B * T / (mu * m_p)) ** 0.5
h = c_s / omega_k
rho_g = sigma_g / (np.sqrt(2 * pi) * h)
@@ -36,28 +38,31 @@ def tvisc_fct(T, sigma_g, tirr, omega_k, alpha, Z, mu):
opacity = opacity_fct(rho_g, T, Z=Z) * 100.0 * Z
- qvisc = (9.0 / 4.0) * sigma_g * nu * omega_k ** 2.0
+ qvisc = (9.0 / 4.0) * sigma_g * nu * omega_k**2.0
tau = sigma_g * opacity # both sides of the disk
- tvisc4 = qvisc * 3.0 * tau / (16.0 * sr) # factor 0.5 because tau is defined on complete disk
+ tvisc4 = (
+ qvisc * 3.0 * tau / (16.0 * sr)
+ ) # factor 0.5 because tau is defined on complete disk
- T_new = np.clip((tvisc4 + tirr ** 4) ** 0.25, tirr, 1.5e5)
+ T_new = np.clip((tvisc4 + tirr**4) ** 0.25, tirr, 1.5e5)
diff = T - T_new
return diff
+
def solve_viscous_heating_globally(
- sigma_g,
- tirr,
- T,
- omega_k,
- alpha,
- Z,
- mu,
- r,
- xtol=1e-1,
- nitermax=100,
- interpolation_idx=None
+ sigma_g,
+ tirr,
+ T,
+ omega_k,
+ alpha,
+ Z,
+ mu,
+ r,
+ xtol=1e-1,
+ nitermax=100,
+ interpolation_idx=None,
):
"""
Solve the midplane tempeature due to viscous heating (only)
@@ -68,44 +73,65 @@ def solve_viscous_heating_globally(
sample = np.linspace(0.1, 300.0, 10000) * AU
- T = np.array([optimize.root_scalar(
- tvisc_fct,
- bracket=(tirr[i], 1.5e4),
- x0=T[i], # estimate T from last iteration
- args=(sigma_g[i], tirr[i], omega_k[i], alpha[i] if hasattr(alpha, "__iter__") else alpha, Z[i] if hasattr(Z, "__iter__") else Z, mu[i] if hasattr(mu, "__iter__") else mu),
- options={"maxiter": nitermax, "xtol": xtol}
- ).root for i in range(len(T)) if interpolation_idx[i]])
- T_sample = savgol_filter(interp1d(r[interpolation_idx], T, fill_value="extrapolate", assume_sorted=True)(sample), 5, 2)
+ T = np.array(
+ [
+ optimize.root_scalar(
+ tvisc_fct,
+ bracket=(tirr[i], 1.5e4),
+ x0=T[i], # estimate T from last iteration
+ args=(
+ sigma_g[i],
+ tirr[i],
+ omega_k[i],
+ alpha[i] if hasattr(alpha, "__iter__") else alpha,
+ Z[i] if hasattr(Z, "__iter__") else Z,
+ mu[i] if hasattr(mu, "__iter__") else mu,
+ ),
+ options={"maxiter": nitermax, "xtol": xtol},
+ ).root
+ for i in range(len(T))
+ if interpolation_idx[i]
+ ]
+ )
+ T_sample = savgol_filter(
+ interp1d(r[interpolation_idx], T, fill_value="extrapolate", assume_sorted=True)(
+ sample
+ ),
+ 5,
+ 2,
+ )
T = interp1d(sample, T_sample, fill_value="extrapolate", assume_sorted=True)(r)
T = np.clip(T, tirr, 1.5e5)
return T
class Disk(object):
- """ Baseclass of the disk. """
+ """Baseclass of the disk."""
def __init__(
- self,
- defaults,
- M_STAR,
- ALPHA,
- chemistry,
- init_grid=True,
- time=0 * u.yr,
- **kwargs
+ self,
+ defaults,
+ M_STAR,
+ ALPHA,
+ chemistry,
+ init_grid=True,
+ time=0 * u.yr,
+ **kwargs,
):
- self.sigmin = 1e-60 # sigma_floor value
- self.defaults = defaults # dictionary containing the defaults section of the config
+ self.sigmin = 1e-60 # sigma_floor value
+ self.defaults = (
+ defaults # dictionary containing the defaults section of the config
+ )
self.M_s = M_STAR.cgs.value # stellar mass
self.alpha = float(ALPHA) # alpha viscosity, fix yaml issue
- self._chem = chemistry # reference to chemistry object
+ self._chem = chemistry # reference to chemistry object
self.t_0 = 1 * time.cgs.value # starting time of the disk
self.t = 1 * time.cgs.value # time in the disk (initialised with starting time)
# init timestep
self._timestep_input = defaults.get("DEF_TIMESTEP", None)
if self._timestep_input is not None:
- print(f'WARNING: We are using a custom timestep of {self._timestep_input}')
+ print(f"WARNING: We are using a custom timestep of {self._timestep_input}")
self._timestep_input = self._timestep_input.cgs.value
# init t_end
@@ -116,20 +142,27 @@ def __init__(
r_in = defaults.get("DEF_R_IN", 0.2 * u.au) # inner grid edge
r_out = defaults.get("DEF_R_OUT", 100.0 * u.au) # outer grid edge
- self.gridsize = defaults.get("DEF_GRIDSIZE", 1000) # number of grid cells
+ self.gridsize = defaults.get("DEF_GRIDSIZE", 1000) # number of grid cells
if defaults.get("DEF_LIN_SPACING", True):
- self.r = np.linspace(r_in, r_out, self.gridsize).to(u.au).value # linear grid
+ self.r = (
+ np.linspace(r_in, r_out, self.gridsize).to(u.au).value
+ ) # linear grid
else:
- self.r = ( # log grid
- np.logspace(np.log10(r_in.to(u.au).value), np.log10(r_out.to(u.au).value), self.gridsize)
+ self.r = np.logspace( # log grid
+ np.log10(r_in.to(u.au).value),
+ np.log10(r_out.to(u.au).value),
+ self.gridsize,
)
- self.r_i = np.array([self.r[0], *((self.r[1:] + self.r[:-1]) / 2.0), self.r[-1]]) * u.au
+ self.r_i = (
+ np.array([self.r[0], *((self.r[1:] + self.r[:-1]) / 2.0), self.r[-1]])
+ * u.au
+ )
self.r = self.r * u.au
self.r_i = self.r_i.cgs.value # interfacepositions of the radial grid
- self.r = self.r.cgs.value # cellcenter of the radial grid
+ self.r = self.r.cgs.value # cellcenter of the radial grid
self._init_params(**kwargs)
@@ -149,45 +182,79 @@ def _init_params(self, **kwargs):
# array holding the dust to gas ratio
self.DTG = {"total": eval_kwargs(kwargs.pop("DTG_total", 0.015))}
- self.DTG_small_grains = self.DTG["total"] * (1.0-0.75) # DTG small grains is only used for the temperature
- self.conf_static = eval_kwargs(kwargs.get("static", False)) # config for static disk
- self.conf_static_stokes = eval_kwargs(kwargs.get("static_stokes", False)) # config for static stokes number
- self.conf_evaporation = eval_kwargs(kwargs.get("evaporation", True)) # config for static stokes number
- self.conf_temp_evol = eval_kwargs(kwargs.get("temp_evol", False)) # use the temperature evolution. Advise: dont do!
-
- self._evap_width = eval_kwargs(kwargs.get("evap_width", 0.1*u.au)).cgs.value # config for static stokes number
-
- self.alpha_height = eval_kwargs(kwargs.get("ALPHAHEIGHT", self.alpha)) # vertical mixing alpha
- self.alpha_frag = eval_kwargs(kwargs.get("ALPHAFRAG", self.alpha)) # fragmentation alpha, not used in Schneider & Bitsch 2021
- self.alpha_mig = eval_kwargs(kwargs.get("ALPHAMIG", self.alpha)) # migration alpha, not used in Schneider & Bitsch 2021
+ self.DTG_small_grains = self.DTG["total"] * (
+ 1.0 - 0.75
+ ) # DTG small grains is only used for the temperature
+ self.conf_static = eval_kwargs(
+ kwargs.get("static", False)
+ ) # config for static disk
+ self.conf_static_stokes = eval_kwargs(
+ kwargs.get("static_stokes", False)
+ ) # config for static stokes number
+ self.conf_evaporation = eval_kwargs(
+ kwargs.get("evaporation", True)
+ ) # config for static stokes number
+ self.conf_temp_evol = eval_kwargs(
+ kwargs.get("temp_evol", False)
+ ) # use the temperature evolution. Advise: dont do!
+
+ self._evap_width = eval_kwargs(
+ kwargs.get("evap_width", 0.1 * u.au)
+ ).cgs.value # config for static stokes number
+
+ self.alpha_height = eval_kwargs(
+ kwargs.get("ALPHAHEIGHT", self.alpha)
+ ) # vertical mixing alpha
+ self.alpha_frag = eval_kwargs(
+ kwargs.get("ALPHAFRAG", self.alpha)
+ ) # fragmentation alpha, not used in Schneider & Bitsch 2021
+ self.alpha_mig = eval_kwargs(
+ kwargs.get("ALPHAMIG", self.alpha)
+ ) # migration alpha, not used in Schneider & Bitsch 2021
self._chem.set_z(self.DTG.get("total", 0.015)) # set the solid to gas ratio
if self.conf_static_stokes:
self.f_m = np.ones_like(self.r) # Massdistribution factor for twopop model
- self.stokes_number_small = np.zeros_like(self.r) # stokes number of the small population for twopop model
- self.stokes_number_df = np.zeros_like(self.r) # stokes number of the drift limited fragmentation for twopop model
- self.stokes_number_drift = np.zeros_like(self.r) # stokes number of the drift limited case for twopop model
- self.stokes_number_frag = np.zeros_like(self.r) # stokes number of the fragmentation limited case for twopop model
+ self.stokes_number_small = np.zeros_like(
+ self.r
+ ) # stokes number of the small population for twopop model
+ self.stokes_number_df = np.zeros_like(
+ self.r
+ ) # stokes number of the drift limited fragmentation for twopop model
+ self.stokes_number_drift = np.zeros_like(
+ self.r
+ ) # stokes number of the drift limited case for twopop model
+ self.stokes_number_frag = np.zeros_like(
+ self.r
+ ) # stokes number of the fragmentation limited case for twopop model
self.const_f_f = 0.37 # fudge fitting factors from Birnstiel2012
self.const_f_d = 0.55 # fudge fitting factors from Birnstiel2012
- self.const_N = 0.5 # N of Birnstiel2012
- self.a_0 = eval_kwargs(kwargs.get("a_0", 1e-4 * u.cm)) # grainsize of the small population
- self.a_0 = self.a_0.cgs.value # grainsize of the small population
- self.a_1 = eval_kwargs(kwargs.get("a_0", 1e-4 * u.cm)) * np.ones_like(self.r) # grainsize of the large population
- self.a_1 = self.a_1.cgs.value # grainsize of the large population
-
- self.omega_k = np.sqrt(G * self.M_s / self.r ** 3) # kepler orbitalperiod
- self.v_k = self.omega_k * self.r # kepler velocity
- self.alpha_factor = np.ones_like(self.r) # gap factor that applies the gap to the alphaviscosity
+ self.const_N = 0.5 # N of Birnstiel2012
+ self.a_0 = eval_kwargs(
+ kwargs.get("a_0", 1e-4 * u.cm)
+ ) # grainsize of the small population
+ self.a_0 = self.a_0.cgs.value # grainsize of the small population
+ self.a_1 = eval_kwargs(kwargs.get("a_0", 1e-4 * u.cm)) * np.ones_like(
+ self.r
+ ) # grainsize of the large population
+ self.a_1 = self.a_1.cgs.value # grainsize of the large population
+
+ self.omega_k = np.sqrt(G * self.M_s / self.r**3) # kepler orbitalperiod
+ self.v_k = self.omega_k * self.r # kepler velocity
+ self.alpha_factor = np.ones_like(
+ self.r
+ ) # gap factor that applies the gap to the alphaviscosity
# Init for output
self.vr_dust = np.zeros_like(self.r[1:]) # dust velocity
- self.vr_gas = np.zeros_like(self.r[1:]) # gas velocity
+ self.vr_gas = np.zeros_like(self.r[1:]) # gas velocity
- tau_disk = eval_kwargs(kwargs.get("tau_disk", None)) # disk dispersal time
- self.begin_photevap = eval_kwargs(kwargs.get("begin_photoevap", 0 * u.Myr).cgs.value) # time until disk disappears
+ tau_disk = eval_kwargs(kwargs.get("tau_disk", None)) # disk dispersal time
+ self.begin_photevap = eval_kwargs(
+ kwargs.get("begin_photoevap", 0 * u.Myr).cgs.value
+ ) # time until disk disappears
if tau_disk is not None:
self.conf_photoevaporation = True
@@ -199,45 +266,59 @@ def _init_params(self, **kwargs):
# init physical quantities
# calculate the chemical composition of sigma:
- self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T) # chemistry vectors
+ self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(
+ self.T
+ ) # chemistry vectors
dtg = self._chem.get_solid_heavies(self.T)
- self.sigma_g_components = self.chemistry_gas * self.sigma_g[:,np.newaxis,np.newaxis] # gas surface density as a compositional vector
- self.sigma_g = np.sum(self.sigma_g_components[:, 0, :], axis=1) # gas surface density
-
- self.rho_g = self.sigma_g / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r) # gas volume density
- self.grad_sig = np.gradient( # gradient of the surfacedensity
+ self.sigma_g_components = (
+ self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ ) # gas surface density as a compositional vector
+ self.sigma_g = np.sum(
+ self.sigma_g_components[:, 0, :], axis=1
+ ) # gas surface density
+
+ self.rho_g = self.sigma_g / (
+ np.sqrt(2 * np.pi) * self.aspect_ratio * self.r
+ ) # gas volume density
+ self.grad_sig = np.gradient( # gradient of the surfacedensity
np.log10(self.sigma_g), np.log10(self.r)
) # -1.5 for MMSN
- self.c_s = self.aspect_ratio * self.r * self.omega_k # soundspeed
- self.P = ( # pressure
- self.c_s ** 2
- * self.sigma_g
- / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r)
+ self.c_s = self.aspect_ratio * self.r * self.omega_k # soundspeed
+ self.P = ( # pressure
+ self.c_s**2
+ * self.sigma_g
+ / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r)
)
- self.grad_p = np.gradient( # pressuregradient
+ self.grad_p = np.gradient( # pressuregradient
np.log10(self.P), np.log10(self.r)
) # -2.75 for MMSN, checked
- self.T = 2.35 * self.c_s ** 2.0 * m_p / k_B # Temperature
+ self.T = 2.35 * self.c_s**2.0 * m_p / k_B # Temperature
if not hasattr(self, "T_irr"):
- self.T_irr = 1.0 * self.T # effective Temperature from irradiation
+ self.T_irr = 1.0 * self.T # effective Temperature from irradiation
if not hasattr(self, "T_visc"):
- self.T_visc = np.zeros_like(self.T) # viscous Temperature from viscous heating
+ self.T_visc = np.zeros_like(
+ self.T
+ ) # viscous Temperature from viscous heating
self.compute_cs_and_hp()
self.compute_nu()
self.calc_opacity()
- self.grad_T = np.gradient( # Temperature gradient
+ self.grad_T = np.gradient( # Temperature gradient
np.log10(self.T), np.log10(self.r)
) # -0.5 for MMSN
- self.eta = -0.5 * self.aspect_ratio ** 2 * self.grad_p
+ self.eta = -0.5 * self.aspect_ratio**2 * self.grad_p
enclosed_mask = self.r < 200.0 * AU
- self.disk_mass = np.trapz(2 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask],
- self.r[enclosed_mask])
+ self.disk_mass = np.trapz(
+ 2 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask],
+ self.r[enclosed_mask],
+ )
self.disk_mass_dust = 0
- self.m_dot_components = self.compute_mdot_at_interfaces(self.sigma_g_components)
+ self.m_dot_components = self.compute_mdot_at_interfaces(
+ self.sigma_g_components
+ )
self.m_dot = np.sum(self.m_dot_components[:, 0, :], axis=1)
self.pebble_flux = np.zeros_like(self.r)
@@ -248,8 +329,10 @@ def _init_params(self, **kwargs):
print("caution: Unstable disk: Q_min = {}".format(np.min(self.Q)))
def calc_pebble_iso(self, diffusion=False):
- """Calculate the pebble isolation mass. """
- f_fit = (self.aspect_ratio / 0.05) ** 3.0 * (0.34 * (np.log10(1e-3) / np.log10(self.alpha)) ** 4.0 + 0.66)
+ """Calculate the pebble isolation mass."""
+ f_fit = (self.aspect_ratio / 0.05) ** 3.0 * (
+ 0.34 * (np.log10(1e-3) / np.log10(self.alpha)) ** 4.0 + 0.66
+ )
if diffusion:
self.peb_iso_pi_crit = self.alpha / (2.0 * self.stokes_number_pebbles)
@@ -257,7 +340,9 @@ def calc_pebble_iso(self, diffusion=False):
else:
self.peb_iso = 25.0 * ME * f_fit
- def init_pebble_parameters(self, epsilon_p, rho_solid, t_f, u_frag, calculate_rho_solid, calculate_u_frag):
+ def init_pebble_parameters(
+ self, epsilon_p, rho_solid, t_f, u_frag, calculate_rho_solid, calculate_u_frag
+ ):
"""
Sets the pebble parameters needed for pebble evaporation model
all in cgs
@@ -275,7 +360,7 @@ def init_pebble_parameters(self, epsilon_p, rho_solid, t_f, u_frag, calculate_rh
self.epsilon_p = epsilon_p
self.rho_solid = rho_solid
if self.rho_solid is not None:
- self.rho_solid = np.ones_like(self.sigma_dust)*self.rho_solid
+ self.rho_solid = np.ones_like(self.sigma_dust) * self.rho_solid
self._calculate_rho_solid = calculate_rho_solid
self._calculate_u_frag = calculate_u_frag
@@ -313,23 +398,62 @@ def calc_sigma_dust(self):
self.solid_fraction = self._chem.get_solid_heavies(self.T)
self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T)
- sigma_g0 = self.sigma_g / (1.0 + (self._chem.dtg0-self.solid_fraction)) # get the fraction of the background gas
-
- self.sigma_dust_components = self.chemistry_solid * self.solid_fraction[:, np.newaxis, np.newaxis] * sigma_g0[:, np.newaxis,
- np.newaxis]
- self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ sigma_g0 = self.sigma_g / (
+ 1.0 + (self._chem.dtg0 - self.solid_fraction)
+ ) # get the fraction of the background gas
+
+ self.sigma_dust_components = (
+ self.chemistry_solid
+ * self.solid_fraction[:, np.newaxis, np.newaxis]
+ * sigma_g0[:, np.newaxis, np.newaxis]
+ )
+ self.sigma_g_components = (
+ self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ )
self.sigma_dust = np.sum(self.sigma_dust_components[:, 1, :], axis=1)
self.sigma_dust_components[:, 1, 0] = np.zeros_like(self.sigma_g)
- assert np.all(np.isclose((self.sigma_g + self.sigma_dust), sigma_g0 * (1.0 + self._chem.dtg0))), "error in chemistry: The total masses do not sum up correctly"
- assert np.all(np.isclose(((self.sigma_dust_components[:,1,:]+self.sigma_g_components[:,1,:])/(self.sigma_g+self.sigma_dust)[:,np.newaxis]).std(axis=0),0)), "error in chemistry: The individual masses do not sum up to a constant value"
- assert np.all(np.isclose(np.sum(self.sigma_dust_components[:, 1, :] + self.sigma_g_components[:, 1, :], axis=1) / sigma_g0,
- (1+self._chem.dtg0))), "error in chemistry: the dust to gas ratio is not correct"
- assert np.all(np.isclose(
- np.sum(self.sigma_dust_components[:, 1, 1:] + self.sigma_g_components[:, 1, 1:], axis=1) / sigma_g0,
- self._chem.dtg0)), "error in chemistry: the dust to gas ratio is not correct"
+ assert np.all(
+ np.isclose(
+ (self.sigma_g + self.sigma_dust), sigma_g0 * (1.0 + self._chem.dtg0)
+ )
+ ), "error in chemistry: The total masses do not sum up correctly"
+ assert np.all(
+ np.isclose(
+ (
+ (
+ self.sigma_dust_components[:, 1, :]
+ + self.sigma_g_components[:, 1, :]
+ )
+ / (self.sigma_g + self.sigma_dust)[:, np.newaxis]
+ ).std(axis=0),
+ 0,
+ )
+ ), "error in chemistry: The individual masses do not sum up to a constant value"
+ assert np.all(
+ np.isclose(
+ np.sum(
+ self.sigma_dust_components[:, 1, :]
+ + self.sigma_g_components[:, 1, :],
+ axis=1,
+ )
+ / sigma_g0,
+ (1 + self._chem.dtg0),
+ )
+ ), "error in chemistry: the dust to gas ratio is not correct"
+ assert np.all(
+ np.isclose(
+ np.sum(
+ self.sigma_dust_components[:, 1, 1:]
+ + self.sigma_g_components[:, 1, 1:],
+ axis=1,
+ )
+ / sigma_g0,
+ self._chem.dtg0,
+ )
+ ), "error in chemistry: the dust to gas ratio is not correct"
self.sigmin_peb = 1e-60
@@ -349,18 +473,25 @@ def remove_from_gas(self, mdot, f_red, idx):
"""
try:
sum_f_red = np.sum(f_red)
- delta_sigma = f_red / (
- (self.r_i[1:][idx] - self.r_i[:-1][idx]) * sum_f_red)
+ delta_sigma = f_red / ((self.r_i[1:][idx] - self.r_i[:-1][idx]) * sum_f_red)
if hasattr(self, "sigdot"): # remove using visc disk evolution
self.sigdot_gas_acc = np.zeros_like(self.sigdot)
self.sigdot_gas_acc[idx] = delta_sigma[:, np.newaxis] * mdot[1]
else: # there is no visc disk evolution, remove directly instead
sum_mdot_el = np.sum(mdot, axis=1)[0]
- self.sigma_g[idx] -= sum_mdot_el * delta_sigma / (2 * np.pi * self.r[idx]) * self.dt
- self.sigma_g_components[idx] -= delta_sigma[:, np.newaxis, np.newaxis] / (
- 2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis]) * mdot * self.dt
+ self.sigma_g[idx] -= (
+ sum_mdot_el * delta_sigma / (2 * np.pi * self.r[idx]) * self.dt
+ )
+ self.sigma_g_components[idx] -= (
+ delta_sigma[:, np.newaxis, np.newaxis]
+ / (2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis])
+ * mdot
+ * self.dt
+ )
self.sigma_g[self.sigma_g < self.sigmin] = self.sigmin
- self.sigma_g_components[self.sigma_g_components < self.sigmin] = self.sigmin
+ self.sigma_g_components[
+ self.sigma_g_components < self.sigmin
+ ] = self.sigmin
return
except IndexError:
return
@@ -382,23 +513,30 @@ def remove_from_peb(self, mdot, idx):
self.sigdot_peb_acc[idx] = delta_sigma * mdot[1] # fill sigdot
else:
mdot_sum = np.sum(mdot, axis=1)[0]
- self.sigma_dust[idx] -= delta_sigma / (2.0 * np.pi * self.r[idx]) * mdot_sum * self.dt
- self.sigma_dust_components[idx] -= delta_sigma / (
- 2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis]) * mdot * self.dt
+ self.sigma_dust[idx] -= (
+ delta_sigma / (2.0 * np.pi * self.r[idx]) * mdot_sum * self.dt
+ )
+ self.sigma_dust_components[idx] -= (
+ delta_sigma
+ / (2.0 * np.pi * self.r[idx, np.newaxis, np.newaxis])
+ * mdot
+ * self.dt
+ )
self.sigma_dust[self.sigma_dust < self.sigmin_peb] = self.sigmin_peb
self.sigma_dust_components[
- self.sigma_dust_components < self.sigmin_peb] = self.sigmin_peb
+ self.sigma_dust_components < self.sigmin_peb
+ ] = self.sigmin_peb
return
except IndexError:
return
def apply_gap_profile(self, a_p, fsigma, FWHM):
- #fsigma = np.maximum(fsigma, 1e-7)
+ # fsigma = np.maximum(fsigma, 1e-7)
self.alpha_factor = 1 - gaussian(self.r, 1.0 - fsigma, a_p, FWHM)
- self.interpolation_idx = np.abs(self.alpha_factor - 1.0 ) < 0.01
+ self.interpolation_idx = np.abs(self.alpha_factor - 1.0) < 0.01
self.interpolation_idx[0] = True
- alpha_factor_interface = 0.5*(self.alpha_factor[1:]+self.alpha_factor[:-1])
+ alpha_factor_interface = 0.5 * (self.alpha_factor[1:] + self.alpha_factor[:-1])
midinterface = np.abs(alpha_factor_interface - 1) < 0.01
self.interpolation_idx_interface = np.array([True, *midinterface, True])
@@ -428,23 +566,38 @@ def compute_sizes(self):
msil = np.sum(self.sigma_dust_components[:, 1, 8:], axis=-1)
mice = np.sum(self.sigma_dust_components[:, 1, 1:8], axis=-1)
# assert np.all(np.isclose(msil + mice, self.sigma_dust)), 'The internal ice densities do not add up correctly'
- self.rho_solid = (msil+mice) / (msil/DENSITY_SILICATES + mice/DENSITY_ICE + 1e-90)
+ self.rho_solid = (msil + mice) / (
+ msil / DENSITY_SILICATES + mice / DENSITY_ICE + 1e-90
+ )
self.rho_solid = np.maximum(self.rho_solid, DENSITY_ICE)
if self._calculate_u_frag:
# Calculate the fragmentation velocity according to the Drazkowska 2017 approach: 1m/s if silicate, 10m/s if icy
mice = np.sum(self.sigma_dust_components[:, 1, 1:8], axis=-1)
- self.u_frag = np.where(mice > 0.01*self.sigma_dust, 1000, 100)
-
- a_frag = self.const_f_f * 2.0 / (3.0 * np.pi) * self.sigma_g / (
- self.alpha_frag * self.rho_solid) * self.u_frag ** 2.0 / (
- self.c_s ** 2.0)
+ self.u_frag = np.where(mice > 0.01 * self.sigma_dust, 1000, 100)
+
+ a_frag = (
+ self.const_f_f
+ * 2.0
+ / (3.0 * np.pi)
+ * self.sigma_g
+ / (self.alpha_frag * self.rho_solid)
+ * self.u_frag**2.0
+ / (self.c_s**2.0)
+ )
tau_grow = self.sigma_g / ((self.omega_k * self.sigma_dust) + 1e-300)
dlnpdlnr = np.abs(self.grad_p)
- a_drift = self.const_f_d * 2 * self.sigma_dust / (
- np.pi * self.rho_solid) * self.v_k ** 2 / self.c_s ** 2 / dlnpdlnr
- St_df = self.u_frag * self.v_k / (dlnpdlnr * self.c_s ** 2 * (1 - self.const_N))
+ a_drift = (
+ self.const_f_d
+ * 2
+ * self.sigma_dust
+ / (np.pi * self.rho_solid)
+ * self.v_k**2
+ / self.c_s**2
+ / dlnpdlnr
+ )
+ St_df = self.u_frag * self.v_k / (dlnpdlnr * self.c_s**2 * (1 - self.const_N))
a_df = St_df * 2 * self.sigma_g / (np.pi * self.rho_solid)
a_1 = np.min([a_drift, a_frag, a_df], axis=0)
@@ -461,9 +614,11 @@ def compute_sizes(self):
self.stokes_number_pebbles = a_1 * stokes_factor
self.stokes_number_small = self.a_0 * stokes_factor
- self.f_m = np.where(np.argmin([a_drift, a_frag, a_df], axis=0) == 0,
- 0.97 * np.ones_like(self.stokes_number_pebbles),
- 0.75 * np.ones_like(self.stokes_number_pebbles))
+ self.f_m = np.where(
+ np.argmin([a_drift, a_frag, a_df], axis=0) == 0,
+ 0.97 * np.ones_like(self.stokes_number_pebbles),
+ 0.75 * np.ones_like(self.stokes_number_pebbles),
+ )
self.stokes_number_df = 1.0 * St_df
self.stokes_number_drift = stokes_factor * a_drift
@@ -494,8 +649,11 @@ def compute_viscous_evolution(self):
if self.conf_photoevaporation:
self.calc_sigdot_photoevap()
- sigma_g_components_mol = self.get_viscous_evolution_next_timestep_components(self.sigma_g_components,
- self.dt)
+ sigma_g_components_mol = (
+ self.get_viscous_evolution_next_timestep_components(
+ self.sigma_g_components, self.dt
+ )
+ )
self.compute_dust_next_timestep(self.dt)
self.compute_sigma_g_from_mol(sigma_g_components_mol)
@@ -503,19 +661,28 @@ def compute_viscous_evolution(self):
else: # init sigma_g_components
self.compute_disktmid()
dtg = self._chem.get_solid_heavies(self.T)
- self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T)
+ self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(
+ self.T
+ )
self.mu = self._chem.mu
- self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
-
+ self.sigma_g_components = (
+ self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ )
self.compute_disktmid()
if not self.conf_evaporation:
- self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(self.T)
+ self.chemistry_solid, self.chemistry_gas = self._chem.get_composition(
+ self.T
+ )
self.mu = self._chem.mu
- self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis]
- self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ self.sigma_dust_components = (
+ self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis]
+ )
+ self.sigma_g_components = (
+ self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ )
if hasattr(self, "vr_peb"):
self.calc_pebble_flux(noupdate=False)
@@ -525,32 +692,32 @@ def compute_viscous_evolution(self):
self._old_disk_mass = 1.0 * self.disk_mass
self._old_disk_mass_dust = 1.0 * self.disk_mass_dust
- self.disk_mass = np.trapz(2.0 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask],
- self.r[enclosed_mask])
- self.disk_mass_dust = np.trapz(2 * np.pi * self.sigma_dust[enclosed_mask] * self.r[enclosed_mask],
- self.r[enclosed_mask])
+ self.disk_mass = np.trapz(
+ 2.0 * np.pi * self.sigma_g[enclosed_mask] * self.r[enclosed_mask],
+ self.r[enclosed_mask],
+ )
+ self.disk_mass_dust = np.trapz(
+ 2 * np.pi * self.sigma_dust[enclosed_mask] * self.r[enclosed_mask],
+ self.r[enclosed_mask],
+ )
- #assert (self.disk_mass + self.disk_mass_dust) <= (
+ # assert (self.disk_mass + self.disk_mass_dust) <= (
# self._old_disk_mass + self._old_disk_mass_dust), "Disk mass is growing"
self.m_dot_components = self.compute_mdot_at_interfaces(self.sigma_g_components)
self.m_dot = np.sum(self.m_dot_components[:, 0, :], axis=1)
self.P = (
- self.c_s ** 2.0
- * self.sigma_g
- / (np.sqrt(2.0 * np.pi) * self.aspect_ratio * self.r)
+ self.c_s**2.0
+ * self.sigma_g
+ / (np.sqrt(2.0 * np.pi) * self.aspect_ratio * self.r)
)
- self.grad_sig = np.gradient(
- np.log10(self.sigma_g), np.log10(self.r)
- )
- self.grad_T = np.gradient(
- np.log10(self.T), np.log10(self.r)
- ) # -0.5 for MMSN
+ self.grad_sig = np.gradient(np.log10(self.sigma_g), np.log10(self.r))
+ self.grad_T = np.gradient(np.log10(self.T), np.log10(self.r)) # -0.5 for MMSN
self.grad_p = np.gradient(np.log10(self.P), np.log10(self.r))
- self.eta = -0.5 * self.aspect_ratio ** 2 * self.grad_p
+ self.eta = -0.5 * self.aspect_ratio**2 * self.grad_p
self.solid_fraction = self.sigma_dust / self.sigma_g
@@ -571,7 +738,9 @@ def get_viscous_evolution_next_timestep_components(self, sigma_g_components, dt)
d = 3.0 * nu[:, np.newaxis]
v = np.zeros(y.shape)
- s = 1 * self.sigdot # if self.t < self.begin_photevap else np.zeros_like(self.sigdot)
+ s = (
+ 1 * self.sigdot
+ ) # if self.t < self.begin_photevap else np.zeros_like(self.sigdot)
if hasattr(self, "sigdot_gas_acc"):
s -= self.sigdot_gas_acc
@@ -581,7 +750,15 @@ def get_viscous_evolution_next_timestep_components(self, sigma_g_components, dt)
s -= self.sigdot_photoevap
# checke das vorzeichen!
- clip_val = (y - 2.0 * np.pi * 0.9 * x[:, np.newaxis] * self.sigmin * self.chemistry_gas[:, 1, :]) / dt
+ clip_val = (
+ y
+ - 2.0
+ * np.pi
+ * 0.9
+ * x[:, np.newaxis]
+ * self.sigmin
+ * self.chemistry_gas[:, 1, :]
+ ) / dt
s = np.clip(s, -clip_val, clip_val)
#
# Set boundary conditions
@@ -612,8 +789,12 @@ def calc_pebble_flux(self, noupdate=False):
-------
"""
- vr_peb = interp1d(self.r_i[1:-1], self.vr_peb, assume_sorted=True, fill_value="extrapolate")(self.r)
- self.pebble_flux = np.abs(2.0 * np.pi * self.r * vr_peb * self.sigma_dust * self.f_m)
+ vr_peb = interp1d(
+ self.r_i[1:-1], self.vr_peb, assume_sorted=True, fill_value="extrapolate"
+ )(self.r)
+ self.pebble_flux = np.abs(
+ 2.0 * np.pi * self.r * vr_peb * self.sigma_dust * self.f_m
+ )
if not noupdate:
self.cum_pebble_flux += self.pebble_flux * self.dt
@@ -637,32 +818,56 @@ def calc_sigdot(self, icelines_mol):
y = 2.0 * np.pi * x[:, np.newaxis] * self.sigma_g_components[:, 1, :]
self.sigdot = np.zeros(y.shape)
# interpolate back to cell centers
- vr_dust = interp1d(self.r_i[1:-1], self.vr_dust, assume_sorted=True, fill_value="extrapolate")(self.r)
+ vr_dust = interp1d(
+ self.r_i[1:-1], self.vr_dust, assume_sorted=True, fill_value="extrapolate"
+ )(self.r)
if not self.conf_evaporation:
included_icelines = []
else:
# curate icelines: (only use icelines that are well enough inside of the grid)
- #included_icelines = np.arange(icelines_mol.shape[0])[
+ # included_icelines = np.arange(icelines_mol.shape[0])[
# np.logical_and(self.r[icelines_mol] > self.r[5], self.r[icelines_mol] < self.r[-5])]
included_icelines = np.arange(1, icelines_mol.shape[0]) # exclude rest_mol
- sigdot_cond_const = -3.0 * self.epsilon_p / (
- 2.0 * np.pi * self.rho_solid[:, np.newaxis]) * y * self.sigma_dust[:, np.newaxis] * self.omega_k[:,
- np.newaxis] * np.sqrt(
- self._chem.mu[:,np.newaxis])
+ sigdot_cond_const = (
+ -3.0
+ * self.epsilon_p
+ / (2.0 * np.pi * self.rho_solid[:, np.newaxis])
+ * y
+ * self.sigma_dust[:, np.newaxis]
+ * self.omega_k[:, np.newaxis]
+ * np.sqrt(self._chem.mu[:, np.newaxis])
+ )
sigdot_cond_const *= (
- self.f_m[:, np.newaxis] / self.a_1[:, np.newaxis] + (1.0 - self.f_m[:, np.newaxis]) / self.a_0)
+ self.f_m[:, np.newaxis] / self.a_1[:, np.newaxis]
+ + (1.0 - self.f_m[:, np.newaxis]) / self.a_0
+ )
for i in included_icelines: # exclude rest_mol
- self.sigdot[icelines_mol[i]:-1, i] = sigdot_cond_const[icelines_mol[i]:-1, i] / np.sqrt(np.squeeze(self._chem.M_array)[i])
- self.sigdot[:icelines_mol[i], i] = 2 * np.pi * x[:icelines_mol[i]] * self.sigma_dust_components[
- :icelines_mol[i], 1, i] * np.abs(
- vr_dust[:icelines_mol[i]]) / self._evap_width
+ self.sigdot[icelines_mol[i] : -1, i] = sigdot_cond_const[
+ icelines_mol[i] : -1, i
+ ] / np.sqrt(np.squeeze(self._chem.M_array)[i])
+ self.sigdot[: icelines_mol[i], i] = (
+ 2
+ * np.pi
+ * x[: icelines_mol[i]]
+ * self.sigma_dust_components[: icelines_mol[i], 1, i]
+ * np.abs(vr_dust[: icelines_mol[i]])
+ / self._evap_width
+ )
- min_density = 2.0 * np.pi * x[:, np.newaxis] * 0.9 * np.minimum(self.sigma_dust_components[:, 1, :],
- self.sigma_g_components[:, 1, :])
+ min_density = (
+ 2.0
+ * np.pi
+ * x[:, np.newaxis]
+ * 0.9
+ * np.minimum(
+ self.sigma_dust_components[:, 1, :], self.sigma_g_components[:, 1, :]
+ )
+ )
- self.sigdot = np.clip(self.sigdot, -min_density / self.dt,
- min_density / self.dt) # do not transfer more then available in one timestep
+ self.sigdot = np.clip(
+ self.sigdot, -min_density / self.dt, min_density / self.dt
+ ) # do not transfer more then available in one timestep
return
def calc_sigdot_photoevap(self):
@@ -675,10 +880,18 @@ def calc_sigdot_photoevap(self):
"""
if self.t > self.begin_photevap:
if hasattr(self, "tau_disk"):
- self.sigdot_photoevap = 1.0 / self.tau_disk * self.sigma_g_components[:, 1, :] * 2.0 * np.pi * self.r[:,
- np.newaxis]
+ self.sigdot_photoevap = (
+ 1.0
+ / self.tau_disk
+ * self.sigma_g_components[:, 1, :]
+ * 2.0
+ * np.pi
+ * self.r[:, np.newaxis]
+ )
else:
- self.sigdot_photoevap = np.zeros_like(self.r) * self.chemistry_gas[:, 1, :]
+ self.sigdot_photoevap = (
+ np.zeros_like(self.r) * self.chemistry_gas[:, 1, :]
+ )
def get_dust_radial_drift_next_timestep(self, dt):
"""
@@ -706,20 +919,31 @@ def get_dust_radial_drift_next_timestep(self, dt):
# lmax = np.ones_like(self.r, dtype="bool")
lmax = self.T < 2000.0
x = self.r[lmax]
- y = 2.0 * np.pi * x[:, np.newaxis] * self.sigma_dust_components[lmax, 1, 1:] # Dust
+ y = (
+ 2.0 * np.pi * x[:, np.newaxis] * self.sigma_dust_components[lmax, 1, 1:]
+ ) # Dust
g = 2.0 * np.pi * x * self.sigma_g[lmax] # Gas
g = g[:, np.newaxis]
d = self.nu[lmax]
di = 0.5 * (d[1:] + d[:-1])[:, np.newaxis]
vi = self.vr_dust[lmax[:-1], np.newaxis]
- s = - self.sigdot[lmax, 1:] # if self.t < self.begin_photevap else np.zeros_like(self.sigdot[:, 1:])
+ s = -self.sigdot[
+ lmax, 1:
+ ] # if self.t < self.begin_photevap else np.zeros_like(self.sigdot[:, 1:])
if hasattr(self, "sigdot_peb_acc"):
s -= self.sigdot_peb_acc[lmax, 1:]
- with np.errstate(under='ignore'):
- clip_val = (y - 2 * np.pi * x[:, np.newaxis] * self.sigmin_peb * self.chemistry_solid[lmax, 1, 1:]) / dt
+ with np.errstate(under="ignore"):
+ clip_val = (
+ y
+ - 2
+ * np.pi
+ * x[:, np.newaxis]
+ * self.sigmin_peb
+ * self.chemistry_solid[lmax, 1, 1:]
+ ) / dt
s = np.clip(s, -clip_val, clip_val)
#
@@ -728,21 +952,39 @@ def get_dust_radial_drift_next_timestep(self, dt):
bcl = (1, 0, 0, 1)
# bcr = (1, 0, 0, 1)
# bcl = (0, 1, 2 * np.pi * x[0] * self.chemistry_solid[0, 1, 1:] * self.sigmin_peb, 0)
- bcr = (0, 1, 2 * np.pi * x[-1] * self.chemistry_solid[-1, 1, 1:] * self.sigmin_peb, 0)
+ bcr = (
+ 0,
+ 1,
+ 2 * np.pi * x[-1] * self.chemistry_solid[-1, 1, 1:] * self.sigmin_peb,
+ 0,
+ )
#
# Get the new value of y after one time step dt
#
- y = solvediffonedee_components(x, y, vi, di, g, s, bcl, bcr, dt=dt, int=True, upwind=True)
- y = np.maximum(y, 2 * np.pi * x[:, np.newaxis] * self.chemistry_solid[lmax, 1, 1:] * self.sigmin_peb)
+ y = solvediffonedee_components(
+ x, y, vi, di, g, s, bcl, bcr, dt=dt, int=True, upwind=True
+ )
+ y = np.maximum(
+ y,
+ 2
+ * np.pi
+ * x[:, np.newaxis]
+ * self.chemistry_solid[lmax, 1, 1:]
+ * self.sigmin_peb,
+ )
# y = np.maximum(y, 2 * np.pi * x[:, np.newaxis] * 1e-60)
#
# Obtain new sigdust
#
- new_sigma_dust = y / (2 * np.pi * x[:, np.newaxis]) # add offset for numerical reasons
+ new_sigma_dust = y / (
+ 2 * np.pi * x[:, np.newaxis]
+ ) # add offset for numerical reasons
#
# Return
#
- sigma_dust = np.ones_like(self.sigma_dust_components[:, 1, 1:]) * self.sigmin_peb
+ sigma_dust = (
+ np.ones_like(self.sigma_dust_components[:, 1, 1:]) * self.sigmin_peb
+ )
sigma_dust[lmax] = new_sigma_dust
return sigma_dust
@@ -770,10 +1012,15 @@ def compute_sigma_g_from_mol(self, sigma_g_components_mol):
def compute_sigma_dust_from_mol(self, sigma_dust_components_mol):
if self.conf_evaporation:
- self.chemistry_solid = self._chem.get_solid_composition(sigma_dust_components_mol, self.T)
+ self.chemistry_solid = self._chem.get_solid_composition(
+ sigma_dust_components_mol, self.T
+ )
- self.sigma_dust = np.where(np.sum(self.chemistry_solid[:, 0], axis=1) > 0.0,
- np.sum(sigma_dust_components_mol, axis=1), np.ones_like(self.r) * self.sigmin_peb)
+ self.sigma_dust = np.where(
+ np.sum(self.chemistry_solid[:, 0], axis=1) > 0.0,
+ np.sum(sigma_dust_components_mol, axis=1),
+ np.ones_like(self.r) * self.sigmin_peb,
+ )
assert np.all(self.sigma_dust > 0), "negative surface densities!"
def compute_mdot_at_interfaces(self, sigma_g, oned=False):
@@ -788,14 +1035,16 @@ def compute_mdot_at_interfaces(self, sigma_g, oned=False):
y = 2.0 * np.pi * self.r[:, np.newaxis, np.newaxis] * sigma_g
else:
y = 2.0 * np.pi * self.r * sigma_g
- g = self.r ** 0.5 / (self.nu / self.alpha_factor)
+ g = self.r**0.5 / (self.nu / self.alpha_factor)
d = 3.0 * self.nu / self.alpha_factor
v = np.zeros(len(x))
#
# Compute the mdot
#
- flux = -getfluxonedee(x, y, v, d, g, int=False, oned=oned) # returns flux at interfaces
+ flux = -getfluxonedee(
+ x, y, v, d, g, int=False, oned=oned
+ ) # returns flux at interfaces
return np.array([flux[0], *flux, flux[-1]])
def compute_vr_at_interfaces(self):
@@ -816,7 +1065,7 @@ def compute_vr_at_interfaces(self):
#
# Now compute the radial velocity from the mdot
#
- self.vr_gas = - self.m_dot[1:-1] / (q + 1e-90)
+ self.vr_gas = -self.m_dot[1:-1] / (q + 1e-90)
def calc_velocity_from_stokes(self, stokes, dlnpdlnr, csi, vki, vrgas):
"""
@@ -834,8 +1083,9 @@ def calc_velocity_from_stokes(self, stokes, dlnpdlnr, csi, vki, vrgas):
"""
stokes_i = (stokes[1:] + stokes[:-1]) / 2.0 + 1e-60
- return vrgas / (1.0 + stokes_i ** 2) + dlnpdlnr * csi ** 2 / vki / (
- stokes_i + 1.0 / stokes_i)
+ return vrgas / (1.0 + stokes_i**2) + dlnpdlnr * csi**2 / vki / (
+ stokes_i + 1.0 / stokes_i
+ )
def compute_interface_velocities(self):
"""
@@ -856,9 +1106,9 @@ def compute_interface_velocities(self):
#
self.P = (
- self.c_s ** 2.0
- * self.sigma_g
- / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r)
+ self.c_s**2.0
+ * self.sigma_g
+ / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r)
)
self.grad_p = np.gradient(
np.log10(self.P), np.log10(self.r)
@@ -871,22 +1121,28 @@ def compute_interface_velocities(self):
#
if not self.conf_static_stokes:
f_m_i = (self.f_m[1:] + self.f_m[:-1]) / 2.0
- self.vr_peb = self.calc_velocity_from_stokes(self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas)
- self.vr_small = self.calc_velocity_from_stokes(self.stokes_number_small, dlnpdlnr, csi, vki, vrgas)
+ self.vr_peb = self.calc_velocity_from_stokes(
+ self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas
+ )
+ self.vr_small = self.calc_velocity_from_stokes(
+ self.stokes_number_small, dlnpdlnr, csi, vki, vrgas
+ )
self.vr_dust = f_m_i * self.vr_peb + (1.0 - f_m_i) * self.vr_small
else:
- self.vr_dust = self.calc_velocity_from_stokes(self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas)
+ self.vr_dust = self.calc_velocity_from_stokes(
+ self.stokes_number_pebbles, dlnpdlnr, csi, vki, vrgas
+ )
self.vr_small = 0.0
self.vr_peb = 1.0 * self.vr_dust
def init_T(self):
- flux = self.lstar / (4.0 * np.pi * self.r ** 2.0)
+ flux = self.lstar / (4.0 * np.pi * self.r**2.0)
self.qirr = (
- flux * self.flang
+ flux * self.flang
) # Factor 2 for two sides, factor 0.5 for half irradiated down
self.T_irr = (
- 0.5 * self.qirr / sr
- ) ** 0.25 # Factor 0.5 because cooling is two-sided
+ 0.5 * self.qirr / sr
+ ) ** 0.25 # Factor 0.5 because cooling is two-sided
self.T_irr = np.maximum(self.T_irr, 10.0)
self.T = 1.0 * self.T_irr
@@ -922,13 +1178,13 @@ def compute_disktmid(self):
self.qirr = Irradiative heating rate of the disk [erg/cm^2/s]
self.qvisc = (if vischeat==True) Viscous heating rate [erg/cm^2/s]
"""
- if hasattr(self,"alpha_factor"):
+ if hasattr(self, "alpha_factor"):
alpha = self.alpha / self.alpha_factor
else:
alpha = self.alpha
if hasattr(self, "T_irr"):
- if (self.conf_temp_evol or not self._T_init_finished):
+ if self.conf_temp_evol or not self._T_init_finished:
self.T = solve_viscous_heating_globally(
sigma_g=self.sigma_g,
tirr=self.T_irr,
@@ -938,9 +1194,12 @@ def compute_disktmid(self):
Z=self.DTG_small_grains,
r=self.r,
mu=self._chem.mu,
- interpolation_idx=None if not hasattr(self,"interpolation_idx") else self.interpolation_idx)
+ interpolation_idx=None
+ if not hasattr(self, "interpolation_idx")
+ else self.interpolation_idx,
+ )
- self.T_visc = (self.T ** 4.0 - self.T_irr ** 4.0) ** 0.25
+ self.T_visc = (self.T**4.0 - self.T_irr**4.0) ** 0.25
else:
self.init_T()
@@ -957,14 +1216,14 @@ def calc_opacity(self):
self.opacity = opacity_fct(self.rho_g, self.T, Z=self.DTG_small_grains)
def compute_cs_and_hp(self):
- """ update all quantities related to T """
+ """update all quantities related to T"""
self.c_s = (k_B * self.T / (self._chem.mu * m_p)) ** 0.5
- self.aspect_ratio = self.c_s / (
- self.omega_k * self.r
- )
+ self.aspect_ratio = self.c_s / (self.omega_k * self.r)
if hasattr(self, "sigma_g"):
- self.rho_g = self.sigma_g / (np.sqrt(2 * np.pi) * self.aspect_ratio * self.r)
+ self.rho_g = self.sigma_g / (
+ np.sqrt(2 * np.pi) * self.aspect_ratio * self.r
+ )
def compute_nu(self):
- """ compute the viscosity. """
+ """compute the viscosity."""
self.nu = self.alpha * self.c_s * self.c_s / self.omega_k
diff --git a/chemcomp/disks/kees.py b/chemcomp/disks/kees.py
index 5d27ebc..04868b4 100755
--- a/chemcomp/disks/kees.py
+++ b/chemcomp/disks/kees.py
@@ -11,18 +11,18 @@
class DiskLabDisk(Disk):
- """ LBP, viscous disk evolution and viscous heating. Disk of Schneider & Bitsch (2021)a,b."""
+ """LBP, viscous disk evolution and viscous heating. Disk of Schneider & Bitsch (2021)a,b."""
def __init__(
- self,
- defaults,
- chemistry,
- M_STAR=1.0 * u.M_sun,
- M0=1e-2 * u.M_sun,
- R0=200 * u.au,
- time=1.0 * u.yr,
- ALPHA=1e-3,
- **kwargs
+ self,
+ defaults,
+ chemistry,
+ M_STAR=1.0 * u.M_sun,
+ M0=1e-2 * u.M_sun,
+ R0=200 * u.au,
+ time=1.0 * u.yr,
+ ALPHA=1e-3,
+ **kwargs
):
super().__init__(defaults, M_STAR, ALPHA, chemistry, time=time, **kwargs)
@@ -35,12 +35,17 @@ def __init__(
self.init_T()
t_iter = 0
while True:
- self.make_disk_from_lbp_alpha(M0=M0.cgs.value, R0=R0.cgs.value, alpha=float(ALPHA), time=time.cgs.value)
+ self.make_disk_from_lbp_alpha(
+ M0=M0.cgs.value,
+ R0=R0.cgs.value,
+ alpha=float(ALPHA),
+ time=time.cgs.value,
+ )
T0 = self.T.copy()
self.compute_disktmid()
self._chem.get_composition(self.T) # updates mu
self._init_params(**kwargs)
- if np.all(np.abs((T0-self.T)/T0) < 1e-2):
+ if np.all(np.abs((T0 - self.T) / T0) < 1e-2):
break
assert t_iter < 1000, "We have trouble to converge the initial T-profile"
@@ -63,33 +68,37 @@ def make_disk_from_simplified_lbp(self, Sigc, Rc, gam):
gam = The gamma coefficient of the LBP model
"""
r = self.r / AU
- sigma = Sigc * (Rc / r) ** gam * np.exp(-(r / Rc) ** (2. - gam))
+ sigma = Sigc * (Rc / r) ** gam * np.exp(-((r / Rc) ** (2.0 - gam)))
sigma[sigma < self.sigmin] = self.sigmin
self.sigma_g = sigma
def make_disk_from_lyndenbellpringle(self, M0, R0, nu0, gam, time):
"""
- Make a model of the gas surface density of the disk.
- Here the model is set by the Lynden-Bell & Pringle solution.
- I follow here Lodato, Scardoni, Manara & Testi (2017)
- MNRAS 472, 4700, their Eqs. 2, 3 and 4.
-
- ARGUMENTS:
- M0 = Initial disk mass [g]
- R0 = Initial disk radius [cm]
- nu0 = Viscosity nu at R0 [cm^2/s]
- gam = The gamma coefficient of the LBP model
- time = Time after start [s]
- """
+ Make a model of the gas surface density of the disk.
+ Here the model is set by the Lynden-Bell & Pringle solution.
+ I follow here Lodato, Scardoni, Manara & Testi (2017)
+ MNRAS 472, 4700, their Eqs. 2, 3 and 4.
+
+ ARGUMENTS:
+ M0 = Initial disk mass [g]
+ R0 = Initial disk radius [cm]
+ nu0 = Viscosity nu at R0 [cm^2/s]
+ gam = The gamma coefficient of the LBP model
+ time = Time after start [s]
+ """
r = self.r
nu = nu0 * (r / R0) ** gam # Eq. 2
- tnu = R0 ** 2 / (3.0 * (2. - gam) ** 2 * nu0) # Eq. 4
+ tnu = R0**2 / (3.0 * (2.0 - gam) ** 2 * nu0) # Eq. 4
T = 1.0 + time / tnu # Above Eq. 4
- eta = (2.5 - gam) / (2. - gam) # Below Eq. 3
- self.sigma_g = (M0 / (2 * np.pi * R0 ** 2.0)) * (2. - gam) * \
- (R0 / r) ** gam * T ** (-eta) * \
- np.exp(-(r / R0) ** (2. - gam) / T) # Eq. 3
+ eta = (2.5 - gam) / (2.0 - gam) # Below Eq. 3
+ self.sigma_g = (
+ (M0 / (2 * np.pi * R0**2.0))
+ * (2.0 - gam)
+ * (R0 / r) ** gam
+ * T ** (-eta)
+ * np.exp(-((r / R0) ** (2.0 - gam)) / T)
+ ) # Eq. 3
self.nu = nu
self.gam = gam
self.sigma_g[self.sigma_g < self.sigmin] = self.sigmin
@@ -110,7 +119,7 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time):
"""
assert hasattr(self, "aspect_ratio")
- self.omega_k = np.sqrt(G * self.M_s / self.r ** 3.0)
+ self.omega_k = np.sqrt(G * self.M_s / self.r**3.0)
self.alpha = alpha
r = self.r
aspect_ratio = 1 * self.aspect_ratio
@@ -123,5 +132,5 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time):
self.make_disk_from_lyndenbellpringle(M0, R0, nu0, gam, time)
def update_parameters(self):
- """ Things that should be called whenever the update_parameters() function is called in the planet class. """
+ """Things that should be called whenever the update_parameters() function is called in the planet class."""
self.compute_viscous_evolution()
diff --git a/chemcomp/disks/mmsn.py b/chemcomp/disks/mmsn.py
index 3cc91b7..e8ec6f8 100755
--- a/chemcomp/disks/mmsn.py
+++ b/chemcomp/disks/mmsn.py
@@ -18,27 +18,33 @@
class MMSN(Disk):
- """ docstring for MMSN. Hayashi 1981, Weidenschilling 1977. Likely outdated"""
-
- def __init__(self,
- defaults,
- chemistry,
- M_STAR=1 * u.M_sun,
- R0=80.0 * u.au,
- ALPHA=1e-3,
- SIGMA_POWER=-15 / 14,
- ASPECT_POWER=2 / 7,
- SIGMA_0=1500 * u.g / (u.cm ** 2),
- ASPECT_0=0.033,
- **kwargs):
+ """docstring for MMSN. Hayashi 1981, Weidenschilling 1977. Likely outdated"""
+
+ def __init__(
+ self,
+ defaults,
+ chemistry,
+ M_STAR=1 * u.M_sun,
+ R0=80.0 * u.au,
+ ALPHA=1e-3,
+ SIGMA_POWER=-15 / 14,
+ ASPECT_POWER=2 / 7,
+ SIGMA_0=1500 * u.g / (u.cm**2),
+ ASPECT_0=0.033,
+ **kwargs
+ ):
super().__init__(defaults, M_STAR, ALPHA, chemistry, **kwargs)
- self.sigma_init = SIGMA_0.cgs.value * (self.r / AU) ** SIGMA_POWER * np.exp(- self.r / R0.cgs.value)
+ self.sigma_init = (
+ SIGMA_0.cgs.value
+ * (self.r / AU) ** SIGMA_POWER
+ * np.exp(-self.r / R0.cgs.value)
+ )
self.aspect_ratio = ASPECT_0 * (self.r / AU) ** ASPECT_POWER # Hayashi 1981
self._init_params(**kwargs)
self.calc_opacity()
def compute_disktmid(self, viscmodel=None):
- """calculate the disk temperature. """
+ """calculate the disk temperature."""
self.compute_cs_and_hp()
self.compute_nu()
self.calc_opacity()
@@ -49,14 +55,21 @@ def calc_sigma_dust_static(self):
This is equal to the pebble surface density, if DTG: peb has been set. Otherwise it contains the whole dust (including also the small grains)
"""
- vr = 2.0 * self.stokes_number_pebbles / (self.stokes_number_pebbles ** 2. + 1.0) * self.eta * self.v_k
+ vr = (
+ 2.0
+ * self.stokes_number_pebbles
+ / (self.stokes_number_pebbles**2.0 + 1.0)
+ * self.eta
+ * self.v_k
+ )
Pebbleflux = 10e-4 * ME / year
- self.sigma_dust = Pebbleflux / (2. * np.pi * self.r * vr)
+ self.sigma_dust = Pebbleflux / (2.0 * np.pi * self.r * vr)
self.chemistry_solid, _ = self._chem.get_composition(self.T)
- self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis,
- np.newaxis]
+ self.sigma_dust_components = (
+ self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis]
+ )
self.sigma_dust_components[:, 1, 0] = np.zeros_like(self.sigma_g)
@@ -64,7 +77,7 @@ def calc_sigma_dust_static(self):
return
def update_parameters(self):
- """ compute the viscous evolution. This function is called whenever the planet calls the update."""
+ """compute the viscous evolution. This function is called whenever the planet calls the update."""
self.compute_viscous_evolution()
def calc_opacity_oplin(self):
diff --git a/chemcomp/disks/twopop.py b/chemcomp/disks/twopop.py
index 959db67..be366a1 100755
--- a/chemcomp/disks/twopop.py
+++ b/chemcomp/disks/twopop.py
@@ -11,20 +11,21 @@
G = const.G.cgs.value
+
class TwoPopDisk(Disk):
- """ LBP, Disk of the twopoppy model."""
+ """LBP, Disk of the twopoppy model."""
def __init__(
- self,
- defaults,
- chemistry,
- M_STAR=0.7 * u.M_sun,
- M0=1e-1 * u.M_sun,
- R0=20 * u.au,
- time=1 * u.yr,
- ALPHA=1e-3,
- output={},
- **kwargs
+ self,
+ defaults,
+ chemistry,
+ M_STAR=0.7 * u.M_sun,
+ M0=1e-1 * u.M_sun,
+ R0=20 * u.au,
+ time=1 * u.yr,
+ ALPHA=1e-3,
+ output={},
+ **kwargs
):
super().__init__(defaults, M_STAR, ALPHA, chemistry, time=time, **kwargs)
@@ -34,9 +35,15 @@ def __init__(
self.tstar = self.tstar.cgs.value
self.compute_disktmid()
- self.make_disk_from_lbp_alpha(M0=M0.cgs.value, R0=R0.cgs.value, alpha=ALPHA, time=time.cgs.value)
+ self.make_disk_from_lbp_alpha(
+ M0=M0.cgs.value, R0=R0.cgs.value, alpha=ALPHA, time=time.cgs.value
+ )
- self.sigma_g = self.sigma_g / np.trapz(2 * np.pi * self.r * self.sigma_g, x=self.r) * M0.cgs.value
+ self.sigma_g = (
+ self.sigma_g
+ / np.trapz(2 * np.pi * self.r * self.sigma_g, x=self.r)
+ * M0.cgs.value
+ )
self.compute_disktmid()
self.compute_disktmid()
@@ -46,16 +53,19 @@ def __init__(
def evolve_init(self):
pass
# DEBUGGING
- #mask = np.where(self.r < (5 * u.au).cgs.value)[0]
- #self.sigma_g_components[mask, 1, 7] = 1e-60
- #self.sigma_dust_components[mask, :, :] = 1e-60
- #self.compute_sigma_g_from_mol(self.sigma_g_components[:, 1, :])
- #self.compute_sigma_dust_from_mol(self.sigma_dust_components[:, 1, :])
- #self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis]
- #self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
+ # mask = np.where(self.r < (5 * u.au).cgs.value)[0]
+ # self.sigma_g_components[mask, 1, 7] = 1e-60
+ # self.sigma_dust_components[mask, :, :] = 1e-60
+ # self.compute_sigma_g_from_mol(self.sigma_g_components[:, 1, :])
+ # self.compute_sigma_dust_from_mol(self.sigma_dust_components[:, 1, :])
+ # self.sigma_dust_components = self.chemistry_solid * self.sigma_dust[:, np.newaxis, np.newaxis]
+ # self.sigma_g_components = self.chemistry_gas * self.sigma_g[:, np.newaxis, np.newaxis]
def compute_disktmid(self, viscmodel=None):
- self.T = ((0.05 ** 0.25 * self.tstar * (self.r / self.rstar) ** -0.5) ** 4 + (7.) ** 4.0) ** 0.25
+ self.T = (
+ (0.05**0.25 * self.tstar * (self.r / self.rstar) ** -0.5) ** 4
+ + (7.0) ** 4.0
+ ) ** 0.25
self.tvisc = np.zeros_like(self.T)
self.tirr = np.zeros_like(self.T)
self.compute_cs_and_hp()
@@ -64,31 +74,34 @@ def compute_disktmid(self, viscmodel=None):
def make_disk_from_lyndenbellpringle(self, M0, R0, nu0, gam, time):
"""
- Make a model of the gas surface density of the disk.
- Here the model is set by the Lynden-Bell & Pringle solution.
- I follow here Lodato, Scardoni, Manara & Testi (2017)
- MNRAS 472, 4700, their Eqs. 2, 3 and 4.
-
- ARGUMENTS:
- M0 = Initial disk mass [g]
- R0 = Initial disk radius [cm]
- nu0 = Viscosity nu at R0 [cm^2/s]
- gam = The gamma coefficient of the LBP model
- time = Time after start [s]
- """
+ Make a model of the gas surface density of the disk.
+ Here the model is set by the Lynden-Bell & Pringle solution.
+ I follow here Lodato, Scardoni, Manara & Testi (2017)
+ MNRAS 472, 4700, their Eqs. 2, 3 and 4.
+
+ ARGUMENTS:
+ M0 = Initial disk mass [g]
+ R0 = Initial disk radius [cm]
+ nu0 = Viscosity nu at R0 [cm^2/s]
+ gam = The gamma coefficient of the LBP model
+ time = Time after start [s]
+ """
r = self.r
nu = nu0 * (r / R0) ** gam # Eq. 2
- tnu = R0 ** 2.0 / (3.0 * (2. - gam) ** 2.0 * nu0) # Eq. 4
+ tnu = R0**2.0 / (3.0 * (2.0 - gam) ** 2.0 * nu0) # Eq. 4
T = 1.0 + time / tnu # Above Eq. 4
- eta = (2.5 - gam) / (2. - gam) # Below Eq. 3
- self.sigma_g = (M0 / (2.0 * np.pi * R0 ** 2)) * (2. - gam) * \
- (R0 / r) ** gam * T ** (-eta) * \
- np.exp(-(r / R0) ** (2. - gam) / T) # Eq. 3
+ eta = (2.5 - gam) / (2.0 - gam) # Below Eq. 3
+ self.sigma_g = (
+ (M0 / (2.0 * np.pi * R0**2))
+ * (2.0 - gam)
+ * (R0 / r) ** gam
+ * T ** (-eta)
+ * np.exp(-((r / R0) ** (2.0 - gam)) / T)
+ ) # Eq. 3
self.nu = nu
self.gam = gam
self.sigma_g[self.sigma_g < self.sigmin] = self.sigmin
-
def make_disk_from_lbp_alpha(self, M0, R0, alpha, time):
"""
Same as make_disk_from_lyndenbellpringle(), but now assuming a constant
@@ -105,7 +118,7 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time):
"""
assert hasattr(self, "aspect_ratio")
- self.omega_k = np.sqrt(G * self.M_s / self.r ** 3.0)
+ self.omega_k = np.sqrt(G * self.M_s / self.r**3.0)
self.alpha = alpha
r = self.r
aspect_ratio = 1.0 * self.aspect_ratio
@@ -114,10 +127,9 @@ def make_disk_from_lbp_alpha(self, M0, R0, alpha, time):
gam = 1.0
om1 = self.omega_k[0]
cs1 = cs[0]
- nu0 = self.alpha * cs1 ** 2.0 / om1
+ nu0 = self.alpha * cs1**2.0 / om1
self.make_disk_from_lyndenbellpringle(M0, R0, nu0, gam, time)
def update_parameters(self):
self.compute_viscous_evolution()
-
diff --git a/chemcomp/helpers/__init__.py b/chemcomp/helpers/__init__.py
index be4bf33..8d4fb3f 100755
--- a/chemcomp/helpers/__init__.py
+++ b/chemcomp/helpers/__init__.py
@@ -1,10 +1,17 @@
import os
from ast import literal_eval
-CONFIG_PATH = os.path.join(os.getcwd(), "config") # set default path of config directory
+CONFIG_PATH = os.path.join(
+ os.getcwd(), "config"
+) # set default path of config directory
JOB_PATH = os.path.join(os.getcwd(), "jobs") # set default path of config directory
-OUTPUT_PATH = os.path.join(os.getcwd(), "output") # set default path of output directory
-ZIP_OUTPUT_PATH = os.path.join(os.getcwd(), "zip_output") # set default path of output directory
+OUTPUT_PATH = os.path.join(
+ os.getcwd(), "output"
+) # set default path of output directory
+ZIP_OUTPUT_PATH = os.path.join(
+ os.getcwd(), "zip_output"
+) # set default path of output directory
+
def eval_kwargs(inp):
try:
@@ -20,4 +27,3 @@ def eval_kwargs(inp):
from .main_helpers import run_setup, run_pipeline
from .solvediffonedee import solvediffonedee_components, getfluxonedee
from .analysis_helper import GrowthPlot, molecules, elements
-
diff --git a/chemcomp/helpers/analysis_helper.py b/chemcomp/helpers/analysis_helper.py
index b6ba31f..8e60629 100755
--- a/chemcomp/helpers/analysis_helper.py
+++ b/chemcomp/helpers/analysis_helper.py
@@ -12,9 +12,27 @@
np.seterr(divide="ignore", invalid="ignore")
elements = ["C", "O", "Fe", "S", "Mg", "Si", "Na", "K", "N", "Al", "Ti", "V", "H", "He"]
-molecules = ["rest", "CO", "N2", "CH4", "CO2", "NH3", "H2S", "H2O", "Fe3O4", "C", "FeS", "NaAlSi3O8", "KAlSi3O8",
- "Mg2SiO4",
- "Fe2O3", "VO", "MgSiO3", "Al2O3", "TiO"]
+molecules = [
+ "rest",
+ "CO",
+ "N2",
+ "CH4",
+ "CO2",
+ "NH3",
+ "H2S",
+ "H2O",
+ "Fe3O4",
+ "C",
+ "FeS",
+ "NaAlSi3O8",
+ "KAlSi3O8",
+ "Mg2SiO4",
+ "Fe2O3",
+ "VO",
+ "MgSiO3",
+ "Al2O3",
+ "TiO",
+]
FeH_dict = {
"FeH": [-0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4],
@@ -54,8 +72,22 @@
MassC = 12.0
MassV = 50.9415
-el_masses = {"C": MassC, "O": MassO, "H": MassH, "Fe": MassFe, "Mg": MassMg, "Si": MassSi, "N": MassN, "S": MassS,
- "He": MassHe, "Al": MassAl, "Ti": MassTi, "K": MassK, "Na": MassNa, "V": MassV}
+el_masses = {
+ "C": MassC,
+ "O": MassO,
+ "H": MassH,
+ "Fe": MassFe,
+ "Mg": MassMg,
+ "Si": MassSi,
+ "N": MassN,
+ "S": MassS,
+ "He": MassHe,
+ "Al": MassAl,
+ "Ti": MassTi,
+ "K": MassK,
+ "Na": MassNa,
+ "V": MassV,
+}
DOT_SPACE = 1e5 * u.yr # time in years between dots
DOT_SPACE_LARGE = 5e5 * u.yr # time in years between dots
@@ -75,31 +107,63 @@
("densely dashdotdotted", (0, (3, 1, 1, 1, 1, 1))),
]
-dict_of_planet_read_names = {"t": r"t", "M": r"M", "M_a": r"M_a",
- "M_c": r"M_c", "a_p": r"a_p", "r_p": r"R_p", "tau_m": r"\tau_m", "T": "T",
- "comp_a": "Atmospheric Content", "comp_c": "Core Content", "sigma_g": r"\Sigma_g",
- "gamma_tot": r"\Gamma", "gamma_norm": r"\Gamma / \Gamma_0", "fSigma": "fSigma"}
-dict_of_acc_names_template = {r"m_dot": r"\dot M_{\mathrm{", r"m_dot_a": r"$\dot M_{\mathrm{atm,",
- r"m_dot_c": r"\dot M_{\mathrm{core,",
- r"m_dot_a_chem": r"\dot M_{\mathrm{atm,", r"m_dot_c_chem": r"\dot M_{\mathrm{core,",
- r"M_z": r"M_{Z,\mathrm{", r"regime":"regime_"}
-dict_of_acc_names = {f"{k}_{acc}": r"{}{}".format(v, acc) + r"}}" for k, v in dict_of_acc_names_template.items() for
- acc in ["peb", "gas"]}
-
-dict_of_extra_names = {"Z_vs_M": r"\frac{M_{Z}}{M}}",
- "Z_vs_M_a": r"\frac{M_{a, Z}}{M_a}}", "Z_vs_M_c": r"\frac{M_{c, Z}}{M_c}}"}
+dict_of_planet_read_names = {
+ "t": r"t",
+ "M": r"M",
+ "M_a": r"M_a",
+ "M_c": r"M_c",
+ "a_p": r"a_p",
+ "r_p": r"R_p",
+ "tau_m": r"\tau_m",
+ "T": "T",
+ "comp_a": "Atmospheric Content",
+ "comp_c": "Core Content",
+ "sigma_g": r"\Sigma_g",
+ "gamma_tot": r"\Gamma",
+ "gamma_norm": r"\Gamma / \Gamma_0",
+ "fSigma": "fSigma",
+}
+dict_of_acc_names_template = {
+ r"m_dot": r"\dot M_{\mathrm{",
+ r"m_dot_a": r"$\dot M_{\mathrm{atm,",
+ r"m_dot_c": r"\dot M_{\mathrm{core,",
+ r"m_dot_a_chem": r"\dot M_{\mathrm{atm,",
+ r"m_dot_c_chem": r"\dot M_{\mathrm{core,",
+ r"M_z": r"M_{Z,\mathrm{",
+ r"regime": "regime_",
+}
+dict_of_acc_names = {
+ f"{k}_{acc}": r"{}{}".format(v, acc) + r"}}"
+ for k, v in dict_of_acc_names_template.items()
+ for acc in ["peb", "gas"]
+}
+
+dict_of_extra_names = {
+ "Z_vs_M": r"\frac{M_{Z}}{M}}",
+ "Z_vs_M_a": r"\frac{M_{a, Z}}{M_a}}",
+ "Z_vs_M_c": r"\frac{M_{c, Z}}{M_c}}",
+}
NAME = {**dict_of_planet_read_names, **dict_of_acc_names, **dict_of_extra_names}
class GrowthPlot:
- def __init__(self, files, out, plotlabels=None, FeH=0, close_plots=False, quantities=None, dry_run=False):
+ def __init__(
+ self,
+ files,
+ out,
+ plotlabels=None,
+ FeH=0,
+ close_plots=False,
+ quantities=None,
+ dry_run=False,
+ ):
self.files = files
self.out = out
self._dry_run = dry_run
if quantities is not None:
self._used_quantities = set(quantities)
- self._used_quantities.update(["M","t"])
+ self._used_quantities.update(["M", "t"])
else:
self._used_quantities = quantities
self.import_data(files)
@@ -130,15 +194,27 @@ def solar_from_FeH(FeH):
NaHabu = 1 * NaHabu0
VHabu = 1 * VHabu0
- return {"O": OHabu * MassO, "Si": SiHabu * MassSi, "C": CHabu * MassC, "S": SHabu * MassS,
- "Fe": FeHabu * MassFe, "Mg": MgHabu * MassMg, "H": 1, "He": 0.085 * MassHe, "Na": NaHabu * MassNa,
- "N": NHabu * MassN, "Al": AlHabu * MassAl, "Ti": TiHabu * MassTi, "K": KHabu * MassK,
- "V": MassV * VHabu}
+ return {
+ "O": OHabu * MassO,
+ "Si": SiHabu * MassSi,
+ "C": CHabu * MassC,
+ "S": SHabu * MassS,
+ "Fe": FeHabu * MassFe,
+ "Mg": MgHabu * MassMg,
+ "H": 1,
+ "He": 0.085 * MassHe,
+ "Na": NaHabu * MassNa,
+ "N": NHabu * MassN,
+ "Al": AlHabu * MassAl,
+ "Ti": TiHabu * MassTi,
+ "K": KHabu * MassK,
+ "V": MassV * VHabu,
+ }
self.FeH = {}
for o in self.out:
if type(FeH) == dict:
- """ dict needs to have out as keys """
+ """dict needs to have out as keys"""
self.data[o]["solar"] = solar_from_FeH(FeH[o])
self.FeH[o] = FeH[o]
else:
@@ -154,8 +230,11 @@ def import_data(self, files):
with h5py.File(files[i], "r") as f:
units = dict(f["/planet/units"].attrs)
keys = list(f["/planet"].keys())
- keys = [k for k in keys if
- k in self._used_quantities] if self._used_quantities is not None else keys
+ keys = (
+ [k for k in keys if k in self._used_quantities]
+ if self._used_quantities is not None
+ else keys
+ )
if "units" in keys:
keys.remove("units")
break
@@ -175,18 +254,25 @@ def import_data(self, files):
data[k] = u.Quantity(d, units.get(k))
else:
data[k] = d
- data = {q: data[q][np.where(~np.isnan(data["M"]))] for q in data.keys()}
+ data = {
+ q: data[q][np.where(~np.isnan(data["M"]))]
+ for q in data.keys()
+ }
self.data[self.out[i]] = data
except OSError:
self.faulty_list.append(i)
if len(self.faulty_list) > 0:
- print(f"WARNING: There are in total {len(self.faulty_list)} corrupt files!")
+ print(
+ f"WARNING: There are in total {len(self.faulty_list)} corrupt files!"
+ )
else:
with h5py.File(files[0], "r") as f:
fill_data = {}
for k in keys:
- fill_data[k] = np.squeeze(u.Quantity(np.array(f["/planet"][k]), units.get(k)))
+ fill_data[k] = np.squeeze(
+ u.Quantity(np.array(f["/planet"][k]), units.get(k))
+ )
for i in range(len(self.out)):
data = {}
@@ -195,8 +281,10 @@ def import_data(self, files):
data = {q: data[q] for q in data.keys()}
self.data[self.out[i]] = data
- self.out = [self.out[i] for i in range(len(self.out)) if i not in self.faulty_list]
- self.dt = (self.data[self.out[0]]["t"][-1] - self.data[self.out[0]]["t"][-2])
+ self.out = [
+ self.out[i] for i in range(len(self.out)) if i not in self.faulty_list
+ ]
+ self.dt = self.data[self.out[0]]["t"][-1] - self.data[self.out[0]]["t"][-2]
def composite_data(self):
def add_to_data(fct, quantity, comp_NAME, **kwargs):
@@ -206,23 +294,44 @@ def add_to_data(fct, quantity, comp_NAME, **kwargs):
NAME[quantity] = comp_NAME
def M_Z_vs_M_a(o):
- return np.sum(self.data[o]["comp_a"][:, 0, :-7], axis=1) / self.data[o]["M_a"] * u.dimensionless_unscaled
+ return (
+ np.sum(self.data[o]["comp_a"][:, 0, :-7], axis=1)
+ / self.data[o]["M_a"]
+ * u.dimensionless_unscaled
+ )
def M_Z_vs_M(o):
- return np.sum(self.data[o]["comp_a"][:, 0, :-7] + self.data[o]["comp_c"][:, 0, :-7], axis=1) / self.data[o][
- "M"] * u.dimensionless_unscaled
+ return (
+ np.sum(
+ self.data[o]["comp_a"][:, 0, :-7]
+ + self.data[o]["comp_c"][:, 0, :-7],
+ axis=1,
+ )
+ / self.data[o]["M"]
+ * u.dimensionless_unscaled
+ )
def M_Z(o):
- return np.sum(self.data[o]["comp_a"][:, 0, :-7] + self.data[o]["comp_c"][:, 0, :-7], axis=1)
+ return np.sum(
+ self.data[o]["comp_a"][:, 0, :-7] + self.data[o]["comp_c"][:, 0, :-7],
+ axis=1,
+ )
def M_Z_vs_M_c(o):
- return np.sum(self.data[o]["comp_c"][:, 0, :-7], axis=1) / self.data[o]["M_c"] * u.dimensionless_unscaled
+ return (
+ np.sum(self.data[o]["comp_c"][:, 0, :-7], axis=1)
+ / self.data[o]["M_c"]
+ * u.dimensionless_unscaled
+ )
def elem_a(o):
return dict(
zip(
elements,
- (self.data[o]["comp_a"][:, 0, :] / self.data[o]["M_a"][:, np.newaxis]).T
+ (
+ self.data[o]["comp_a"][:, 0, :]
+ / self.data[o]["M_a"][:, np.newaxis]
+ ).T,
)
)
@@ -230,7 +339,10 @@ def elem_c(o):
return dict(
zip(
elements,
- (self.data[o]["comp_c"][:, 0, :] / self.data[o]["M_c"][:, np.newaxis]).T
+ (
+ self.data[o]["comp_c"][:, 0, :]
+ / self.data[o]["M_c"][:, np.newaxis]
+ ).T,
)
)
@@ -238,8 +350,13 @@ def elem_tot(o):
return dict(
zip(
elements,
- ((self.data[o]["comp_c"][:, 0, :] + self.data[o]["comp_a"][:, 0, :]) / self.data[o]["M"][:,
- np.newaxis]).T
+ (
+ (
+ self.data[o]["comp_c"][:, 0, :]
+ + self.data[o]["comp_a"][:, 0, :]
+ )
+ / self.data[o]["M"][:, np.newaxis]
+ ).T,
)
)
@@ -247,7 +364,10 @@ def mol_a(o):
return dict(
zip(
molecules,
- (self.data[o]["comp_a"][:, 1, :] / self.data[o]["M_a"][:, np.newaxis]).T
+ (
+ self.data[o]["comp_a"][:, 1, :]
+ / self.data[o]["M_a"][:, np.newaxis]
+ ).T,
)
)
@@ -255,7 +375,10 @@ def mol_c(o):
return dict(
zip(
molecules,
- (self.data[o]["comp_c"][:, 1, :] / self.data[o]["M_c"][:, np.newaxis]).T
+ (
+ self.data[o]["comp_c"][:, 1, :]
+ / self.data[o]["M_c"][:, np.newaxis]
+ ).T,
)
)
@@ -263,28 +386,53 @@ def mol_tot(o):
return dict(
zip(
molecules,
- ((self.data[o]["comp_c"][:, 1, :] + self.data[o]["comp_a"][:, 1, :]) / self.data[o]["M"][:,
- np.newaxis]).T
+ (
+ (
+ self.data[o]["comp_c"][:, 1, :]
+ + self.data[o]["comp_a"][:, 1, :]
+ )
+ / self.data[o]["M"][:, np.newaxis]
+ ).T,
)
)
def XH(o, X, orig):
- return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"]["H"] / self.data[o]["solar"][X]
+ return (
+ self.data[o][f"elem_{orig}"][X]
+ / self.data[o][f"elem_{orig}"]["H"]
+ / self.data[o]["solar"][X]
+ )
def XY(o, X, Y, orig):
- return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"][Y] / self.data[o]["solar"][X] * \
- self.data[o]["solar"][Y]
+ return (
+ self.data[o][f"elem_{orig}"][X]
+ / self.data[o][f"elem_{orig}"][Y]
+ / self.data[o]["solar"][X]
+ * self.data[o]["solar"][Y]
+ )
def XH_not_normed(o, X, orig):
- return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"]["H"] / el_masses[X]
+ return (
+ self.data[o][f"elem_{orig}"][X]
+ / self.data[o][f"elem_{orig}"]["H"]
+ / el_masses[X]
+ )
def XY_not_normed(o, X, Y, orig):
- return self.data[o][f"elem_{orig}"][X] / self.data[o][f"elem_{orig}"][Y] * el_masses[Y] / el_masses[X]
+ return (
+ self.data[o][f"elem_{orig}"][X]
+ / self.data[o][f"elem_{orig}"][Y]
+ * el_masses[Y]
+ / el_masses[X]
+ )
def pebiso(o):
- iso = np.where(self.data[o]["m_dot_gas"].value > 0, np.ones_like(self.data[o]["t"].value),
- np.zeros_like(self.data[o]["t"].value))
- iso[np.argmax(iso>0):] = 1
+ iso = np.where(
+ self.data[o]["m_dot_gas"].value > 0,
+ np.ones_like(self.data[o]["t"].value),
+ np.zeros_like(self.data[o]["t"].value),
+ )
+ iso[np.argmax(iso > 0) :] = 1
iso = iso * u.dimensionless_unscaled
return iso
@@ -341,8 +489,10 @@ def convert_units(self):
def get_label(self, n, quantity):
if hasattr(quantity, "unit") and quantity.unit != u.dimensionless_unscaled:
x = r"${}$ / {:latex}".format(n, quantity.unit)
- elif hasattr(self.data[self.out[0]].get(n, None), "unit") and self.data[self.out[0]][
- n].unit != u.dimensionless_unscaled:
+ elif (
+ hasattr(self.data[self.out[0]].get(n, None), "unit")
+ and self.data[self.out[0]][n].unit != u.dimensionless_unscaled
+ ):
quantity = self.data[self.out[0]][n]
x = r"${}$ / {:latex}".format(n, quantity.unit)
elif n == "a_p":
@@ -360,20 +510,20 @@ def get_label(self, n, quantity):
return x
def plot_data(
- self,
- quantities,
- out=None,
- sub_quant=None,
- t_dot=None,
- label=None,
- ulabel=True,
- fig=None,
- units=None,
- ax=None,
- plt=None,
- color_dict=None,
- use_ls=True,
- **kwargs
+ self,
+ quantities,
+ out=None,
+ sub_quant=None,
+ t_dot=None,
+ label=None,
+ ulabel=True,
+ fig=None,
+ units=None,
+ ax=None,
+ plt=None,
+ color_dict=None,
+ use_ls=True,
+ **kwargs,
):
snapshot = kwargs.get("snapshot", None)
out = self.out if out is None else out
@@ -395,17 +545,25 @@ def plot_data(
use_t_dot = False
if snapshot is None:
+
def get_dots(x_data, y_data, t):
- t_at_dot = np.arange(dotspace.cgs.value, np.max(t[:-1]).cgs.value, dotspace.cgs.value)
+ t_at_dot = np.arange(
+ dotspace.cgs.value, np.max(t[:-1]).cgs.value, dotspace.cgs.value
+ )
x = np.interp(t_at_dot, t.cgs.value, x_data)
y = np.interp(t_at_dot, t.cgs.value, y_data)
- t_at_dot = np.arange(dotspace_large.cgs.value, np.max(t[:-1]).cgs.value, dotspace_large.cgs.value)
+ t_at_dot = np.arange(
+ dotspace_large.cgs.value,
+ np.max(t[:-1]).cgs.value,
+ dotspace_large.cgs.value,
+ )
x_large = np.interp(t_at_dot, t.cgs.value, x_data)
y_large = np.interp(t_at_dot, t.cgs.value, y_data)
- x, y = u.Quantity([xi for xi in x if xi not in x_large], x.unit), u.Quantity(
- [yi for yi in y if yi not in y_large], y.unit)
+ x, y = u.Quantity(
+ [xi for xi in x if xi not in x_large], x.unit
+ ), u.Quantity([yi for yi in y if yi not in y_large], y.unit)
return x, y, x_large, y_large
if not hasattr(self, "_plot_x_data"):
@@ -434,63 +592,117 @@ def get_dots(x_data, y_data, t):
legend = label
if units is not None:
- self._plot_x_data[ax][o], self._plot_y_data[ax][o] = planet[quantities[0]].to(units[0]), planet[
- quantities[1]].to(units[1])
+ self._plot_x_data[ax][o], self._plot_y_data[ax][o] = planet[
+ quantities[0]
+ ].to(units[0]), planet[quantities[1]].to(units[1])
else:
- self._plot_x_data[ax][o], self._plot_y_data[ax][o] = planet[quantities[0]], planet[quantities[1]]
+ self._plot_x_data[ax][o], self._plot_y_data[ax][o] = (
+ planet[quantities[0]],
+ planet[quantities[1]],
+ )
if ax is not None:
-
if use_ls:
- ls_mask = self.data[o]["pebiso"]==0
+ ls_mask = self.data[o]["pebiso"] == 0
not_ls_mask = np.logical_not(ls_mask)
- not_ls_mask[np.argmin(ls_mask)-1] = True
- self._plot[ax][o], = ax.loglog(self._plot_x_data[ax][o][not_ls_mask],
- self._plot_y_data[ax][o][not_ls_mask], label=legend, ls=":",
- **kwargs)
+ not_ls_mask[np.argmin(ls_mask) - 1] = True
+ (self._plot[ax][o],) = ax.loglog(
+ self._plot_x_data[ax][o][not_ls_mask],
+ self._plot_y_data[ax][o][not_ls_mask],
+ label=legend,
+ ls=":",
+ **kwargs,
+ )
else:
ls_mask = np.arange(len(self._plot_x_data[ax][o]))
- self._plot[ax][o], = ax.loglog(self._plot_x_data[ax][o][ls_mask], self._plot_y_data[ax][o][ls_mask], label=legend,
- **kwargs)
-
-
+ (self._plot[ax][o],) = ax.loglog(
+ self._plot_x_data[ax][o][ls_mask],
+ self._plot_y_data[ax][o][ls_mask],
+ label=legend,
+ **kwargs,
+ )
- ax.set_xlabel(self.get_label(NAME.get(quantities[0], quantities[0]), self._plot_x_data[ax][o]))
- ax.set_ylabel(self.get_label(NAME.get(quantities[1], quantities[1]), self._plot_y_data[ax][o]))
+ ax.set_xlabel(
+ self.get_label(
+ NAME.get(quantities[0], quantities[0]),
+ self._plot_x_data[ax][o],
+ )
+ )
+ ax.set_ylabel(
+ self.get_label(
+ NAME.get(quantities[1], quantities[1]),
+ self._plot_y_data[ax][o],
+ )
+ )
if use_t_dot:
x, y, x_large, y_large = get_dots(
planet[quantities[0]], planet[quantities[1]], planet["t"]
)
if len(x) > 0:
- ax.scatter(x, y, marker="o", s=10, color=color_t_dot_small[:len(x)], alpha=1)
+ ax.scatter(
+ x,
+ y,
+ marker="o",
+ s=10,
+ color=color_t_dot_small[: len(x)],
+ alpha=1,
+ )
if len(x_large) > 0:
- ax.scatter(x_large, y_large, s=25, marker="o", color=color_t_dot[:len(x_large)], alpha=1)
+ ax.scatter(
+ x_large,
+ y_large,
+ s=25,
+ marker="o",
+ color=color_t_dot[: len(x_large)],
+ alpha=1,
+ )
elif plt is not None:
-
- self._plot[o], = plt.loglog(
- self._plot_x_data[ax][o], self._plot_y_data[ax][o], label=legend, **kwargs
+ (self._plot[o],) = plt.loglog(
+ self._plot_x_data[ax][o],
+ self._plot_y_data[ax][o],
+ label=legend,
+ **kwargs,
)
- plt.xlabel(self.get_label(NAME.get(quantities[0], quantities[0]), self._plot_x_data[ax][o]))
- plt.ylabel(self.get_label(NAME.get(quantities[1], quantities[1]), self._plot_y_data[ax][o]))
+ plt.xlabel(
+ self.get_label(
+ NAME.get(quantities[0], quantities[0]),
+ self._plot_x_data[ax][o],
+ )
+ )
+ plt.ylabel(
+ self.get_label(
+ NAME.get(quantities[1], quantities[1]),
+ self._plot_y_data[ax][o],
+ )
+ )
plt.title(
- self.get_label(NAME.get(quantities[1], quantities[1]), NAME.get(quantities[0], quantities[0])))
+ self.get_label(
+ NAME.get(quantities[1], quantities[1]),
+ NAME.get(quantities[0], quantities[0]),
+ )
+ )
if use_t_dot:
x, y, x_large, y_large = get_dots(
planet[quantities[0]], planet[quantities[1]], planet["t"]
)
plt.plot(x, y, "ko", markersize=4, color="gray", alpha=1)
- plt.plot(x_large, y_large, "ko", markersize=8, color="gray", alpha=1)
+ plt.plot(
+ x_large, y_large, "ko", markersize=8, color="gray", alpha=1
+ )
else:
for o in out:
if o in self.out:
snapshot = np.minimum(snapshot, len(self._plot_x_data[ax][o]))
- self._plot[ax][o].set_data(self._plot_x_data[ax][o][:snapshot], self._plot_y_data[ax][o][:snapshot])
+ self._plot[ax][o].set_data(
+ self._plot_x_data[ax][o][:snapshot],
+ self._plot_y_data[ax][o][:snapshot],
+ )
return self._plot[ax]
@@ -515,24 +727,70 @@ def apply_mask(self, mask):
self.out = mask
self.plotlabels = {k: v for (k, v) in self.plotlabels.items() if k in mask}
- def pqplot_video(self, plt, quantity, frames, lines=[], sub_qaunt=None, log_data=False, units=None, **kwargs):
- """ not working anymore, needs some more fixes"""
+ def pqplot_video(
+ self,
+ plt,
+ quantity,
+ frames,
+ lines=[],
+ sub_qaunt=None,
+ log_data=False,
+ units=None,
+ **kwargs,
+ ):
+ """not working anymore, needs some more fixes"""
fig, ax = plt.subplots(sharex=True, sharey=True)
div = make_axes_locatable(ax)
- cax = div.append_axes('right', '5%', '5%')
+ cax = div.append_axes("right", "5%", "5%")
def animate(i):
[ax.clear() for ax in fig.axes]
- self.pqplot(plt=plt, snapshot=i, fig=fig, quantity=quantity, lines=lines, sub_qaunt=sub_qaunt,
- log_data=log_data, cax=cax, units=units, **kwargs)
+ self.pqplot(
+ plt=plt,
+ snapshot=i,
+ fig=fig,
+ quantity=quantity,
+ lines=lines,
+ sub_qaunt=sub_qaunt,
+ log_data=log_data,
+ cax=cax,
+ units=units,
+ **kwargs,
+ )
return animation.FuncAnimation(fig, animate, frames=frames)
- def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=False, symlog_data=False, snapshot=-1,
- units=None, icelines=None, ax=None, shape=None, cbar_location=None,
- **kwargs):
- data_all, extend_r, p_data, p_name, q_data, q_name, r_list, r_name, x_data_all, y_data_all, z_data, z_data_all = self.extract_r0t0_data(
- quantity, snapshot, sub_quant, units)
+ def pqplot(
+ self,
+ fig,
+ quantity,
+ lines=[],
+ cax=None,
+ sub_quant=None,
+ log_data=False,
+ symlog_data=False,
+ snapshot=-1,
+ units=None,
+ icelines=None,
+ ax=None,
+ shape=None,
+ cbar_location=None,
+ **kwargs,
+ ):
+ (
+ data_all,
+ extend_r,
+ p_data,
+ p_name,
+ q_data,
+ q_name,
+ r_list,
+ r_name,
+ x_data_all,
+ y_data_all,
+ z_data,
+ z_data_all,
+ ) = self.extract_r0t0_data(quantity, snapshot, sub_quant, units)
if shape is None:
if len(extend_r) > 2:
@@ -569,20 +827,32 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal
for l in lines:
if isinstance(l[0], str):
Z = np.array(
- [d[l[0]].value[snapshot if snapshot < len(d["t"]) else -1] for d in data_all[i].values()])
+ [
+ d[l[0]].value[snapshot if snapshot < len(d["t"]) else -1]
+ for d in data_all[i].values()
+ ]
+ )
else:
Z = l[0]
- zi, yi, xi = np.histogram2d(p_data, q_data, bins=(len(np.unique(p_data)), len(np.unique(q_data))),
- weights=Z,
- normed=False)
+ zi, yi, xi = np.histogram2d(
+ p_data,
+ q_data,
+ bins=(len(np.unique(p_data)), len(np.unique(q_data))),
+ weights=Z,
+ normed=False,
+ )
# zi = gaussian_filter(zi.T, 0.15)
x_step = l[1].pop("x_step", 1)
y_step = l[1].pop("y_step", 1)
zi = zi.T
zi = np.ma.masked_equal(zi, np.inf)
- CS = ax_temp.contour(x_data_all[i][::x_step, ::y_step], y_data_all[i][::x_step, ::y_step],
- zi[::x_step, ::y_step], **l[1])
+ CS = ax_temp.contour(
+ x_data_all[i][::x_step, ::y_step],
+ y_data_all[i][::x_step, ::y_step],
+ zi[::x_step, ::y_step],
+ **l[1],
+ )
ax_temp.clabel(CS, l[2])
l[1]["x_step"] = x_step
l[1]["y_step"] = y_step
@@ -593,11 +863,9 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal
time = 3 * u.Myr if snapshot == -1 else snapshot * self.dt
if r_name is not None:
if r_name == "case":
- ax_temp.set_title(
- r"{}".format(r_list[i]))
+ ax_temp.set_title(r"{}".format(r_list[i]))
else:
- ax_temp.set_title(
- r"{}={}".format(r_name, r_list[i]))
+ ax_temp.set_title(r"{}={}".format(r_name, r_list[i]))
shrink = 0.4 if cbar_location == "bottom" else 0.6
aspect = 6 if cbar_location == "bottom" else 20
@@ -607,7 +875,9 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal
cbar = fig.colorbar(im, ax=ax[:-1], cax=cax, shrink=shrink, aspect=aspect)
ax.flat[-1].remove()
else:
- cbar = fig.colorbar(im, ax=ax, cax=cax, location=cbar_location, shrink=shrink, aspect=aspect)
+ cbar = fig.colorbar(
+ im, ax=ax, cax=cax, location=cbar_location, shrink=shrink, aspect=aspect
+ )
try:
for temp_ax in ax.flat:
@@ -615,15 +885,17 @@ def pqplot(self, fig, quantity, lines=[], cax=None, sub_quant=None, log_data=Fal
except Exception as e:
print(f"no label outer possible {e}")
-
cbar.ax.set_title(self.get_label(quantity, z_data), y=1.03)
fig.suptitle(
- r"${}$-${}$ Plot for {}, t = {}".format(p_name, q_name, quantity, time.to(u.Myr)))
+ r"${}$-${}$ Plot for {}, t = {}".format(
+ p_name, q_name, quantity, time.to(u.Myr)
+ )
+ )
return im
def extract_r0t0_data(self, quantity, snapshot, sub_quant, units):
- """ data in pipline should look like: r0 = [20, 10, 5, 20, 10, 5]), t0 = [1e5, 1e5, 1e5, 1e6, 1e6, 1e6]) - in other words: flattened meshgrids"""
+ """data in pipline should look like: r0 = [20, 10, 5, 20, 10, 5]), t0 = [1e5, 1e5, 1e5, 1e6, 1e6, 1e6]) - in other words: flattened meshgrids"""
format_string = "{}:{}, {}:{}"
format_string_long = "{}:{}, {}:{}, {}:{}"
if parse.parse(format_string_long, self.plotlabels[self.out[0]]) is not None:
@@ -637,7 +909,10 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units):
if np.all(unique == np.array(["evaporation", "plain"])):
unique = unique[::-1]
- extend_r = [[self.out[i] for i in range(len(self.out)) if pa[i] == lsta] for lsta in unique]
+ extend_r = [
+ [self.out[i] for i in range(len(self.out)) if pa[i] == lsta]
+ for lsta in unique
+ ]
elif parse.parse(format_string, self.plotlabels[self.out[0]]) is not None:
# -> case of only par_1, par_2
extend_r = [self.out]
@@ -650,9 +925,13 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units):
p_data, q_data = [], []
for o in sub:
if len(extend_r) > 1:
- r_name, r, p_name, p, q_name, q = parse.parse(format_string_long, self.plotlabels[o])
+ r_name, r, p_name, p, q_name, q = parse.parse(
+ format_string_long, self.plotlabels[o]
+ )
else:
- p_name, p, q_name, q = parse.parse(format_string, self.plotlabels[o])
+ p_name, p, q_name, q = parse.parse(
+ format_string, self.plotlabels[o]
+ )
r_name = None
r = 0
p_data.append(float(p))
@@ -665,20 +944,35 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units):
q_data = (q_data * units[1]).cgs.to(units[1])
if sub_quant is not None:
- z_data = np.array([data[quantity][sub_quant][snapshot if snapshot < len(data["t"]) else -1] for data in
- data.values()])
+ z_data = np.array(
+ [
+ data[quantity][sub_quant][
+ snapshot if snapshot < len(data["t"]) else -1
+ ]
+ for data in data.values()
+ ]
+ )
else:
z_data = np.array(
- [data[quantity].value[snapshot if snapshot < len(data["t"]) else -1] for data in
- data.values()])
+ [
+ data[quantity].value[
+ snapshot if snapshot < len(data["t"]) else -1
+ ]
+ for data in data.values()
+ ]
+ )
if units is not None:
z_data = z_data.to(units[2])
- zi, yi, xi = np.histogram2d(p_data, q_data, bins=(len(np.unique(p_data)), len(np.unique(q_data))),
- weights=z_data,
- normed=False)
+ zi, yi, xi = np.histogram2d(
+ p_data,
+ q_data,
+ bins=(len(np.unique(p_data)), len(np.unique(q_data))),
+ weights=z_data,
+ normed=False,
+ )
zi = np.ma.masked_equal(zi, 0)
zi = np.ma.masked_equal(zi, np.inf)
X, Y = np.meshgrid(np.unique(p_data), np.unique(q_data))
@@ -689,35 +983,86 @@ def extract_r0t0_data(self, quantity, snapshot, sub_quant, units):
data_all.append(data)
x_data_all = np.array(x_data_all)
y_data_all = np.array(y_data_all)
- return data_all, extend_r, p_data, p_name, q_data, q_name, r_list, r_name, x_data_all, y_data_all, z_data, z_data_all
+ return (
+ data_all,
+ extend_r,
+ p_data,
+ p_name,
+ q_data,
+ q_name,
+ r_list,
+ r_name,
+ x_data_all,
+ y_data_all,
+ z_data,
+ z_data_all,
+ )
def get_iceline_background(self, file, ax):
def get_data_iceline():
data = {}
- with tables.open_file(file, mode='r') as f:
+ with tables.open_file(file, mode="r") as f:
data["T"] = np.array(f.root.disk.T)
data["r"] = np.squeeze(np.array(f.root.disk.r) / const.au.cgs.value)
- data["t"] = np.squeeze(np.array(f.root.disk.t)) / 1e6 / 365.25 / 24 / 3600
+ data["t"] = (
+ np.squeeze(np.array(f.root.disk.t)) / 1e6 / 365.25 / 24 / 3600
+ )
- data["icelines"] = np.array([20, 30, 70, 150, 371, 970, 1354, 1423, 1500, 1653, 2000])
+ data["icelines"] = np.array(
+ [20, 30, 70, 150, 371, 970, 1354, 1423, 1500, 1653, 2000]
+ )
data["iceline_pos"] = np.array(
- [np.searchsorted(-data["T"][i], -data["icelines"]) for i in range(len(data["T"]))])
- data["iceline_r"] = np.array([data["r"][data["iceline_pos"][i]] for i in range(len(data["t"]))])
- data["iceline_name"] = ["CO & N2", "CH4", "CO2", "H2O & H2S", "Fe3O4", "NaAlSi3O8 & KAlSi3O8",
- "Mg2SiO4 & Fe2O3",
- "VO", "MgSiO3", "Al2O3", "TiO"]
+ [
+ np.searchsorted(-data["T"][i], -data["icelines"])
+ for i in range(len(data["T"]))
+ ]
+ )
+ data["iceline_r"] = np.array(
+ [data["r"][data["iceline_pos"][i]] for i in range(len(data["t"]))]
+ )
+ data["iceline_name"] = [
+ "CO & N2",
+ "CH4",
+ "CO2",
+ "H2O & H2S",
+ "Fe3O4",
+ "NaAlSi3O8 & KAlSi3O8",
+ "Mg2SiO4 & Fe2O3",
+ "VO",
+ "MgSiO3",
+ "Al2O3",
+ "TiO",
+ ]
return data
data = get_data_iceline()
xlim, ylim = ax.get_xlim(), ax.get_ylim()
texts = []
for i in range(len(data["icelines"])):
- if np.any(np.logical_and(data["iceline_r"][:, i] > xlim[0], data["iceline_r"][:, i] < xlim[1])):
+ if np.any(
+ np.logical_and(
+ data["iceline_r"][:, i] > xlim[0], data["iceline_r"][:, i] < xlim[1]
+ )
+ ):
ymask = np.logical_and(data["t"] > ylim[0], data["t"] < ylim[1])
- ax.plot(data["iceline_r"][ymask, i], data["t"][ymask], ls="--", color="gray", alpha=0.8)
- t = ax.text(np.mean(data["iceline_r"][ymask, i]), 0.5 * (np.mean(ax.get_ylim())),
- data["iceline_name"][i], fontsize=10,
- rotation=90, ha="center", va="center", color="white", backgroundcolor="gray")
+ ax.plot(
+ data["iceline_r"][ymask, i],
+ data["t"][ymask],
+ ls="--",
+ color="gray",
+ alpha=0.8,
+ )
+ t = ax.text(
+ np.mean(data["iceline_r"][ymask, i]),
+ 0.5 * (np.mean(ax.get_ylim())),
+ data["iceline_name"][i],
+ fontsize=10,
+ rotation=90,
+ ha="center",
+ va="center",
+ color="white",
+ backgroundcolor="gray",
+ )
t.set_bbox(dict(alpha=0.0))
texts.append(t)
diff --git a/chemcomp/helpers/import_config.py b/chemcomp/helpers/import_config.py
index e72d76d..bb64703 100755
--- a/chemcomp/helpers/import_config.py
+++ b/chemcomp/helpers/import_config.py
@@ -1,5 +1,6 @@
import astropy.units as u
from yaml import load, dump
+
try:
from yaml import CLoader as Loader, CDumper as Dumper
except ImportError:
@@ -7,7 +8,7 @@
def import_config(config_file):
- """Wrapper that imports the config. """
+ """Wrapper that imports the config."""
stream = open(config_file, "r")
config = load(stream, Loader=Loader)
@@ -17,7 +18,7 @@ def import_config(config_file):
def scale(dictionary, quantity, unit=1):
- """ Helper function to scale dictionary with unit """
+ """Helper function to scale dictionary with unit"""
if dictionary is not None:
if isinstance(dictionary, list):
for item in dictionary:
@@ -31,15 +32,15 @@ def scale(dictionary, quantity, unit=1):
def scale_units(config):
- """Add a unit to all unit based quantities from the configuration file. """
+ """Add a unit to all unit based quantities from the configuration file."""
planets = config.get("config_planet", None)
scale(planets, "M_c", u.earthMass)
scale(planets, "M_a", u.earthMass)
scale(planets, "a_p", u.au)
scale(planets, "t_0", u.Myr)
scale(planets, "r_p", u.jupiterRad)
- scale(planets, "rho_c", u.g / (u.cm ** 3))
- scale(planets, "rho_0", u.g / (u.cm ** 3))
+ scale(planets, "rho_c", u.g / (u.cm**3))
+ scale(planets, "rho_0", u.g / (u.cm**3))
scale(planets, "dt", u.yr)
scale(planets, "M_end", u.jupiterMass)
scale(planets, "r_in", u.au)
@@ -50,7 +51,7 @@ def scale_units(config):
scale(disk, "M0", u.M_sun)
scale(disk, "R0", u.au)
scale(disk, "time", u.Myr)
- scale(disk, "SIGMA_0", u.g / (u.cm ** 2))
+ scale(disk, "SIGMA_0", u.g / (u.cm**2))
scale(disk, "a_0", u.cm)
scale(disk, "FUV_MDOT", u.M_sun / u.yr)
scale(disk, "tau_disk", u.yr)
@@ -60,11 +61,11 @@ def scale_units(config):
pebble_accretion = config.get("config_pebble_accretion", None)
scale(pebble_accretion, "M_DOT_P", u.earthMass / (u.Myr))
scale(pebble_accretion, "R_peb", u.cm)
- scale(pebble_accretion, "rho_solid", u.g / (u.cm ** 3))
+ scale(pebble_accretion, "rho_solid", u.g / (u.cm**3))
scale(pebble_accretion, "u_frag", u.m / u.s)
gas_accretion = config.get("config_gas_accretion", None)
- scale(gas_accretion, "kappa_env", u.cm ** 2 / u.g)
+ scale(gas_accretion, "kappa_env", u.cm**2 / u.g)
defaults = config.get("defaults", None)
scale(defaults, "DEF_R_IN", u.au)
diff --git a/chemcomp/helpers/main_helpers.py b/chemcomp/helpers/main_helpers.py
index bf408b1..21ca67a 100755
--- a/chemcomp/helpers/main_helpers.py
+++ b/chemcomp/helpers/main_helpers.py
@@ -4,6 +4,7 @@
import numpy as np
from slugify import slugify
+
try:
from yaml import CLoader as Loader, CDumper as Dumper
except ImportError:
@@ -39,7 +40,7 @@ def file_checks_out(config_file, continue_job):
if not continue_job:
return True
- output = import_config(config_file).get("output",{})
+ output = import_config(config_file).get("output", {})
name = output.get("name", "planet")
filename = output.get("directory", os.path.join(OUTPUT_PATH, "{}.h5".format(name)))
@@ -94,11 +95,23 @@ def create_pipeline(job, config):
values = np.array([mesh.flatten() for mesh in parameter_mesh])
if double_disks:
- files, names = double_disks_wrapper(template, values, parameter_info, name=name,
- save_disk=save_disk, save_interval=save_interval)
+ files, names = double_disks_wrapper(
+ template,
+ values,
+ parameter_info,
+ name=name,
+ save_disk=save_disk,
+ save_interval=save_interval,
+ )
else:
- files = create_yaml_test_range(template, values, parameter_info, name=name, save_disk=save_disk,
- save_interval=save_interval)
+ files = create_yaml_test_range(
+ template,
+ values,
+ parameter_info,
+ name=name,
+ save_disk=save_disk,
+ save_interval=save_interval,
+ )
names = [name]
return files, names
@@ -107,7 +120,7 @@ def create_pipeline(job, config):
def run_setup(config_file):
"""load config and start modelling"""
if not os.path.isdir(OUTPUT_PATH):
- print(f'creating Folder {OUTPUT_PATH}')
+ print(f"creating Folder {OUTPUT_PATH}")
os.mkdir(OUTPUT_PATH)
# Load config
@@ -131,7 +144,9 @@ def run_setup(config_file):
planet_model = "BertPlanet"
print("No planetmodel specified, using BertPlanet")
- disk = getattr(chemcomp.disks, disk_model)(defaults, chemistry=chemistry, **config_disk)
+ disk = getattr(chemcomp.disks, disk_model)(
+ defaults, chemistry=chemistry, **config_disk
+ )
accretion_models = get_acc_models(
config_pebble_accretion,
@@ -152,12 +167,16 @@ def run_setup(config_file):
def zip_all(names=[], delete=False):
if not os.path.isdir(ZIP_OUTPUT_PATH):
- print(f'creating Folder {ZIP_OUTPUT_PATH}')
+ print(f"creating Folder {ZIP_OUTPUT_PATH}")
os.mkdir(ZIP_OUTPUT_PATH)
zipf = zipfile.ZipFile(
- os.path.join(ZIP_OUTPUT_PATH,
- "{}_{}.zip".format(os.path.commonprefix(names), datetime.now().strftime("%Y%m%d-%H%M%S"))),
+ os.path.join(
+ ZIP_OUTPUT_PATH,
+ "{}_{}.zip".format(
+ os.path.commonprefix(names), datetime.now().strftime("%Y%m%d-%H%M%S")
+ ),
+ ),
"w",
zipfile.ZIP_DEFLATED,
)
@@ -198,7 +217,15 @@ def is_defined(conf, acc):
]
-def write_yaml(template, parameter_info, values, template_name, name="", save_disk=False, save_interval=None):
+def write_yaml(
+ template,
+ parameter_info,
+ values,
+ template_name,
+ name="",
+ save_disk=False,
+ save_interval=None,
+):
if isinstance(template, str):
with open(os.path.join(CONFIG_PATH, "{}.yaml".format(template))) as f:
config_template = load(f, Loader=Loader)
@@ -235,7 +262,9 @@ def write_yaml(template, parameter_info, values, template_name, name="", save_di
return new_file_name, slug, plot_name
-def create_yaml_test_range(template, values, parameter_info, name="", save_disk=False, save_interval=None):
+def create_yaml_test_range(
+ template, values, parameter_info, name="", save_disk=False, save_interval=None
+):
files, slugs, plot_names = [], [], []
if not isinstance(template, str):
@@ -244,7 +273,7 @@ def create_yaml_test_range(template, values, parameter_info, name="", save_disk=
template_name = template
if not os.path.isdir(CONFIG_PATH):
- print(f'creating Folder {CONFIG_PATH}')
+ print(f"creating Folder {CONFIG_PATH}")
os.mkdir(CONFIG_PATH)
path = os.path.join(CONFIG_PATH, "p_{}".format(name))
@@ -255,9 +284,15 @@ def create_yaml_test_range(template, values, parameter_info, name="", save_disk=
os.mkdir(path)
for i in range(len(values.T)):
- file, slug, plot_name = write_yaml(template, parameter_info, values.T[i],
- template_name=template_name, name=name,
- save_disk=save_disk, save_interval=save_interval)
+ file, slug, plot_name = write_yaml(
+ template,
+ parameter_info,
+ values.T[i],
+ template_name=template_name,
+ name=name,
+ save_disk=save_disk,
+ save_interval=save_interval,
+ )
files.append(file)
slugs.append(slug)
plot_names.append(plot_name)
@@ -266,14 +301,16 @@ def create_yaml_test_range(template, values, parameter_info, name="", save_disk=
slugfile.write("\n".join(slugs))
with open(
- os.path.join(OUTPUT_PATH, "plot_names_{}.dat".format(name)), "w"
+ os.path.join(OUTPUT_PATH, "plot_names_{}.dat".format(name)), "w"
) as plot_name_file:
plot_name_file.write("\n".join(plot_names))
return files
-def double_disks_wrapper(template, values, parameter_info, name, save_disk, save_interval):
+def double_disks_wrapper(
+ template, values, parameter_info, name, save_disk, save_interval
+):
"""
change template and create four different disks: With Pla, With evap, With Pla+evap, plain
Note: currently only two disks. Thank you referee!
@@ -288,14 +325,22 @@ def double_disks_wrapper(template, values, parameter_info, name, save_disk, save
template_None = config_template.copy()
template_None["config_disk"]["DTG_pla"] = 0.0
template_None["config_disk"]["evaporation"] = False
- files.append(create_yaml_test_range(template_None, values, parameter_info, name, save_disk, save_interval))
+ files.append(
+ create_yaml_test_range(
+ template_None, values, parameter_info, name, save_disk, save_interval
+ )
+ )
# evap:
template_evap = config_template.copy()
template_evap["config_disk"]["DTG_pla"] = 0.0
template_evap["config_disk"]["evaporation"] = True
name_evap = f"{name}_evap"
- files.append(create_yaml_test_range(template_evap, values, parameter_info, name_evap, save_disk, save_interval))
+ files.append(
+ create_yaml_test_range(
+ template_evap, values, parameter_info, name_evap, save_disk, save_interval
+ )
+ )
return [item for sublist in files for item in sublist], [name_evap, name]
diff --git a/chemcomp/helpers/solvediffonedee.py b/chemcomp/helpers/solvediffonedee.py
index 66b1e8d..2bdf9b3 100755
--- a/chemcomp/helpers/solvediffonedee.py
+++ b/chemcomp/helpers/solvediffonedee.py
@@ -1,8 +1,21 @@
from .tridag import *
-def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False,
- extracond=None, retall=False):
+def solvediffonedee_components(
+ x,
+ y,
+ v,
+ d,
+ g,
+ s,
+ bcl,
+ bcr,
+ dt=None,
+ int=False,
+ upwind=False,
+ extracond=None,
+ retall=False,
+):
"""
Vectorised version of solvediffonedee (works with sigma_g_components_mol)
@@ -121,7 +134,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u
#
# Apply it to y
#
- a[0] = 0.
+ a[0] = 0.0
b[0] = bcl[1] - bcl[0] / dxi[0]
c[0] = bcl[0] / dxi[0]
r[0] = bcl[2]
@@ -129,7 +142,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u
#
# Apply it to y/g
#
- a[0] = 0.
+ a[0] = 0.0
b[0] = (bcl[1] - bcl[0] / dxi[0]) / g[0]
c[0] = (bcl[0] / dxi[0]) / g[1]
r[0] = bcl[2]
@@ -142,7 +155,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u
#
a[-1] = -bcr[0] / dxi[-1]
b[-1] = bcr[1] + bcr[0] / dxi[-1]
- c[-1] = 0.
+ c[-1] = 0.0
r[-1] = bcr[2]
else:
#
@@ -150,7 +163,7 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u
#
a[-1] = (-bcr[0] / dxi[-1]) / g[-2]
b[-1] = (bcr[1] + bcr[0] / dxi[-1]) / g[-1]
- c[-1] = 0.
+ c[-1] = 0.0
r[-1] = bcr[2]
#
# Internal conditions, if present
@@ -196,7 +209,21 @@ def solvediffonedee_components(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, u
return sol
-def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False, extracond=None, retall=False):
+def solvediffonedee(
+ x,
+ y,
+ v,
+ d,
+ g,
+ s,
+ bcl,
+ bcr,
+ dt=None,
+ int=False,
+ upwind=False,
+ extracond=None,
+ retall=False,
+):
"""
Solve the 1-D linear advection-diffusion equation with boundary conditions,
either in stationary state, or as an implicit time step.
@@ -312,7 +339,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False
#
# Apply it to y
#
- a[0] = 0.
+ a[0] = 0.0
b[0] = bcl[1] - bcl[0] / dxi[0]
c[0] = bcl[0] / dxi[0]
r[0] = bcl[2]
@@ -320,7 +347,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False
#
# Apply it to y/g
#
- a[0] = 0.
+ a[0] = 0.0
b[0] = (bcl[1] - bcl[0] / dxi[0]) / g[0]
c[0] = (bcl[0] / dxi[0]) / g[1]
r[0] = bcl[2]
@@ -333,7 +360,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False
#
a[-1] = -bcr[0] / dxi[-1]
b[-1] = bcr[1] + bcr[0] / dxi[-1]
- c[-1] = 0.
+ c[-1] = 0.0
r[-1] = bcr[2]
else:
#
@@ -341,7 +368,7 @@ def solvediffonedee(x, y, v, d, g, s, bcl, bcr, dt=None, int=False, upwind=False
#
a[-1] = (-bcr[0] / dxi[-1]) / g[-2]
b[-1] = (bcr[1] + bcr[0] / dxi[-1]) / g[-1]
- c[-1] = 0.
+ c[-1] = 0.0
r[-1] = bcr[2]
#
# Internal conditions, if present
@@ -439,7 +466,9 @@ def getfluxonedee(x, y, v, d, g, int=False, upwind=False, oned=False):
# Now compute the flux
#
if not oned:
- flux = b[:, np.newaxis, np.newaxis] * y[:-1] + c[:, np.newaxis, np.newaxis] * y[1:]
+ flux = (
+ b[:, np.newaxis, np.newaxis] * y[:-1] + c[:, np.newaxis, np.newaxis] * y[1:]
+ )
else:
flux = b * y[:-1] + c * y[1:]
#
diff --git a/chemcomp/helpers/tridag.py b/chemcomp/helpers/tridag.py
index 24854fe..caef300 100755
--- a/chemcomp/helpers/tridag.py
+++ b/chemcomp/helpers/tridag.py
@@ -1,10 +1,11 @@
import numpy as np
from scipy.linalg import solve_banded
+
def tridag_components(a, b, c, r):
"""
- The Numerical Recipes tridag() routine for Python, using the SciPy library.
- This is just a wrapper around the solve_banded function of SciPy, to make it
+ The Numerical Recipes tridag() routine for Python, using the SciPy library.
+ This is just a wrapper around the solve_banded function of SciPy, to make it
work like the tridag routine.
"""
n = r.shape
@@ -13,8 +14,7 @@ def tridag_components(a, b, c, r):
ab[1, :] = b[:]
ab[2, :-1] = a[1:]
- u = np.transpose(
- [solve_banded((1, 1), ab[:, :, i], r[:, i]) for i in range(n[1])])
+ u = np.transpose([solve_banded((1, 1), ab[:, :, i], r[:, i]) for i in range(n[1])])
return u
@@ -43,8 +43,7 @@ def pentdag(aa, a, b, c, cc, r):
ab[0, 2:] = cc[:-2]
ab[1, 1:] = c[:-1]
ab[2, :] = b[:]
- ab[3,:-1] = a[1:]
- ab[4,:-2] = aa[2:]
- u = solve_banded((2,2),ab,r)
+ ab[3, :-1] = a[1:]
+ ab[4, :-2] = aa[2:]
+ u = solve_banded((2, 2), ab, r)
return u
-
diff --git a/chemcomp/helpers/units/chemistry_const.py b/chemcomp/helpers/units/chemistry_const.py
index 62f415e..0653bfb 100755
--- a/chemcomp/helpers/units/chemistry_const.py
+++ b/chemcomp/helpers/units/chemistry_const.py
@@ -66,9 +66,9 @@
MassC = 12.0 # C in terms of H
# Masses in terms of H (atomic unit)
-MassCO = MassC + MassO #28.0 # CO mass in terms of H
+MassCO = MassC + MassO # 28.0 # CO mass in terms of H
MassCH4 = MassC + 4 * MassH # CH4 in terms of H
-MassCO2 = MassC + 2 * MassO # CO2 in terms of H
+MassCO2 = MassC + 2 * MassO # CO2 in terms of H
MassH2O = MassO + 2 * MassH # H20 in terms of H
MassFe2O3 = 2 * MassFe + 3 * MassO
@@ -78,7 +78,7 @@
MassMg2SiO4 = 2 * MassMg + MassSi + 4 * MassO # MgSiO3 in terms of H
MassNH3 = 3 * MassH + MassN
MassN2 = 2 * MassN
-MassH2S = 2 * MassH + MassS
+MassH2S = 2 * MassH + MassS
MassNaAlSi3O8 = MassNa + MassAl + 3 * MassSi + 8 * MassO
MassKAlSi3O8 = MassK + MassAl + 3 * MassSi + 8 * MassO
MassTiO = MassTi + MassO
diff --git a/chemcomp/planets/_planet_class.py b/chemcomp/planets/_planet_class.py
index 1449cf5..236a6f4 100755
--- a/chemcomp/planets/_planet_class.py
+++ b/chemcomp/planets/_planet_class.py
@@ -25,43 +25,115 @@ def __init__(self, planet, output):
self.conf_save_disk = output.get("save_disk", False)
self.conf_save_disk_interval = output.get("save_disk_interval", 100)
self._save_counter = 0
- filename = output.get("directory", os.path.join(OUTPUT_PATH, "{}.h5".format(name)))
+ filename = output.get(
+ "directory", os.path.join(OUTPUT_PATH, "{}.h5".format(name))
+ )
if os.path.exists(filename):
os.remove(filename)
- self._h5file = tables.open_file(filename, mode='w', title="Output of {}".format(name))
+ self._h5file = tables.open_file(
+ filename, mode="w", title="Output of {}".format(name)
+ )
self._h5file.set_node_attr("/", "finished", 0)
- self.output_disk = self._h5file.create_group("/", 'disk', 'Output of the Disk')
- self.output_planet = self._h5file.create_group("/", 'planet', 'Output of the Planet')
- self.output_acc = self._h5file.create_group("/", 'acc', 'Output of accretion')
+ self.output_disk = self._h5file.create_group("/", "disk", "Output of the Disk")
+ self.output_planet = self._h5file.create_group(
+ "/", "planet", "Output of the Planet"
+ )
+ self.output_acc = self._h5file.create_group("/", "acc", "Output of accretion")
self.planet = planet
self.disk = planet.disk
- list_of_disk_quantities = output.get("quantities",
- ["a_1", "t", "sigma_g", "sigma_dust", "sigma_dust_components", "T",
- "sigma_g_components",
- "f_m", "stokes_number_pebbles", "stokes_number_df", "stokes_number_drift",
- "stokes_number_frag", "stokes_number_small", "pebble_flux",
- "cum_pebble_flux",
- "peb_iso", "peb_iso_pi_crit", "T_visc", "T_irr", "mu", "m_dot", "vr_dust", "vr_gas", "m_dot_components"])
-
- self.dict_of_planet_units = {"t": "u.s", "M": "u.g", "M_a": "u.g", "M_c": "u.g", "a_p": "u.cm",
- "T": "u.K", "comp_a": "u.g", "comp_c": "u.g", "tau_m": "u.s",
- "sigma_g": "u.g/u.cm**2",
- "gamma_tot": "u.g*u.cm**2/u.s**2", "sigma_peb": "u.g/u.cm**2",
- "pebble_flux": "u.g/u.s", "peb_iso": "u.g"}
- self.dict_of_acc_units = {"m_dot": "u.g/u.s", "m_dot_a": "u.g/u.s", "m_dot_c": "u.g/u.s",
- "m_dot_a_chem": "u.g/u.s", "m_dot_c_chem": "u.g/u.s", "M_z": "u.g"}
-
- self.list_of_planet_quantities = ["t", "M", "M_a", "M_c", "a_p", "T", "comp_a", "comp_c", "tau_m",
- "gamma_tot", "sigma_g", "gamma_norm", "fSigma", "sigma_peb", "pebble_flux",
- "peb_iso"]
- self.list_of_acc_quantities = ["m_dot", "m_dot_a", "m_dot_c", "m_dot_a_chem", "m_dot_c_chem", "M_z", "regime"]
+ list_of_disk_quantities = output.get(
+ "quantities",
+ [
+ "a_1",
+ "t",
+ "sigma_g",
+ "sigma_dust",
+ "sigma_dust_components",
+ "T",
+ "sigma_g_components",
+ "f_m",
+ "stokes_number_pebbles",
+ "stokes_number_df",
+ "stokes_number_drift",
+ "stokes_number_frag",
+ "stokes_number_small",
+ "pebble_flux",
+ "cum_pebble_flux",
+ "peb_iso",
+ "peb_iso_pi_crit",
+ "T_visc",
+ "T_irr",
+ "mu",
+ "m_dot",
+ "vr_dust",
+ "vr_gas",
+ "m_dot_components",
+ ],
+ )
+
+ self.dict_of_planet_units = {
+ "t": "u.s",
+ "M": "u.g",
+ "M_a": "u.g",
+ "M_c": "u.g",
+ "a_p": "u.cm",
+ "T": "u.K",
+ "comp_a": "u.g",
+ "comp_c": "u.g",
+ "tau_m": "u.s",
+ "sigma_g": "u.g/u.cm**2",
+ "gamma_tot": "u.g*u.cm**2/u.s**2",
+ "sigma_peb": "u.g/u.cm**2",
+ "pebble_flux": "u.g/u.s",
+ "peb_iso": "u.g",
+ }
+ self.dict_of_acc_units = {
+ "m_dot": "u.g/u.s",
+ "m_dot_a": "u.g/u.s",
+ "m_dot_c": "u.g/u.s",
+ "m_dot_a_chem": "u.g/u.s",
+ "m_dot_c_chem": "u.g/u.s",
+ "M_z": "u.g",
+ }
+
+ self.list_of_planet_quantities = [
+ "t",
+ "M",
+ "M_a",
+ "M_c",
+ "a_p",
+ "T",
+ "comp_a",
+ "comp_c",
+ "tau_m",
+ "gamma_tot",
+ "sigma_g",
+ "gamma_norm",
+ "fSigma",
+ "sigma_peb",
+ "pebble_flux",
+ "peb_iso",
+ ]
+ self.list_of_acc_quantities = [
+ "m_dot",
+ "m_dot_a",
+ "m_dot_c",
+ "m_dot_a_chem",
+ "m_dot_c_chem",
+ "M_z",
+ "regime",
+ ]
# trim list of quantities to existing:
_, self.list_of_disk_quantities = [], []
- [self.list_of_disk_quantities.append(quant) if hasattr(self.disk, f"{quant}") else _.append(quant) for quant in
- list_of_disk_quantities]
+ [
+ self.list_of_disk_quantities.append(quant)
+ if hasattr(self.disk, f"{quant}")
+ else _.append(quant)
+ for quant in list_of_disk_quantities
+ ]
self._h5file.create_earray(self.output_disk, "r", obj=self.disk.r)
self._h5file.create_earray(self.output_disk, "r_i", obj=self.disk.r_i)
@@ -69,12 +141,16 @@ def __init__(self, planet, output):
self.list_of_disk_pointer = []
for quantity in self.list_of_disk_quantities:
obj = np.array(getattr(self.disk, quantity))[np.newaxis]
- self.list_of_disk_pointer.append(self._h5file.create_earray(self.output_disk, quantity, obj=obj))
+ self.list_of_disk_pointer.append(
+ self._h5file.create_earray(self.output_disk, quantity, obj=obj)
+ )
self.list_of_planet_pointer = []
for quantity in self.list_of_planet_quantities:
obj = np.array(getattr(self.planet, quantity))[np.newaxis]
- self.list_of_planet_pointer.append(self._h5file.create_earray(self.output_planet, quantity, obj=obj))
+ self.list_of_planet_pointer.append(
+ self._h5file.create_earray(self.output_planet, quantity, obj=obj)
+ )
acc_name_dict = {"PebbleAccretion": "peb", "GasAccretion": "gas"}
# drop acc model if not used:
@@ -90,12 +166,18 @@ def __init__(self, planet, output):
for i in self.acc_index_list:
self.list_of_acc_pointer[i] = []
for quantity in self.list_of_acc_quantities:
- obj = np.array(getattr(self.planet.accretion_models[i], quantity))[np.newaxis]
+ obj = np.array(getattr(self.planet.accretion_models[i], quantity))[
+ np.newaxis
+ ]
self.list_of_acc_pointer[i].append(
- self._h5file.create_earray(self.output_planet, "{}_{}".format(quantity, self.acc_models[i]),
- obj=obj))
-
- self._planet_units = self._h5file.create_group("/planet", 'units', 'units')
+ self._h5file.create_earray(
+ self.output_planet,
+ "{}_{}".format(quantity, self.acc_models[i]),
+ obj=obj,
+ )
+ )
+
+ self._planet_units = self._h5file.create_group("/planet", "units", "units")
for k, v in self.dict_of_planet_units.items():
self._h5file.set_node_attr("/planet/units", k, v)
@@ -119,13 +201,19 @@ def save_acc(self):
for j in self.acc_index_list:
ptr = self.list_of_acc_pointer[j]
for i in range(len(ptr)):
- array = getattr(self.planet.accretion_models[j], self.list_of_acc_quantities[i])
+ array = getattr(
+ self.planet.accretion_models[j], self.list_of_acc_quantities[i]
+ )
array = np.array(array)
self.list_of_acc_pointer[j][i].append(array[np.newaxis])
def save(self, force=False):
if self.planet.t >= self.planet.save_interval * self._save_counter or force:
- if self.conf_save_disk and self._save_counter % self.conf_save_disk_interval == 0 or force:
+ if (
+ self.conf_save_disk
+ and self._save_counter % self.conf_save_disk_interval == 0
+ or force
+ ):
self.save_disk()
self.save_planet()
self.save_acc()
@@ -138,10 +226,19 @@ def finish(self, error=None):
class Planet(object):
- """ docstring for Planet. """
+ """docstring for Planet."""
def __init__(
- self, output, rho_c, a_p, t_0, disk, accretion_models=[], M_c=None, M_a=None, **kwargs
+ self,
+ output,
+ rho_c,
+ a_p,
+ t_0,
+ disk,
+ accretion_models=[],
+ M_c=None,
+ M_a=None,
+ **kwargs,
):
self.save_interval = output.get("save_interval", 50000 * u.yr).cgs.value
self._print_number = 0
@@ -161,53 +258,77 @@ def __init__(
if M_c is None:
# use Lambrecht2012 transitionmass
# -> ignored if pebiso_start
- M_c = (1-self.solid_frac)* (
- np.sqrt(1.0 / 3.0)
- * (disk.eta * disk.v_k) ** 3.0
- / (const.G.cgs.value * disk.omega_k)
+ M_c = (1 - self.solid_frac) * (
+ np.sqrt(1.0 / 3.0)
+ * (disk.eta * disk.v_k) ** 3.0
+ / (const.G.cgs.value * disk.omega_k)
)
M_c = M0_fact * np.interp(a_p.cgs.value, disk.r, M_c) * u.g
if M_a is None:
M_a = 0.0 * u.g
if self.solid_frac != 0.0:
- M_a = self.solid_frac / (1.0 - self.solid_frac) * M_c
+ M_a = self.solid_frac / (1.0 - self.solid_frac) * M_c
self.M_c = 1.0 * M_c.cgs.value # core mass
self.M_a = 1.0 * M_a.cgs.value # atmospheric mass
self.M = self.M_a + self.M_c
- self.accretion_models = accretion_models # list containing the accretion objects
+ self.accretion_models = (
+ accretion_models # list containing the accretion objects
+ )
self._accretion_models_dict = {acc.name: acc for acc in self.accretion_models}
self._accrete = eval_kwargs(kwargs.get("accrete", True))
self.matter_removal = eval_kwargs(kwargs.get("matter_removal", True))
self.force_insitu = eval_kwargs(kwargs.get("force_insitu", False))
- self.r_in = eval_kwargs(kwargs.get("r_in", self.disk.r_i[0]*u.cm).cgs.value) # where to stop the planet formation
+ self.r_in = eval_kwargs(
+ kwargs.get("r_in", self.disk.r_i[0] * u.cm).cgs.value
+ ) # where to stop the planet formation
self.name = output.get("name", "planet")
- self.output_file = output.get("directory", os.path.join(OUTPUT_PATH, "{}.h5").format(self.name))
+ self.output_file = output.get(
+ "directory", os.path.join(OUTPUT_PATH, "{}.h5").format(self.name)
+ )
self.tau_m = np.inf # migration rate
- self.fSigma = 1.0 # gap depth
+ self.fSigma = 1.0 # gap depth
- self.M_end = eval_kwargs(kwargs.get("M_end", 1.0 * u.M_sun).cgs.value) # final mass at which we interrupt the formation of the planet
+ self.M_end = eval_kwargs(
+ kwargs.get("M_end", 1.0 * u.M_sun).cgs.value
+ ) # final mass at which we interrupt the formation of the planet
- self.past_pebble_iso = False # boolean flag that stores if we passed the pebiso mass
- self._keep_peb_iso = eval_kwargs(kwargs.get("keep_peb_iso", True)) # keep pebble isolation mass constant when pebble isolation mass is reached (otherwise we might have multiple phases of pebble accretion)
- self._use_pebiso_diffusion = eval_kwargs(kwargs.get("use_pebiso_diffusion", True)) # use the diffusion part of the pebble isolation mass
+ self.past_pebble_iso = (
+ False # boolean flag that stores if we passed the pebiso mass
+ )
+ self._keep_peb_iso = eval_kwargs(
+ kwargs.get("keep_peb_iso", True)
+ ) # keep pebble isolation mass constant when pebble isolation mass is reached (otherwise we might have multiple phases of pebble accretion)
+ self._use_pebiso_diffusion = eval_kwargs(
+ kwargs.get("use_pebiso_diffusion", True)
+ ) # use the diffusion part of the pebble isolation mass
self._init_accretion() # initialise the accretion models
self.update(second=True) # first update of all quantities
- self.chem_mask_array = self.disk._chem.mask_array != 0 # see mask_array in chemistry class (needed since there are less elements than mol. species)
- self.comp_a = self.chemistry_solid * self.M_a # molecular and elemental composition of the planetary envelope
- self.comp_c = self.chemistry_solid * self.M_c # molecular and elemental composition of the planetary core
-
- self.output = DataObject(self, output) # Initialize the IO class
- self.evolve_disk_init() # start evolving the disk to the point where the planet is placed into the disk
+ self.chem_mask_array = (
+ self.disk._chem.mask_array != 0
+ ) # see mask_array in chemistry class (needed since there are less elements than mol. species)
+ self.comp_a = (
+ self.chemistry_solid * self.M_a
+ ) # molecular and elemental composition of the planetary envelope
+ self.comp_c = (
+ self.chemistry_solid * self.M_c
+ ) # molecular and elemental composition of the planetary core
+
+ self.output = DataObject(self, output) # Initialize the IO class
+ self.evolve_disk_init() # start evolving the disk to the point where the planet is placed into the disk
if pebiso_start: # if planet is seeded with M=M_iso
self.calc_pebble_iso()
- print("{}: starting at peb_iso with {:.2e} M_e".format(self.name,(self.peb_iso*u.g).to(u.earthMass)))
- self.M_c = (1.0-self.solid_frac) * self.peb_iso
+ print(
+ "{}: starting at peb_iso with {:.2e} M_e".format(
+ self.name, (self.peb_iso * u.g).to(u.earthMass)
+ )
+ )
+ self.M_c = (1.0 - self.solid_frac) * self.peb_iso
self.M_a = self.solid_frac * self.peb_iso
self.M = self.M_a + self.M_c
self.past_pebble_iso = True
@@ -216,18 +337,17 @@ def __init__(
self.comp_a = self.chemistry_solid * self.M_a
self.comp_c = self.chemistry_solid * self.M_c
-
def _init_accretion(self):
- """ Initialise the Accretionmodels. """
+ """Initialise the Accretionmodels."""
for acc in self.accretion_models:
acc.init_accretion(self)
def update_torque(self):
- """ Calculate the Torques"""
+ """Calculate the Torques"""
pass
def update_a_p(self):
- """ Migrate the Planet. """
+ """Migrate the Planet."""
pass
def update(self, second, interpolation_idx=None, interpolation_idx_interface=None):
@@ -239,32 +359,57 @@ def update(self, second, interpolation_idx=None, interpolation_idx_interface=Non
interpolation_idx can be used if gap should be skipped during the calculation of sigma_gap relevant quantities
"""
- if interpolation_idx is None: # indexes as which we want to interpolate the disk quantities to the planetary position
- interpolation_idx = np.ones_like(self.disk.r,dtype=bool)
+ if (
+ interpolation_idx is None
+ ): # indexes as which we want to interpolate the disk quantities to the planetary position
+ interpolation_idx = np.ones_like(self.disk.r, dtype=bool)
if interpolation_idx_interface is None:
- interpolation_idx_interface = np.ones_like(self.disk.r_i,dtype=bool)
-
- self.sigma_g = ( # surface density of the gas phase
- interp1d(self.disk.r[interpolation_idx], self.disk.sigma_g[interpolation_idx], fill_value="extrapolate", assume_sorted=True)(self.a_p)
- )
- self.chemistry_gas = interp1d(self.disk.r, self.disk.chemistry_gas, axis=0, fill_value="extrapolate",
- assume_sorted=True)(self.a_p) # composition of the gas phase
+ interpolation_idx_interface = np.ones_like(self.disk.r_i, dtype=bool)
+
+ self.sigma_g = interp1d( # surface density of the gas phase
+ self.disk.r[interpolation_idx],
+ self.disk.sigma_g[interpolation_idx],
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
+ self.chemistry_gas = interp1d(
+ self.disk.r,
+ self.disk.chemistry_gas,
+ axis=0,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(
+ self.a_p
+ ) # composition of the gas phase
if self.t < self.disk.begin_photevap:
assert np.all(self.chemistry_gas >= 0.0), "gas chemistry is negative!"
elif not np.all(self.chemistry_gas >= 0.0):
- self.chemistry_gas = np.maximum(self.chemistry_gas, np.zeros_like(self.chemistry_gas)) # composition of the gas phase
+ self.chemistry_gas = np.maximum(
+ self.chemistry_gas, np.zeros_like(self.chemistry_gas)
+ ) # composition of the gas phase
if self.t < self.disk.begin_photevap:
assert self.sigma_g >= 0.0, "negative surface density!"
else:
- self.sigma_g = np.maximum(self.sigma_g, 1e-60) # surface density
+ self.sigma_g = np.maximum(self.sigma_g, 1e-60) # surface density
- self.sigma_g_components = self.sigma_g * self.chemistry_gas # surface density as a compositional vector
+ self.sigma_g_components = (
+ self.sigma_g * self.chemistry_gas
+ ) # surface density as a compositional vector
- self.nu = interp1d(self.disk.r, self.disk.nu, fill_value="extrapolate", assume_sorted=True)(self.a_p) # viscosity
+ self.nu = interp1d(
+ self.disk.r, self.disk.nu, fill_value="extrapolate", assume_sorted=True
+ )(
+ self.a_p
+ ) # viscosity
self.m_dot_disk = ( # viscous accretion rate of the gas phase of the disk
- interp1d(self.disk.r_i[interpolation_idx_interface], np.abs(self.disk.m_dot[interpolation_idx_interface]),fill_value="extrapolate", assume_sorted=True)(self.a_p)
+ interp1d(
+ self.disk.r_i[interpolation_idx_interface],
+ np.abs(self.disk.m_dot[interpolation_idx_interface]),
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
)
if self.m_dot_disk <= 0.0:
@@ -274,55 +419,129 @@ def update(self, second, interpolation_idx=None, interpolation_idx_interface=Non
self.m_dot_disk_components = self.chemistry_gas * self.m_dot_disk
if self.t < self.disk.begin_photevap:
- assert self.m_dot_disk >= 0.0, "the absolute value of m_dot_disk is not positive! Interpolation error"
+ assert (
+ self.m_dot_disk >= 0.0
+ ), "the absolute value of m_dot_disk is not positive! Interpolation error"
else:
- self.m_dot_disk_components = np.maximum(self.m_dot_disk_components, np.ones_like(self.chemistry_gas)*1e-60)
+ self.m_dot_disk_components = np.maximum(
+ self.m_dot_disk_components, np.ones_like(self.chemistry_gas) * 1e-60
+ )
self.m_dot_disk = np.maximum(1e-60, self.m_dot_disk)
- self.sigma_g_pert = ( # perturbed surface density of the gap
- interp1d(self.disk.r, self.disk.sigma_g, fill_value="extrapolate", assume_sorted=True)(self.a_p)
- )
- self.sigma_g_components_pert = interp1d(self.disk.r, # perturbed surface density of the gap as a compositional vector
- self.disk.sigma_g_components, axis=0,
- fill_value="extrapolate", assume_sorted=True)(
- self.a_p)
- self.T = interp1d(self.disk.r, self.disk.T, fill_value="extrapolate", assume_sorted=True)(self.a_p) # Temperature
- self.H_g = interp1d(self.disk.r, self.disk.aspect_ratio, fill_value="extrapolate", assume_sorted=True)(self.a_p) * self.a_p # Disk scalehight
- self.c_s = interp1d(self.disk.r, self.disk.c_s, fill_value="extrapolate", assume_sorted=True)(self.a_p) # soundspeed
- self.rho_g = self.sigma_g / (np.sqrt(2.0 * np.pi) * self.H_g) # gas volumedensity
- self.grad_sig = interp1d(self.disk.r[interpolation_idx], self.disk.grad_sig[interpolation_idx], fill_value="extrapolate", assume_sorted=True)(self.a_p) # gradient of the surfacedensity
- self.grad_T = interp1d(self.disk.r, self.disk.grad_T, fill_value="extrapolate", assume_sorted=True)(self.a_p) # gradient of the Temperature
- self.grad_p = interp1d(self.disk.r, self.disk.grad_p, fill_value="extrapolate", assume_sorted=True)(self.a_p) # gradient of the pressure
- self.r_h = self.a_p * (self.M / (3.0 * self.disk.M_s)) ** (1.0 / 3.0) # hill radius
+ self.sigma_g_pert = interp1d( # perturbed surface density of the gap
+ self.disk.r, self.disk.sigma_g, fill_value="extrapolate", assume_sorted=True
+ )(self.a_p)
+ self.sigma_g_components_pert = interp1d(
+ self.disk.r, # perturbed surface density of the gap as a compositional vector
+ self.disk.sigma_g_components,
+ axis=0,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
+ self.T = interp1d(
+ self.disk.r, self.disk.T, fill_value="extrapolate", assume_sorted=True
+ )(
+ self.a_p
+ ) # Temperature
+ self.H_g = (
+ interp1d(
+ self.disk.r,
+ self.disk.aspect_ratio,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
+ * self.a_p
+ ) # Disk scalehight
+ self.c_s = interp1d(
+ self.disk.r, self.disk.c_s, fill_value="extrapolate", assume_sorted=True
+ )(
+ self.a_p
+ ) # soundspeed
+ self.rho_g = self.sigma_g / (
+ np.sqrt(2.0 * np.pi) * self.H_g
+ ) # gas volumedensity
+ self.grad_sig = interp1d(
+ self.disk.r[interpolation_idx],
+ self.disk.grad_sig[interpolation_idx],
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(
+ self.a_p
+ ) # gradient of the surfacedensity
+ self.grad_T = interp1d(
+ self.disk.r, self.disk.grad_T, fill_value="extrapolate", assume_sorted=True
+ )(
+ self.a_p
+ ) # gradient of the Temperature
+ self.grad_p = interp1d(
+ self.disk.r, self.disk.grad_p, fill_value="extrapolate", assume_sorted=True
+ )(
+ self.a_p
+ ) # gradient of the pressure
+ self.r_h = self.a_p * (self.M / (3.0 * self.disk.M_s)) ** (
+ 1.0 / 3.0
+ ) # hill radius
self.v_k = np.sqrt(G * self.disk.M_s / self.a_p) # keppler velocity
self.omega_k = self.v_k / self.a_p # keppler orbital velocity
- self.del_v = interp1d(self.disk.r, self.disk.eta, fill_value="extrapolate", assume_sorted=True)(
- self.a_p) * self.v_k # difference between keppler and real gas orbital speed
- self.r_b = (G * self.M) / (self.del_v ** 2.0) # bondi radius
+ self.del_v = (
+ interp1d(
+ self.disk.r, self.disk.eta, fill_value="extrapolate", assume_sorted=True
+ )(self.a_p)
+ * self.v_k
+ ) # difference between keppler and real gas orbital speed
+ self.r_b = (G * self.M) / (self.del_v**2.0) # bondi radius
self.v_h = self.r_h * self.omega_k # hill velocity
- self.P = interp1d(self.disk.r, self.disk.P, fill_value="extrapolate", assume_sorted=True)(
- self.a_p)
+ self.P = interp1d(
+ self.disk.r, self.disk.P, fill_value="extrapolate", assume_sorted=True
+ )(self.a_p)
if second:
# This may only be called after migration has been applied
# migration calculates the gap parameter according to Kanagawa/Ndugu -> self.fsigma
- self.pebble_flux = interp1d(self.disk.r, self.disk.pebble_flux, fill_value="extrapolate", assume_sorted=True)(self.a_p) # pebble flux
+ self.pebble_flux = interp1d(
+ self.disk.r,
+ self.disk.pebble_flux,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(
+ self.a_p
+ ) # pebble flux
if hasattr(self.disk, "f_m"):
- self.sigma_peb = interp1d(self.disk.r, self.disk.sigma_dust * self.disk.f_m, fill_value="extrapolate", assume_sorted=True)(self.a_p) # pebble surface density
- self.sigma_peb_components = interp1d(self.disk.r, # pebble surface density as a compositional vector
- self.disk.sigma_dust_components * self.disk.f_m[:, np.newaxis,
- np.newaxis], axis=0,
- fill_value="extrapolate",
- assume_sorted=True)(self.a_p)
+ self.sigma_peb = interp1d(
+ self.disk.r,
+ self.disk.sigma_dust * self.disk.f_m,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(
+ self.a_p
+ ) # pebble surface density
+ self.sigma_peb_components = interp1d(
+ self.disk.r, # pebble surface density as a compositional vector
+ self.disk.sigma_dust_components
+ * self.disk.f_m[:, np.newaxis, np.newaxis],
+ axis=0,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
else:
- self.sigma_peb = interp1d(self.disk.r, self.disk.sigma_dust, fill_value="extrapolate", assume_sorted=True)(self.a_p) # checked
- self.sigma_peb_components = interp1d(self.disk.r,
- self.disk.sigma_dust_components, axis=0,
- fill_value="extrapolate",
- assume_sorted=True)(self.a_p)
+ self.sigma_peb = interp1d(
+ self.disk.r,
+ self.disk.sigma_dust,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(
+ self.a_p
+ ) # checked
+ self.sigma_peb_components = interp1d(
+ self.disk.r,
+ self.disk.sigma_dust_components,
+ axis=0,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
self.calc_pebble_iso()
@@ -333,18 +552,30 @@ def update(self, second, interpolation_idx=None, interpolation_idx_interface=Non
# idx: cell center
# tested -> Check
- self.idx_i_r = np.searchsorted(self.disk.r_i, self.a_p,
- side="right") # interface i+1: idx of right cell interface in self.disk.r_i
- self.idx_i_l = max(self.idx_i_r - 1, 0) # interface i: idx of left cell interface in self.disk.r_i
- self.idx = 1 * self.idx_i_l # center i: idx of cell (self.disk.r) planet is located in
-
- self.chemistry_solid = interp1d(self.disk.r, self.disk.chemistry_solid, axis=0, fill_value="extrapolate",
- assume_sorted=True)(self.a_p) # compositional vector of pebbles
+ self.idx_i_r = np.searchsorted(
+ self.disk.r_i, self.a_p, side="right"
+ ) # interface i+1: idx of right cell interface in self.disk.r_i
+ self.idx_i_l = max(
+ self.idx_i_r - 1, 0
+ ) # interface i: idx of left cell interface in self.disk.r_i
+ self.idx = (
+ 1 * self.idx_i_l
+ ) # center i: idx of cell (self.disk.r) planet is located in
+
+ self.chemistry_solid = interp1d(
+ self.disk.r,
+ self.disk.chemistry_solid,
+ axis=0,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(
+ self.a_p
+ ) # compositional vector of pebbles
self.update_parameters()
def update_parameters(self):
- """ Parameters that are dependent on planet Mass need to be kept updated """
+ """Parameters that are dependent on planet Mass need to be kept updated"""
pass
@property
@@ -352,32 +583,44 @@ def r_c(self):
"""
radius of the core
"""
- return (3.0 / 4.0 * self.M_c / (np.pi * self.rho_c)) ** (1.0 / 3.0) # core radius
+ return (3.0 / 4.0 * self.M_c / (np.pi * self.rho_c)) ** (
+ 1.0 / 3.0
+ ) # core radius
def calc_pebble_iso(self):
- """ calculate the pebble isolation mass."""
+ """calculate the pebble isolation mass."""
if self._keep_peb_iso and self.past_pebble_iso:
return
self.disk.calc_pebble_iso(diffusion=self._use_pebiso_diffusion)
- if hasattr(self.disk,"interpolation_idx"):
- self.peb_iso = interp1d(self.disk.r[self.disk.interpolation_idx], self.disk.peb_iso[self.disk.interpolation_idx], fill_value="extrapolate", assume_sorted=True)(self.a_p)
+ if hasattr(self.disk, "interpolation_idx"):
+ self.peb_iso = interp1d(
+ self.disk.r[self.disk.interpolation_idx],
+ self.disk.peb_iso[self.disk.interpolation_idx],
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
else:
- self.peb_iso = interp1d(self.disk.r, self.disk.peb_iso, fill_value="extrapolate", assume_sorted=True)(self.a_p)
+ self.peb_iso = interp1d(
+ self.disk.r,
+ self.disk.peb_iso,
+ fill_value="extrapolate",
+ assume_sorted=True,
+ )(self.a_p)
self.past_pebble_iso = True if self.M > self.peb_iso else False
return
def update_all(self):
- """ Wrapper that handles the complete update of disk and planet as well as migration."""
+ """Wrapper that handles the complete update of disk and planet as well as migration."""
self.disk.update(self.t, self.h)
self.update(second=False)
self.update_a_p()
self.update(second=True) # update all quantities effected by a_p
def evolve_disk_init(self):
- """ Evolve the disk before the protoplanetary seed is placed into the disk. """
+ """Evolve the disk before the protoplanetary seed is placed into the disk."""
self.disk.evolve_init()
# initialise disk and save
@@ -394,7 +637,7 @@ def evolve_disk_init(self):
self.print_params()
def calc_timestep(self):
- """ Fix the timestep and update time. Important: This is a fixed timestep currently!"""
+ """Fix the timestep and update time. Important: This is a fixed timestep currently!"""
if self.disk._timestep_input:
self.h = self.disk._timestep_input
else:
@@ -416,14 +659,14 @@ def keep_running(self):
return True
def grow_mass(self) -> None:
- """ Main loop and supervisor of chemcomp.
+ """Main loop and supervisor of chemcomp.
dt: the total time that the model should be evolved.
Needs to be a astropy Quantity
"""
def increment_mass(acc):
- """ increment the mass of the planet for a certain accretion model """
+ """increment the mass of the planet for a certain accretion model"""
acc.calc_m_dot()
acc.update_z()
self.total_m_dot += acc.m_dot
@@ -435,7 +678,7 @@ def increment_mass(acc):
acc.remove_mass_from_flux()
def _grow():
- """ wrap around increment_mass and update final mass. """
+ """wrap around increment_mass and update final mass."""
self.total_m_dot = 0.0
list(map(increment_mass, self.accretion_models))
self.M = self.M_a + self.M_c
@@ -481,7 +724,7 @@ def _grow():
return
def print_params(self, force=False):
- """" Handle the IO using self.output -> DataObject.
+ """ " Handle the IO using self.output -> DataObject.
Saves !after! accretion applied
"""
self.output.save(force)
@@ -492,8 +735,14 @@ def print_params(self, force=False):
CO = self.comp_a[0, 0] / self.comp_a[0, 1] * 16.0 / 12.0
except FloatingPointError:
CO = 0.0
- print("{}: t={:.1f}Myr, dt={:.1f}yr, M={:.1f}M_e, r={:.1f}AU, CO={:.2f}".format(self.name, self.t / Myr,
- self.h / year,
- self.M / M_earth,
- self.a_p / AU, CO))
+ print(
+ "{}: t={:.1f}Myr, dt={:.1f}yr, M={:.1f}M_e, r={:.1f}AU, CO={:.2f}".format(
+ self.name,
+ self.t / Myr,
+ self.h / year,
+ self.M / M_earth,
+ self.a_p / AU,
+ CO,
+ )
+ )
self._print_number += 1
diff --git a/chemcomp/planets/bertmigration.py b/chemcomp/planets/bertmigration.py
index db868b8..d1b7f6f 100755
--- a/chemcomp/planets/bertmigration.py
+++ b/chemcomp/planets/bertmigration.py
@@ -13,24 +13,32 @@ class BertPlanet(Planet):
"""Migrating planet. torques adapted from Planettrack.F by Bertram Bitsch. Planet used in Schneider & Bitsch (2021a,b)"""
def __init__(
- self,
- output,
- disk,
- accretion_models=[],
- rho_c=5.5 * u.g / (u.cm ** 3),
- a_p=5.2 * u.au,
- t_0=5e5 * u.yr,
- M_c=None,
- M_a=None,
- **kwargs
+ self,
+ output,
+ disk,
+ accretion_models=[],
+ rho_c=5.5 * u.g / (u.cm**3),
+ a_p=5.2 * u.au,
+ t_0=5e5 * u.yr,
+ M_c=None,
+ M_a=None,
+ **kwargs
):
self.gamma_tot = 0.0 # total torque
self.gamma_norm = 0.0 # torque normalized to gamma_0
- self._use_heat_torque = eval_kwargs(kwargs.get("use_heat_torque", True)) # use masset heating torque
- self._use_dynamical_torque = eval_kwargs(kwargs.get("use_dynamical_torque", True)) # use paardekooper dynamical torque
- self._apply_gap = eval_kwargs(kwargs.get("apply_gap_profile", True)) # if the gap should be added to the gas surface density
- self._migrate = eval_kwargs(kwargs.get("migration", True)) # if the planet should migrate
+ self._use_heat_torque = eval_kwargs(
+ kwargs.get("use_heat_torque", True)
+ ) # use masset heating torque
+ self._use_dynamical_torque = eval_kwargs(
+ kwargs.get("use_dynamical_torque", True)
+ ) # use paardekooper dynamical torque
+ self._apply_gap = eval_kwargs(
+ kwargs.get("apply_gap_profile", True)
+ ) # if the gap should be added to the gas surface density
+ self._migrate = eval_kwargs(
+ kwargs.get("migration", True)
+ ) # if the planet should migrate
super().__init__(
a_p=a_p,
@@ -46,11 +54,11 @@ def __init__(
self._old_a_p = 1.0 * self.a_p
self.init_paardekooper_units()
- print(self.name, (self.M*u.g).to(u.earthMass).value)
+ print(self.name, (self.M * u.g).to(u.earthMass).value)
def update_all(self):
- """ Update all quantities including disk evolution, migration and planetary quantities.
- Overwrites parent update_all() method """
+ """Update all quantities including disk evolution, migration and planetary quantities.
+ Overwrites parent update_all() method"""
if self._apply_gap:
self.update(second=False)
self.calc_gammaeff(migphase=False)
@@ -59,25 +67,43 @@ def update_all(self):
if self.fSigma < 1 and self.past_pebble_iso:
self.disk.apply_gap_profile(self.a_p, self.fSigma, gap_width)
- interpolation_idx = self.disk.interpolation_idx if hasattr(self.disk,"interpolation_idx") and self._apply_gap else None
- interpolation_idx_interface = self.disk.interpolation_idx_interface if hasattr(self.disk, "interpolation_idx_interface") and self._apply_gap else None
+ interpolation_idx = (
+ self.disk.interpolation_idx
+ if hasattr(self.disk, "interpolation_idx") and self._apply_gap
+ else None
+ )
+ interpolation_idx_interface = (
+ self.disk.interpolation_idx_interface
+ if hasattr(self.disk, "interpolation_idx_interface") and self._apply_gap
+ else None
+ )
self.disk.update(self.t, self.h)
- self.update(second=False, interpolation_idx=interpolation_idx, interpolation_idx_interface=interpolation_idx_interface)
+ self.update(
+ second=False,
+ interpolation_idx=interpolation_idx,
+ interpolation_idx_interface=interpolation_idx_interface,
+ )
self.update_a_p()
- self.update(second=True, interpolation_idx=interpolation_idx, interpolation_idx_interface=interpolation_idx_interface) # update all quantities effected by a_p
+ self.update(
+ second=True,
+ interpolation_idx=interpolation_idx,
+ interpolation_idx_interface=interpolation_idx_interface,
+ ) # update all quantities effected by a_p
def init_paardekooper_units(self):
- """ Define the paardekooper scaling factors."""
+ """Define the paardekooper scaling factors."""
self._r_0 = 1.0 * self.a_p
self._omega_0 = 1.0 * self.omega_k
self._sigma_0 = 1.0 * self.sigma_g
- self._q_d = np.pi * self._r_0 ** 2.0 * self._sigma_0 / self.disk.M_s
+ self._q_d = np.pi * self._r_0**2.0 * self._sigma_0 / self.disk.M_s
self._nu_0 = 1.0 * self.nu
# altough mig beta and alpha will change during runtime, we keep them constant since a gap would manipulate this even more!
self.mig_beta = -self.grad_T # checked
self.mig_alpha = -self.grad_sig # checked
- self.xi = self.mig_beta - (self.disk._chem.gamma - 1.0) * self.mig_alpha # checked, see equation 5
+ self.xi = (
+ self.mig_beta - (self.disk._chem.gamma - 1.0) * self.mig_alpha
+ ) # checked, see equation 5
def calc_heating_torque(self):
"""
@@ -89,13 +115,24 @@ def calc_heating_torque(self):
if self.M < thermal_crit:
L_c_mod_M = 4.0 * np.pi * G * self.Vchi * self.rho_g / self._gammaeff
- L_mod_M = G * (self.accretion_models[0].m_dot + self.accretion_models[1].m_dot) / self.r_c
+ L_mod_M = (
+ G
+ * (self.accretion_models[0].m_dot + self.accretion_models[1].m_dot)
+ / self.r_c
+ )
L_norm = L_mod_M / L_c_mod_M
lambda_c = np.sqrt(self.Vchi / (1.5 * self.omega_k * self._gammaeff))
eta = self.mig_alpha / 3.0 + (self.mig_beta + 3.0) / 6.0
- gamma_thermal = 1.61 * (self._gammaeff - 1.0) / self._gammaeff * eta * (self.H_g / lambda_c) * (L_norm - 1.0)
+ gamma_thermal = (
+ 1.61
+ * (self._gammaeff - 1.0)
+ / self._gammaeff
+ * eta
+ * (self.H_g / lambda_c)
+ * (L_norm - 1.0)
+ )
self.gamma_tot += self.gamma_zero * gamma_thermal
self.gamma_norm += gamma_thermal
@@ -123,9 +160,15 @@ def update_a_p(self):
self.tau_m = self.get_static_tau_m()
if self.fSigma <= 0.53:
- tau_II = - self.a_p ** 2.0 / self.nu * np.maximum(1.0, self.M/(4*np.pi*self.sigma_g*self.a_p**2.0))
+ tau_II = (
+ -self.a_p**2.0
+ / self.nu
+ * np.maximum(1.0, self.M / (4 * np.pi * self.sigma_g * self.a_p**2.0))
+ )
if self.fSigma > 0.1:
- self.tau_m = (self.tau_m - tau_II) / (0.53 - 0.1) * (self.fSigma-0.1) + tau_II
+ self.tau_m = (self.tau_m - tau_II) / (0.53 - 0.1) * (
+ self.fSigma - 0.1
+ ) + tau_II
else:
self.tau_m = tau_II
@@ -145,12 +188,23 @@ def get_dynamical_tau_m(self):
"""
- k = 8.0 / 3.0 / np.pi * (3.0 / 2.0 - self.mig_alpha) * self.gamma_norm * self._q_d ** 2.0 * self._xs ** 3.0 / (
- self.aspect_ratio ** 2.0) * self._r_0 ** 2.0 * self._omega_0 / self._nu_0 * (
- self.a_p / self._r_0) ** (5.0 - 3.0 * self.mig_alpha)
+ k = (
+ 8.0
+ / 3.0
+ / np.pi
+ * (3.0 / 2.0 - self.mig_alpha)
+ * self.gamma_norm
+ * self._q_d**2.0
+ * self._xs**3.0
+ / (self.aspect_ratio**2.0)
+ * self._r_0**2.0
+ * self._omega_0
+ / self._nu_0
+ * (self.a_p / self._r_0) ** (5.0 - 3.0 * self.mig_alpha)
+ )
k = np.minimum(k, 0.5)
theta = (1.0 - np.sqrt(1.0 - 2.0 * k)) / k
- J_p = self.M * self.a_p ** 2.0 * self.omega_k
+ J_p = self.M * self.a_p**2.0 * self.omega_k
return 0.5 * J_p / self.gamma_tot / theta
def get_static_tau_m(self):
@@ -162,81 +216,93 @@ def get_static_tau_m(self):
tau_m: Migration timescale [s]
"""
- J_p = self.M * self.a_p ** 2.0 * self.omega_k
+ J_p = self.M * self.a_p**2.0 * self.omega_k
return 0.5 * J_p / self.gamma_tot
@property
def _kappa(self):
"""Interpolate opacity to planetary position. Adds minimal offset (for stability)"""
- return (
- np.interp(self.a_p, self.disk.r, self.disk.opacity)
- + 1.0e-300
- )
+ return np.interp(self.a_p, self.disk.r, self.disk.opacity) + 1.0e-300
def update_torque(self):
"""Calculate Paardekooper2011 torques"""
self.calc_gammaeff(migphase=True)
- ang_p = self.a_p ** 2.0 * self.omega_k # spec. Ang.mom of planet, checked
- pnu = 2.0 / 3.0 * np.sqrt((ang_p * self._xs ** 3.0) / (2.0 * np.pi * self.nu)) # checked
- pxi = np.sqrt((ang_p * self._xs ** 3.0) / (2.0 * np.pi * self.Vchi)) # checked
+ ang_p = self.a_p**2.0 * self.omega_k # spec. Ang.mom of planet, checked
+ pnu = (
+ 2.0 / 3.0 * np.sqrt((ang_p * self._xs**3.0) / (2.0 * np.pi * self.nu))
+ ) # checked
+ pxi = np.sqrt((ang_p * self._xs**3.0) / (2.0 * np.pi * self.Vchi)) # checked
Fpnu = 1.0 / (1.0 + (pnu / 1.3) ** 2.0) # checked
Fpxi = 1.0 / (1.0 + (pxi / 1.3) ** 2.0) # checked
if pnu < np.sqrt(8.0 / (45.0 * np.pi)):
- Gpnu = 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pnu ** (1.5) # checked
+ Gpnu = (
+ 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pnu ** (1.5)
+ ) # checked
else:
Gpnu = 1.0 - 9.0 / 25.0 * (8.0 / (45.0 * np.pi)) ** (1.33333) * pnu ** (
- -8.0 / 3.0
+ -8.0 / 3.0
) # checked
if (pxi) < (np.sqrt(8.0 / (45.0 * np.pi))):
- Gpxi = 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pxi ** (1.5) # checked
+ Gpxi = (
+ 16.0 / 25.0 * (45.0 * np.pi / 8.0) ** (0.75) * pxi ** (1.5)
+ ) # checked
else:
Gpxi = 1.0 - 9.0 / 25.0 * (8.0 / (45.0 * np.pi)) ** (1.33333) * pxi ** (
- -8.0 / 3.0
+ -8.0 / 3.0
) # checked
if (pnu) < (np.sqrt(28.0 / (45.0 * np.pi))):
- Kpnu = 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** 0.75 * pnu ** 1.5 # checked
+ Kpnu = 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** 0.75 * pnu**1.5 # checked
else:
Kpnu = 1.0 - 9.0 / 25.0 * (28.0 / (45.0 * np.pi)) ** (1.33333) * pnu ** (
- -8.0 / 3.0
+ -8.0 / 3.0
) # checked
if (pxi) < (np.sqrt(28.0 / (45.0 * np.pi))):
- Kpxi = 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** (0.75) * pxi ** (1.5) # checked
+ Kpxi = (
+ 16.0 / 25.0 * (45.0 * np.pi / 28.0) ** (0.75) * pxi ** (1.5)
+ ) # checked
else:
Kpxi = 1.0 - 9.0 / 25.0 * (28.0 / (45.0 * np.pi)) ** (1.33333) * pxi ** (
- -8.0 / 3.0
+ -8.0 / 3.0
) # checked
# Calculate torque contributions
self.gamma_zero = (
- (self.massratio / self.aspect_ratio) ** 2.0
- * self.sigma_g
- * self.a_p ** 4.0
- * self.omega_k ** 2.0
+ (self.massratio / self.aspect_ratio) ** 2.0
+ * self.sigma_g
+ * self.a_p**4.0
+ * self.omega_k**2.0
) # checked
if self.fSigma < 1.0 and self.fSigma > 0.1:
self.gamma_zero = self.gamma_zero * self.fSigma
- self._m_gap = np.sqrt(self.aspect_ratio ** 5.0 * self.disk.alpha_mig * 1.0 / 0.04) * self.disk.M_s
+ self._m_gap = (
+ np.sqrt(self.aspect_ratio**5.0 * self.disk.alpha_mig * 1.0 / 0.04)
+ * self.disk.M_s
+ )
self._iso_vs_gap = self.peb_iso / self._m_gap
xfact = self.gamma_zero / self._gammaeff # checked
- GammaLIN = -(2.5 + 1.7 * self.mig_beta - 0.1 * self.mig_alpha) * xfact # checked
+ GammaLIN = (
+ -(2.5 + 1.7 * self.mig_beta - 0.1 * self.mig_alpha) * xfact
+ ) # checked
GammaHSENT = self.xi / self._gammaeff * 7.9 * xfact # equation 5, checked
- GammaCLINENT = (2.2 - 1.4 / self._gammaeff) * self.xi * xfact # equation 7, checked
+ GammaCLINENT = (
+ (2.2 - 1.4 / self._gammaeff) * self.xi * xfact
+ ) # equation 7, checked
GammaCLINBARO = 0.7 * (1.5 - self.mig_alpha) * xfact # equation 6, checked
GammaHSBARO = 1.1 * (1.5 - self.mig_alpha) * xfact # equation 4, checked
- GammaCBARO = (GammaHSBARO * Fpnu * Gpnu + (1.0 - Kpnu) * GammaCLINBARO) # checked
+ GammaCBARO = GammaHSBARO * Fpnu * Gpnu + (1.0 - Kpnu) * GammaCLINBARO # checked
GammaCENT = (
- GammaHSENT * Fpnu * Fpxi * np.sqrt(Gpnu * Gpxi)
- + np.sqrt((1.0 - Kpnu) * (1.0 - Kpxi)) * GammaCLINENT
+ GammaHSENT * Fpnu * Fpxi * np.sqrt(Gpnu * Gpxi)
+ + np.sqrt((1.0 - Kpnu) * (1.0 - Kpxi)) * GammaCLINENT
) # checked
# Normalization factor to Gamma_0 as defined in Paardekooper 2011
@@ -244,85 +310,95 @@ def update_torque(self):
self.gamma_norm = self.gamma_tot / self.gamma_zero # checked
def calc_gammaeff(self, migphase):
- """Calculate the effective gamma AND calculate the gapdepth """
+ """Calculate the effective gamma AND calculate the gapdepth"""
self.massratio = self.M / self.disk.M_s # checked
# rename gradients to self.mig_alpha and beta as in Paardekooper 2011
# Calculate Vchi
self.Vchi = ( # eq. 4 in Paardekooper 2011
- 16.0 # faktor of 4 above paardekooper
- * self.disk._chem.gamma
- * (self.disk._chem.gamma - 1.0)
- * sr
- * self.T ** 4.0
- / (
- 3.0
- * self._kappa
- * self.rho_g ** 2.0
- * self.H_g ** 2
- * self.omega_k ** 2
- )
+ 16.0 # faktor of 4 above paardekooper
+ * self.disk._chem.gamma
+ * (self.disk._chem.gamma - 1.0)
+ * sr
+ * self.T**4.0
+ / (
+ 3.0
+ * self._kappa
+ * self.rho_g**2.0
+ * self.H_g**2
+ * self.omega_k**2
+ )
) # checked
# calculate Q
Q = (
- 2.0
- * self.Vchi
- / (3.0 * (self.H_g / self.a_p) ** 3.0)
- / (self.a_p ** 2.0 * self.omega_k)
+ 2.0
+ * self.Vchi
+ / (3.0 * (self.H_g / self.a_p) ** 3.0)
+ / (self.a_p**2.0 * self.omega_k)
) # checked
self._gammaeff = (
- 2.0
- * Q
- * self.disk._chem.gamma
- / (
- self.disk._chem.gamma * Q
- + 0.5
- * np.sqrt(
+ 2.0
+ * Q
+ * self.disk._chem.gamma
+ / (
+ self.disk._chem.gamma * Q
+ + 0.5
+ * np.sqrt(
2.0
* np.sqrt(
- (self.disk._chem.gamma ** 2.0 * Q ** 2.0 + 1.0) ** 2.0 - 16.0 * Q ** 2.0 * (self.disk._chem.gamma - 1.0)
+ (self.disk._chem.gamma**2.0 * Q**2.0 + 1.0) ** 2.0
+ - 16.0 * Q**2.0 * (self.disk._chem.gamma - 1.0)
)
- + 2.0 * self.disk._chem.gamma ** 2.0 * Q ** 2.0
+ + 2.0 * self.disk._chem.gamma**2.0 * Q**2.0
- 2.0
)
- )
+ )
) # checked
# Calculate functions F,G,K
self.aspect_ratio = self.H_g / self.a_p # relative disk thickness
self._xs = (
- 1.11 / self._gammaeff ** 0.25 * np.sqrt(self.massratio / self.aspect_ratio)
+ 1.11 / self._gammaeff**0.25 * np.sqrt(self.massratio / self.aspect_ratio)
) # checked, dimensionless in units of r_p
#######
# GAP #
#######
if not hasattr(self, "_M_HS"):
- self._M_HS = 4 * np.pi * self._xs * self.a_p ** 2 * self.sigma_g
+ self._M_HS = 4 * np.pi * self._xs * self.a_p**2 * self.sigma_g
self._sigma_HS = 1 * self.sigma_g
- K = self.massratio ** 2.0 / (self.aspect_ratio ** 5.0 * self.disk.alpha_mig) # checked
+ K = self.massratio**2.0 / (
+ self.aspect_ratio**5.0 * self.disk.alpha_mig
+ ) # checked
self.fSigma_kanagawa = 1.0 / (1.0 + 0.04 * K) # checked
- R = self.a_p ** 2.0 * self.omega_k / self.nu
+ R = self.a_p**2.0 * self.omega_k / self.nu
Gap = 3.0 / 4.0 * self.H_g / self.r_h + 50.0 / (self.massratio * R)
if Gap <= 2.4646:
- f_P = (Gap - 0.541) / 4.
+ f_P = (Gap - 0.541) / 4.0
else:
- f_P = 1. - np.exp(-Gap ** 0.75 / 3.)
+ f_P = 1.0 - np.exp(-(Gap**0.75) / 3.0)
if self.past_pebble_iso:
_, m_dot_gas = self._accretion_models_dict["GasAccretion"].get_min()
if migphase:
- self._sigma_HS = self._sigma_HS * (self._old_a_p/self.a_p)**1.5
+ self._sigma_HS = self._sigma_HS * (self._old_a_p / self.a_p) ** 1.5
self._old_a_p = 1.0 * self.a_p
- M_HS_check = 4.0 * np.pi * self.a_p ** 2.0 * self._xs * self._sigma_HS
+ M_HS_check = 4.0 * np.pi * self.a_p**2.0 * self._xs * self._sigma_HS
if M_HS_check < self._M_HS:
self._M_HS = self._M_HS + (self.m_dot_disk - m_dot_gas) * self.h
else:
- self._M_HS = self._M_HS + (4.0 * np.pi * self.a_p ** 2.0 * self._xs - self._M_HS / self._sigma_HS) * self.sigma_g
+ self._M_HS = (
+ self._M_HS
+ + (
+ 4.0 * np.pi * self.a_p**2.0 * self._xs
+ - self._M_HS / self._sigma_HS
+ )
+ * self.sigma_g
+ )
- self._sigma_HS = self._M_HS / (4.0 * np.pi * self.a_p ** 2.0 * self._xs)
+ self._sigma_HS = self._M_HS / (4.0 * np.pi * self.a_p**2.0 * self._xs)
if self.t < self.disk.begin_photevap:
assert self._M_HS >= 0.0, "negative mass in horseshoe region"
@@ -335,4 +411,3 @@ def calc_gammaeff(self, migphase):
self.fSigma = np.maximum(f_A * f_P, 1e-7)
return
-
diff --git a/chemcomp/planets/migration_1.py b/chemcomp/planets/migration_1.py
index 9e71af4..d0de876 100755
--- a/chemcomp/planets/migration_1.py
+++ b/chemcomp/planets/migration_1.py
@@ -8,18 +8,17 @@ class LindbladPlanet(Planet):
"""simple migrating planet. torques calculated with viscosity of α = 0.004"""
def __init__(
- self,
- output,
- disk,
- accretion_models=[],
- rho_c=5.5 * u.g / (u.cm ** 3.0),
+ self,
+ output,
+ disk,
+ accretion_models=[],
+ rho_c=5.5 * u.g / (u.cm**3.0),
a_p=5.2 * u.au,
t_0=0.0 * u.yr,
M_c=None,
M_a=None,
**kwargs
):
-
super().__init__(
a_p=a_p,
M_c=M_c,
@@ -34,7 +33,7 @@ def __init__(
def update_a_p(self):
self.update_torque()
- J_p = self.M * self.a_p ** 2.0 * self.omega_k
+ J_p = self.M * self.a_p**2.0 * self.omega_k
self.tau_m = 0.5 * J_p / self._gamma_tot
self.a_p += self.a_p / self.tau_m * self.h
@@ -42,7 +41,7 @@ def update_torque(self):
adiabatic_index = 1.4
q = self.M / self.disk.M_s
h = self.H_g / self.a_p
- self._gamma_0 = (q / h) ** 2 * self.sigma_g * self.a_p ** 4 * self.omega_k ** 2
+ self._gamma_0 = (q / h) ** 2 * self.sigma_g * self.a_p**4 * self.omega_k**2
tgrad = np.interp(self.a_p, self.disk.r, self.disk.grad_T) # to be defined
Siggrad = np.interp(self.a_p, self.disk.r, self.disk.grad_sig) # checked
self._gamma_tot = (
diff --git a/chemcomp/planets/no_accretion.py b/chemcomp/planets/no_accretion.py
index 3044c87..4ff5f55 100755
--- a/chemcomp/planets/no_accretion.py
+++ b/chemcomp/planets/no_accretion.py
@@ -7,16 +7,16 @@ class NoAccretion(Planet):
"""NoAccretion planet. Useful if you are just interested in the disk and not in the planet itself"""
def __init__(
- self,
- output,
- disk,
- accretion_models=[],
- rho_c=5.5 * u.g / (u.cm ** 3),
- a_p=5.2 * u.au,
- t_0=0.0 * u.yr,
- M_c=None,
- M_a=None,
- **kwargs,
+ self,
+ output,
+ disk,
+ accretion_models=[],
+ rho_c=5.5 * u.g / (u.cm**3),
+ a_p=5.2 * u.au,
+ t_0=0.0 * u.yr,
+ M_c=None,
+ M_a=None,
+ **kwargs,
):
self.gamma_tot = 0
self.gamma_norm = 0
diff --git a/chemcomp/planets/simple_planet.py b/chemcomp/planets/simple_planet.py
index 74090a3..ec2a179 100755
--- a/chemcomp/planets/simple_planet.py
+++ b/chemcomp/planets/simple_planet.py
@@ -7,18 +7,17 @@ class SimplePlanet(Planet):
"""A very simple Planet without anything."""
def __init__(
- self,
- output,
- disk,
- accretion_models=[],
- rho_c=5.5 * u.g / (u.cm ** 3),
- a_p=5.2 * u.au,
- t_0=0.0 * u.yr,
- M_c=None,
- M_a=None,
- **kwargs,
+ self,
+ output,
+ disk,
+ accretion_models=[],
+ rho_c=5.5 * u.g / (u.cm**3),
+ a_p=5.2 * u.au,
+ t_0=0.0 * u.yr,
+ M_c=None,
+ M_a=None,
+ **kwargs,
):
-
super().__init__(
a_p=a_p,
M_c=M_c,
@@ -27,6 +26,5 @@ def __init__(
rho_c=rho_c,
disk=disk,
accretion_models=accretion_models,
- output=output
- ** kwargs,
+ output=output**kwargs,
)