Module diskchef.physics.base

Classes

class PhysicsBase (star_mass: Unit("solMass") = <Quantity 1. solMass>, xray_plasma_temperature: Unit("K") = <Quantity 10000000. K>, xray_luminosity: Unit("erg / s") = <Quantity 1.e+31 erg / s>, cr_padovani_use_l: bool = False)

Base class for disk physics models

Expand source code

class PhysicsBase:
    """Base class for disk physics models"""
    star_mass: u.solMass = 1 * u.solMass
    xray_plasma_temperature: u.K = 1e7 * u.K
    xray_luminosity: u.erg / u.s = 1e31 * u.erg / u.s
    cr_padovani_use_l: bool = False

    def __post_init__(self):
        self.logger = logging.getLogger(__name__ + '.' + self.__class__.__qualname__)
        self.logger.info("Creating an instance of %s", self.__class__.__qualname__)
        self.logger.debug("With parameters: %s", self.__dict__)

    @property
    def table(self) -> CTable:
        """Storage for all the position-dependent data of the disk"""
        return self._table

    @table.setter
    def table(self, value: CTable):
        self._table = value

    def plot_density(
            self,
            axes: matplotlib.axes.Axes = None,
            table: CTable = None, folder=".",
            cmap: Union[matplotlib.colors.Colormap, str] = 'PuBuGn',
            **kwargs
    ) -> Plot2D:
        """Plot 2D Gas and Dust density figures"""
        if table is None:
            table = self.table
        return Plot2D(table, axes=axes, data1="Gas density", data2="Dust density", cmap=cmap, **kwargs)

    def plot_temperatures(
            self,
            axes: matplotlib.axes.Axes = None,
            table: CTable = None, folder=".",
            cmap: Union[matplotlib.colors.Colormap, str] = 'afmhot',
            **kwargs
    ) -> Plot2D:
        """Plot 2D Gas and Dust temperature figures"""
        if table is None:
            self.check_temperatures()
            table = self.table
        return Plot2D(table, axes=axes, data1="Gas temperature", data2="Dust temperature", cmap=cmap, **kwargs)

    def plot_column_densities(
            self,
            axes: matplotlib.axes.Axes = None,
            table: CTable = None, folder=".",
            data=("Gas density", "Dust density"),
            **kwargs
    ) -> Plot1D:
        """Plot 1D column plots of Gas and Dust density. Use the code of this method as an example."""
        if table is None:
            table = self.table
        return Plot1D(table, axes=axes, data=data, **kwargs)

    def check_temperatures(self):
        self.table.check_zeros("Gas temperature")
        self.table.check_zeros("Dust temperature")

    @u.quantity_input
    def column_density_to(
            self, r: u.au, z: u.au,
            colname,
            r0: u.au = 0 * u.au, z0: u.au = 0 * u.au,
            steps=500, steps_for_log_fraction=0.9,
            only_gridpoint=False,
            point_by_point=False
    ):
        """
        Calculates column density of the given CTable quantity towards the `(r0, z0)`

        Args:
            r: radial coordinate of the desired point(s)
            z: vertical coordinate of the desired point(s)
            colname: name of self.table column to integrate
            r0: coordinate of the initial point (0 by default)
            z0: coordinate of the initial point (0 by default)
            steps: number of steps
            steps_for_log_fraction: fraction of steps to be taken in log scale
            only_gridpoint: when calculating column density towards star, whether to interpolate, or only take grid points
            point_by_point: force separate integration for each point if True
        Returns:
            column density, in self.table[colname].unit * u.cm
        """
        steps_for_log = int(steps_for_log_fraction * steps)
        steps_for_lin = steps - steps_for_log
        klog = np.geomspace(1 / steps_for_log, 1, steps_for_log)
        klin = np.linspace(0, 1 / steps_for_log, steps_for_lin, endpoint=False)
        k = np.concatenate([klin, klog])
        integrals = np.empty_like(self.table.r.value) << (u.cm * self.table[colname].unit)

        if (self.table.is_in_zr_regular_grid and r0 == z0 == 0 * u.au) and not point_by_point:
            zr_set = set(self.table.zr)
            for _zr in zr_set:
                this_zr_rows = np.where(self.table.zr == _zr)[0]
                int_r = r0 + (np.max(self.table.r[this_zr_rows]) - r0) * k
                int_z = z0 + (np.max(self.table.z[this_zr_rows]) - z0) * k
                if not only_gridpoint:
                    # Mix the grid elements into the int_r array
                    int_r = u.Quantity([*int_r, *self.table.r[this_zr_rows]])
                    int_z = u.Quantity([*int_z, *self.table.z[this_zr_rows]])
                    # The position of the grid elemenets in new int_r array
                    idx = np.where(int_r.argsort() - len(k) >= 0)[0]
                    int_r.sort()
                    int_z.sort()
                else:
                    idx = np.argsort(self.table.r[this_zr_rows])
                    int_r = self.table.r[this_zr_rows][idx]
                    int_z = self.table.z[this_zr_rows][idx]
                k = (int_r - r0) / (np.max(self.table.r[this_zr_rows]) - r0)
                k_length_cm = (
                        ((int_r[-1] - int_r[0]) ** 2 + (int_z[-1] - int_z[0]) ** 2) ** 0.5
                ).to(u.cm)
                value = self.table.interpolate(colname)(int_r, int_z)
                integrals_at_zr = scipy.integrate.cumtrapz(value.value, k, initial=0) * value.unit * k_length_cm
                integrals[this_zr_rows] = integrals_at_zr[idx]
        elif (self.table.is_in_zr_regular_grid and z0 == np.inf * u.au) and not point_by_point:
            radii = set(self.table.r)
            for _r in radii:
                indices = self.table.r == _r
                _z = self.table.z[indices]
                _data = self.table[colname][indices]
                order = np.argsort(_z)[::-1]
                integrals[indices] = np.abs(scipy.integrate.cumtrapz(
                    _data.value[order],
                    _z.value[order],
                    initial=0
                ))[np.argsort(order)] * _data.unit * _z.unit
        else:
            warnings.warn(CHEFSlowDownWarning(
                "Column density calculations are expensive, unless the grid is regular in z/r AND r0 == z0 == 0"
            ))
            for i, (_r, _z) in enumerate(zip(self.table.r, self.table.z)):
                int_r = r0 + (_r - r0) * k
                int_z = z0 + (_z - z0) * k
                k_length_cm = (
                        ((int_r[-1] - int_r[0]) ** 2 + (int_z[-1] - int_z[0]) ** 2) ** 0.5
                ).to(u.cm)
                value = self.table.interpolate(colname)(int_r, int_z)
                integral = scipy.integrate.trapz(value, k) * k_length_cm
                integrals[i] = integral
        return u.Quantity(integrals)

    def ionization(self, xray: str = "xray_bruderer", cosmic_ray: str = "cosmic_ray_padovani18"):
        """Calculate ionization by running x-ray and cosmic ray methods"""

        self.logger.debug("Calculating 'X ray ionization rate' according to %s", xray)
        getattr(self, xray)()

        self.logger.debug("Calculating 'CR ionization rate' according to %s", cosmic_ray)
        getattr(self, cosmic_ray)()

        self.table["Ionization rate"] = self.table["CR ionization rate"] + self.table["X ray ionization rate"]

    def xray_bruderer(self):
        """Calculate X-ray ionization rates using Bruderer+09 Table 3

        Requires `xray_plasma_temperature` and `xray_luminosity` to be set
        """
        if "Total density" not in self.table.colnames:
            self.table["Total density"] = self.table["Gas density"] + self.table["Dust density"]
        self.table["Nucleon column density towards star"] = (self.column_density_to(
            self.table.r, self.table.z,
            f"Total density",
            only_gridpoint=True,
        ) / u.u).to(u.cm ** -2)

        self.table["X ray ionization rate"] = (
                diskchef.physics.ionization.bruderer09(
                    self.table["Nucleon column density towards star"].to_value(u.cm ** -2),
                    self.xray_plasma_temperature.to_value(u.K)
                ) / u.s * (
                        self.xray_luminosity / (4 * np.pi * (self.table.r ** 2 + self.table.z ** 2))
                ).to_value(u.erg / u.s / u.cm ** 2)
        ).to(1 / u.s)

    def cosmic_ray_padovani18(self):
        """Calculate CR ionization rate according to App. F model L of Padovani+2018

        https://www.aanda.org/articles/aa/pdf/2018/06/aa32202-17.pdf

        TODO: what is about magnetic fields?
        """
        if "Total density" not in self.table.colnames:
            self.table["Total density"] = self.table["Gas density"] + self.table["Dust density"]

        if "Nucleon column density upwards" not in self.table.colnames:
            self.table["Nucleon column density upwards"] = (self.column_density_to(
                np.nan * u.au, np.nan * u.au,
                f"Total density",
                r0=np.nan * u.au, z0=np.inf * u.au,
                only_gridpoint=True,
            ) / u.u).to(u.cm ** -2)

        density_upwards = self.table["Nucleon column density upwards"].to_value(u.cm ** -2)

        midplane_i = self.table.zr == 0
        midplane_coldens = self.table["Nucleon column density upwards"][midplane_i]
        midplane_r = self.table.r[midplane_i]
        midplane_dict = dict(zip(midplane_r, midplane_coldens))
        coldens = u.Quantity([midplane_dict[r] for r in self.table.r]).to_value(u.cm ** -2)

        if self.cr_padovani_use_l:
            padovani = diskchef.physics.ionization.padovani18l
        else:
            padovani = diskchef.physics.ionization.padovani18h
        self.table["CR ionization rate"] = 0.5 * (
                padovani(np.log10(density_upwards)) +
                padovani(np.log10(2 * coldens - density_upwards))
        ) / u.s

