black formating

AaronDavidSchneider · Nov 22, 2023 · 1f90761 · 1f90761
1 parent e3c0dc5
commit 1f90761
Show file tree

Hide file tree

Showing 22 changed files with 2,138 additions and 984 deletions.
diff --git 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
+from .pebbles import PebbleAccretion
diff --git 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
@@ -10,25 +10,22 @@
 
 
 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)
         self.f_disk_max = config.get("f_disk_max", 1)
 
     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
+        return
diff --git 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,42 +39,50 @@ 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
         else:
             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,21 +122,26 @@ 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
 
             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
@@ -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