Subclasses

• ReadAndesPhysData
• ReadHurakanData
• PhysicsModel
• Disk2D
• Disk3D

Class variables

var cr_padovani_use_l : bool

var star_mass : Unit("solMass")

var xray_luminosity : Unit("erg / s")

var xray_plasma_temperature : Unit("K")

Instance variables

prop table : CTable

Storage for all the position-dependent data of the disk

Expand source code

@property
def table(self) -> CTable:
    """Storage for all the position-dependent data of the disk"""
    return self._table

Methods

def check_temperatures(self)

def column_density_to(self, r: Unit("AU"), z: Unit("AU"), colname, r0: Unit("AU") = <Quantity 0. AU>, z0: Unit("AU") = <Quantity 0. AU>, steps=500, steps_for_log_fraction=0.9, only_gridpoint=False, point_by_point=False)

Calculates column density of the given CTable quantity towards the (r0, z0)

Args

r

radial coordinate of the desired point(s)

z

vertical coordinate of the desired point(s)

colname

name of self.table column to integrate

r0

coordinate of the initial point (0 by default)

z0

coordinate of the initial point (0 by default)

steps

number of steps

steps_for_log_fraction

fraction of steps to be taken in log scale

only_gridpoint

when calculating column density towards star, whether to interpolate, or only take grid points

point_by_point

force separate integration for each point if True

Returns

column density, in self.table[colname].unit * u.cm

def cosmic_ray_padovani18(self)

Calculate CR ionization rate according to App. F model L of Padovani+2018

https://www.aanda.org/articles/aa/pdf/2018/06/aa32202-17.pdf

TODO: what is about magnetic fields?

def ionization(self, xray: str = 'xray_bruderer', cosmic_ray: str = 'cosmic_ray_padovani18')

Calculate ionization by running x-ray and cosmic ray methods

def plot_column_densities(self, axes: matplotlib.axes._axes.Axes = None, table: CTable = None, folder='.', data=('Gas density', 'Dust density'), **kwargs) ‑> Plot1D

Plot 1D column plots of Gas and Dust density. Use the code of this method as an example.

def plot_density(self, axes: matplotlib.axes._axes.Axes = None, table: CTable = None, folder='.', cmap: Union[matplotlib.colors.Colormap, str] = 'PuBuGn', **kwargs) ‑> Plot2D

Plot 2D Gas and Dust density figures

def plot_temperatures(self, axes: matplotlib.axes._axes.Axes = None, table: CTable = None, folder='.', cmap: Union[matplotlib.colors.Colormap, str] = 'afmhot', **kwargs) ‑> Plot2D

Plot 2D Gas and Dust temperature figures

def xray_bruderer(self)

Calculate X-ray ionization rates using Bruderer+09 Table 3

Requires xray_plasma_temperature and xray_luminosity to be set

class PhysicsModel (star_mass: Unit("solMass") = <Quantity 1. solMass>, xray_plasma_temperature: Unit("K") = <Quantity 10000000. K>, xray_luminosity: Unit("erg / s") = <Quantity 1.e+31 erg / s>, cr_padovani_use_l: bool = False, r_min: Unit("AU") = <Quantity 0.1 AU>, r_max: Unit("AU") = <Quantity 500. AU>, zr_max: float = 0.7, radial_bins: int = 100, vertical_bins: int = 100, dust_to_gas: float = 0.01)

The base class describing the most basic parameters of the disk

Can not be used directly, rather subclasses should be used. See WilliamsBest2014 documentation for more details

Expand source code

class PhysicsModel(PhysicsBase):
    """
    The base class describing the most basic parameters of the disk

    Can not be used directly, rather subclasses should be used. See `WilliamsBest2014` documentation for more details
    """
    r_min: u.au = 0.1 * u.au
    r_max: u.au = 500 * u.au
    zr_max: float = 0.7
    radial_bins: int = 100
    vertical_bins: int = 100
    dust_to_gas: float = 0.01

    @u.quantity_input
    def gas_density(self, r: u.au, z: u.au) -> u.g / u.cm ** 3:
        """Calculates gas density at given r, z"""
        raise CHEFNotImplementedError

    @u.quantity_input
    def dust_temperature(self, r: u.au, z: u.au) -> u.K:
        """Calculates dust temperature at given r, z"""
        raise CHEFNotImplementedError

    @u.quantity_input
    def dust_density(self, r: u.au, z: u.au) -> u.g / u.cm ** 3:
        """Calculates dust density at given r, z"""
        return self.gas_density(r, z) * self.dust_to_gas

    @u.quantity_input
    def gas_temperature(self, r: u.au, z: u.au) -> u.K:
        """Calculates gas temperature at given r, z

        Returns:
            dust temperature
        """
        return self.dust_temperature(r, z)

    @cached_property
    def table(self) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        z2r, r = np.meshgrid(
            np.linspace(0, self.zr_max, self.vertical_bins),
            np.geomspace(self.r_min.to(u.au).value, self.r_max.to(u.au).value, self.radial_bins),

        )
        table = CTable(
            [
                (r * u.au).flatten(),
                (r * z2r * u.au).flatten()
            ],
            names=["Radius", "Height"]
        )
        table["Height to radius"] = u.Quantity(z2r.flatten())
        table["Gas density"] = self.gas_density(table.r, table.z)
        table["Dust density"] = self.dust_density(table.r, table.z)
        table["Gas temperature"] = self.gas_temperature(table.r, table.z)
        table["Dust temperature"] = self.dust_temperature(table.r, table.z)
        return table

Ancestors

• PhysicsBase

Subclasses

• ParametrizedPhysics

Class variables

var dust_to_gas : float

var r_max : Unit("AU")

var r_min : Unit("AU")

var radial_bins : int

var vertical_bins : int

var zr_max : float

Instance variables

var table

Return the associated diskchef.CTable with r, z, and dust and gas properties

Expand source code

def __get__(self, instance, owner=None):
    if instance is None:
        return self
    if self.attrname is None:
        raise TypeError(
            "Cannot use cached_property instance without calling __set_name__ on it.")
    try:
        cache = instance.__dict__
    except AttributeError:  # not all objects have __dict__ (e.g. class defines slots)
        msg = (
            f"No '__dict__' attribute on {type(instance).__name__!r} "
            f"instance to cache {self.attrname!r} property."
        )
        raise TypeError(msg) from None
    val = cache.get(self.attrname, _NOT_FOUND)
    if val is _NOT_FOUND:
        with self.lock:
            # check if another thread filled cache while we awaited lock
            val = cache.get(self.attrname, _NOT_FOUND)
            if val is _NOT_FOUND:
                val = self.func(instance)
                try:
                    cache[self.attrname] = val
                except TypeError:
                    msg = (
                        f"The '__dict__' attribute on {type(instance).__name__!r} instance "
                        f"does not support item assignment for caching {self.attrname!r} property."
                    )
                    raise TypeError(msg) from None
    return val

Methods

def dust_density(self, r: Unit("AU"), z: Unit("AU")) ‑> Unit("g / cm3")

Calculates dust density at given r, z

def dust_temperature(self, r: Unit("AU"), z: Unit("AU")) ‑> Unit("K")

Calculates dust temperature at given r, z

def gas_density(self, r: Unit("AU"), z: Unit("AU")) ‑> Unit("g / cm3")

Calculates gas density at given r, z

def gas_temperature(self, r: Unit("AU"), z: Unit("AU")) ‑> Unit("K")

Calculates gas temperature at given r, z

Returns

dust temperature

Inherited members

• PhysicsBase:
  • column_density_to
  • cosmic_ray_padovani18
  • ionization
  • plot_column_densities
  • plot_density
  • plot_temperatures
  • xray_bruderer