Module diskchef.chemistry.andes

Classes

class Reaction (reac_table, idx)

Expand source code

class Reaction:

    def __init__(self, reac_table, idx):
        mask_reac = reac_table[:,0].astype(int) == idx
        reacstr = reac_table[mask_reac][0]
        ireac, r1, r2, p1, p2, p3, p4, p5, alpha, beta, gamma, add_inf = reacstr
        if int(ireac) != idx:
            raise IndexError("Something was read incorrectly")
        self.ireac = idx
        self.r1 = r1.strip()
        self.r2 = r2.strip()
        self.p1 = p1.strip()
        self.p2 = p2.strip()
        self.p3 = p3.strip()
        self.p4 = p4.strip()
        self.p5 = p5.strip()
        self.alpha = float(alpha)
        self.beta = float(beta)
        self.gamma = float(gamma)
        self.get_rtype()

    def __str__(self):
        reacstr = f"{self.r1} + {self.r2} --> {self.p1}"
        for p in [self.p2, self.p3, self.p4, self.p5]:
            if p != '':
                reacstr += f" + {p}"
        reacstr += f"\n rtype: {self.rtype}"
        return reacstr

    def get_rtype(self):
        rtype = 1

        if self.r2 == 'CRP':
            rtype = 2
        if self.r2 == 'PHOTON':
            rtype = 3
            self.alpha *= 1/u.s
            if (self.r1 == 'H2'):
                rtype = 4
            if (self.r1 == 'CO'):
                rtype = 5
        if self.r2 == 'CRPHOT':
            rtype = 6
        if self.r2 == 'G-':
            rtype = 7
        if self.r2 == 'G0':
            rtype = 8
        if self.r1 == 'G0':
            rtype = 9
        if self.r1 == 'G+':
            rtype = 10
        if self.r2 == 'FREEZE':
            rtype = 11
        if self.r2 == 'DESORB':
            rtype = 12
            self.gamma *= u.K
        if (self.r1[:1] == 'g') & (self.r2[:1] == 'g'):
            rtype = 13
            self.gamma *= u.K
        if (self.r1[:1] == 'g') & (self.r2 == 'PHOTON'):
            rtype = 14
        if (self.r1[:1] == 'a') & (self.r2 == 'ASTEROID'):
            rtype = 15

        if rtype==1:
            self.alpha *= u.cm ** 3 / u.s
            self.gamma *= u.K

        self.rtype = rtype

Methods

def get_rtype(self)

class ReadAndesData (physics: PhysicsBase = None, molar_mass: Unit("g / mol") = <Quantity 2.33 g / mol>, hydrogen_mass_fraction: float = 0.739, folder: Union[str, os.PathLike] = None, index: int = 0, read_uv: bool = False, read_ionization: bool = False, read_av: bool = False, read_envelope: bool = False, read_settling: bool = False, nonset_env: bool = False, ionization_format: str = 'ascii.commented_header', age: Unit("Myr") = <Quantity 1. Myr>)

Class to read ANDES2 output for following usage in diskchef.maps

Expand source code

class ReadAndesData(ChemistryBase):
    """
    Class to read ANDES2 output for following usage in `diskchef.maps`
    """
    folder: PathLike = None
    index: int = 0
    read_uv: bool = False
    read_ionization: bool = False
    read_av: bool = False
    read_envelope: bool = False
    read_settling: bool = False
    nonset_env: bool = False
    ionization_format: str = "ascii.commented_header"
    age: u.Myr = 1 * u.Myr

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        if index is None:
            index = self.index
        reader = astropy.io.ascii.get_reader(Reader=astropy.io.ascii.CommentedHeader)
        reader.header.splitter.delimiter = '|'
        chemistry = reader.read(os.path.join(self.folder, f"Chemistry_{index:05d}"))
        if index == 0:
            physics = reader.read(os.path.join(self.folder, f"physical_structure_{index:05d}"))
        else:
            physics = reader.read(os.path.join(self.folder, f"physical_structure_{index-1:05d}"))
        # if physics[["ir", "iz"]] != chemistry[["ir", "iz"]]:
        #     raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and chemistry tables do not match!")
        self._config = self._config_read(os.path.join(self.folder, "../config.ini"))
        self._network_read(os.path.join(self.folder, "../input/Chemical_Network_ALCHEMIC_KIDAbench.dat"))
        self._time_grid = reader.read(os.path.join(self.folder, "grid_time"))
        self.age = float(self._config['last_time_mark  [Myr]'])*u.Myr
        self.amin = float(self._config["minimal initial grain radius [cm]"]) * u.cm
        self.amax = float(self._config["maximal initial grain radius [cm]"]) * u.cm
        self.rho_mon = float(self._config['grain monomer solid density [g/cm3]']) * u.g/u.cm**3
        table = CTable(
            [
                physics["Radius"] << u.au,
                physics["Height"] << u.au,
            ],
            names=["Radius", "Height"]
        )
        table["Height to radius"] = np.nan * u.dimensionless_unscaled
        for izr in physics["iz"]:
            indices = physics["iz"] == izr
            table["Height to radius"][indices] = np.nanmedian(physics["Height"][indices] / physics["Radius"][indices])


        table["Height"] = table.zr * table.r
        table["Distance to star"] = np.sqrt(table.r**2 + table.z**2)
        table["Gas density"] = physics["Gas density"] << (u.g / u.cm ** 3)
        table["Dust density"] = physics["Dust density"] << (u.g / u.cm ** 3)
        table["Gas temperature"] = physics["Gas temperature"] << (u.K)
        table["Dust temperature"] = physics["Dust temperature"] << (u.K)
        table["Av"] = physics["AV_true"] << (u.dimensionless_unscaled)
        if "H2" not in chemistry.colnames and ("oH2" in chemistry.colnames and "pH2" in chemistry.colnames):
            chemistry["H2"] = chemistry["oH2"] + chemistry["pH2"]
        table["n(H+2H2)"] = (chemistry["H+"] + chemistry["H"] + 2 * chemistry["H2"]) << (u.cm ** (-3))
        for species in chemistry.colnames[3:]:
            table[species] = (chemistry[species] << (u.cm ** (-3))) / table["n(H+2H2)"]



        if self.read_uv:
            try:
                radiation = reader.read(os.path.join(self.folder, f"RadInt_StrengthTotal_{index:05d}"))
            except FileNotFoundError:
                self.logger.info("Radiation strength file %:05d not found! Trying 00000 instead", index)
                radiation = reader.read(os.path.join(self.folder, f"RadInt_StrengthTotal_00000"))
            if np.any(radiation[["ir", "iz"]] != physics[["ir", "iz"]]):
                raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and radiation tables do not match!")
            table["G_UV"] = radiation["G_UV"]
            table["G_H2"] = radiation["G_H2"]
            table["G_VIS"] = radiation["G_VIS"]

        if self.read_av:
            try:
                Av = astropy.table.Table.read(os.path.join(self.folder, f"AV_{index:05d}"), format='ascii')
            except FileNotFoundError:
                self.logger.info("Av file %:05d not found! Trying 00000 instead", index)
                radiation = astropy.table.Table.read(os.path.join(self.folder, f"Av_00000"), format='ascii')
            if np.any(radiation[["ir", "iz"]] != physics[["ir", "iz"]]):
                raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and Av tables do not match!")
            table["Av_true"] = Av["AV_true"]
            table["Av_G"] = Av["AV_from_Gfactor"]

        if self.read_envelope:
            try:
                envelope_mask_table = reader.read(os.path.join(self.folder, f"envelope_{index:05d}"))
                if np.any(envelope_mask_table[["ir", "iz"]] != physics[["ir", "iz"]]):
                    raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and envelope tables do not match!")
            except FileNotFoundError:
                self.logger.info("Envelope mask file at given moment not found. Assuming static envelope")
                try:
                    envelope_mask_table = reader.read(os.path.join(self.folder, f"envelope_{0:05d}"))
                    if np.any(envelope_mask_table[["ir", "iz"]] != physics[["ir", "iz"]]):
                        raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and envelope tables do not match!")
                except FileNotFoundError:
                    self.logger.info("Envelope mask file at zero moment not found. Assuming no envelope")
                    envelope_mask_table = np.zeros(len(table['Gas density']))
            try:
                table["Envelope mask"] = envelope_mask_table['is_envelope']
            except IndexError:
                table["Envelope mask"] = envelope_mask_table

            if not self.read_settling:
                dust_types = np.sort(np.unique(table["Envelope mask"]))
                if len(dust_types) > 1:
                    for dt in dust_types:
                        dtmask = table["Envelope mask"] == dt
                        table[f"Dust density {dt}"] = table[f"Dust density"]
                        table[f"Dust temperature {dt}"] = table[f"Dust temperature"]
                        table[f"Dust density {dt}"][~dtmask] = 0
                        table[f"Dust temperature {dt}"][~dtmask] = 0
        else:
            table["Envelope mask"] = np.zeros(len(table['Gas density']))

        if self.read_settling:
            if self.read_envelope:
                self.logger.info(
                    "Settling reading is on while envelope reading is on too. "
                    "\nThese two additions are very similar and will produce conflicting results."
                    "\nFor now we turn off envelope dust types if settling is present")
            if 'settled_dust' in physics.columns:
                dust_types = np.sort(np.unique(physics["settled_dust"]))
                if len(dust_types) > 1:
                    for dt in dust_types:
                        dtmask = physics["settled_dust"] == dt
                        table[f"Dust density {dt}"] = table[f"Dust density"]
                        table[f"Dust temperature {dt}"] = table[f"Dust temperature"]
                        table[f"Dust density {dt}"][~dtmask] = 0
                        # table[f"Dust temperature {dt}"][~dtmask] = 0
                if not self.read_envelope and self.nonset_env:
                    table["Envelope mask"] = 1 - physics["settled_dust"]
                else:
                    table["Settled mask"] = physics["settled_dust"]

            else:
                self.logger.info("Settling reading is turned on but no settling column found. Proceeding with no settling.")


        if 'Settled mask' in table.columns:
            table['Average dust size'] = np.sqrt(self.amin * self.amax * (1000 * table["Settled mask"] + (1 - table["Settled mask"])))
        else:
            table['Average dust size'] = np.sqrt(self.amin * self.amax)

        table['Dust number density'] = table['Dust density']/(4/3 * np.pi * table['Average dust size']**3 * self.rho_mon)

        if self.read_ionization:
            try:
                ionization = astropy.table.Table.read(os.path.join(self.folder, f"Ionization_Rate_{index:05d}"))
            except FileNotFoundError:
                if os.path.exists(os.path.join(self.folder, f"Ionization_Rate")):
                    self.logger.info("Ionization rate file %:05d not found! Ionization_Rate is found instead", index)
                    ionization = astropy.table.Table.read(
                        os.path.join(self.folder, f"Ionization_Rate"),
                        format=self.ionization_format
                    )
                else:
                    self.logger.info("Ionization rate file %:05d not found! Trying 00000 instead", index)
                    ionization = astropy.table.Table.read(
                        os.path.join(self.folder, f"Ionization_Rate_00000"),
                        format=self.ionization_format
                    )
            if np.any(ionization[["ir", "iz"]] != physics[["ir", "iz"]]):
                raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and ionization tables do not match!")
            table["X ray ionization rate"] = ionization["IR_XR[1/s]"] << u.s ** (-1)
            table["CR ionization rate"] = ionization["IR_CR[1/s]"] << u.s ** (-1)
            table["Ionization rate"] = ionization["IR_Total[1/s]"] << u.s ** (-1)

        table.sort(['Height to radius', 'Radius'])

        self.physics = diskchef.physics.base.PhysicsBase(
            star_mass=float(self._config['stellar mass [MSun]']) * u.solMass,
        )
        self.physics.table = table
        return table

    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()
        self.update_hydrogen_atom_number_density()

    def _config_read(self, path: PathLike) -> dict:
        out = {}
        with open(path, 'r') as configfile:
            for line in configfile.readlines():
                value, key = [entry.strip() for entry in line.split('::')]
                out[key] = value
        return out

    def _network_read(self, path: PathLike):
        nsp = int(np.loadtxt(path, max_rows=1))
        self.species = np.loadtxt(path, max_rows=nsp, skiprows=1,
                              dtype='str')
        nreac = int(np.loadtxt(path, max_rows=1, skiprows=nsp + 1))
        self.reactions = np.genfromtxt(path, skip_header=nsp + 2, dtype='str',
                              max_rows=nreac,  # delimiter='/n'
                              delimiter=(6, 10, 20, 10, 10, 10, 10, 10, 9, 6, 9, 2)
                              )

    def calc_rates(self, reac):
        Cion = 1 + (1.671e-3 * u.cm * u.K) / self.table["Average dust size"] / self.table["Gas temperature"]
        photon_flux_draine = 0.463e8 / u.cm ** 2 / u.s
        stick_H = 1 / (
                1 + 0.4 * np.sqrt(2 * self.table["Gas temperature"] / (100 * u.K)) + 0.2 * self.table["Gas temperature"] / (100 * u.K) + 0.08 * (self.table["Gas temperature"] / (100 * u.K)) ** 2)
        site_size_ice = 2.6e-8 * u.cm
        site_density = 4e14 / u.cm ** 2
        ebed = 0.5
        rdratio = 0.01

        if reac.rtype == 1:
            rate = reac.alpha * (self.table["Gas temperature"] / (300 * u.K)) ** reac.beta * np.exp(-reac.gamma / self.table["Gas temperature"]) * self.table["Dust number density"]
        if reac.rtype == 2:
            rate = reac.alpha * self.table["Ionization rate"]
        if reac.rtype == 3:
            rate = reac.alpha * np.exp(-reac.gamma * self.table["Av_G"]) * self.table["G_UV"]
        if reac.rtype == 4 or reac.rtype==5:
            f_H2, f_CO = self.calc_H2_CO_shielding()
            H2_shield_coeff = np.repeat(f_H2[np.newaxis, :], len(self.table.zr_grid), axis=0)
            H2_shield_coeff = H2_shield_coeff.flatten()
            CO_shield_coeff = np.repeat(f_CO[np.newaxis, :], len(self.table.zr_grid), axis=0)
            CO_shield_coeff = CO_shield_coeff.flatten()
            if reac.rtype == 4:
                rate = self.table["G_H2"] * reac.alpha * H2_shield_coeff
            if reac.rtype == 5:
                rate = self.table["G_VIS"] * reac.alpha * CO_shield_coeff
        if reac.rtype == 6:
            rate = reac.alpha * self.table["Ionization rate"]
            if (reac.r1[:1] == 'g') & (self._config["nocrphotice"] == '.true.'):
                rate = 0
        if reac.rtype == 7:
            weight = mol_weight(reac.r1)
            Vth = np.sqrt(8 * c.k_B * self.table["Gas temperature"] / np.pi / weight)
            rate = reac.alpha * np.pi * self.table["Average dust size"] ** 2 * Vth * Cion * self.table["Dust number density"]
        if reac.rtype == 8:
            weight = mol_weight(reac.r1)
            Vth = np.sqrt(8 * c.k_B * self.table["Gas temperature"] / np.pi / weight)
            rate = reac.alpha * np.pi * self.table["Average dust size"] ** 2 * Vth * self.table["Dust number density"]
        if reac.rtype == 9:
            Vth = np.sqrt(8 * c.k_B * self.table["Gas temperature"] / np.pi / c.m_e)
            rate = np.pi * self.table["Average dust size"] ** 2 * Vth * self.table["Dust number density"]
        if reac.rtype == 10:
            Vth = np.sqrt(8 * c.k_B * self.table["Gas temperature"] / np.pi / c.m_e)
            rate = np.pi * self.table["Average dust size"] ** 2 * Vth * Cion * self.table["Dust number density"]
        if reac.rtype == 11:
            weight = mol_weight(reac.r1)
            Vth = np.sqrt(8 * c.k_B * self.table["Gas temperature"] / np.pi / weight)
            rate = reac.alpha * np.pi * self.table["Average dust size"] ** 2 * Vth * self.table["Dust number density"]

            if (reac.r1 == 'H'):
                rate = reac.alpha * np.pi * self.table["Average dust size"] ** 2 * Vth * self.table["Dust number density"] * stick_H
            if (reac.r1 == 'H2') | (reac.r1 == 'pH2') | (reac.r1 == 'oH2'):
                rate = 0
        if reac.rtype == 12:
            photodesorption_yield = 3.5e-3 + 0.13 * np.exp(-336. * u.K / self.table["Dust temperature"])
            ak_photo = photodesorption_yield * photon_flux_draine * self.table["G_UV"] * site_size_ice ** 2 / 4

            weight = mol_weight(reac.r1)
            anu0 = np.sqrt(2 * site_density * c.k_B * reac.gamma / np.pi ** 2 / weight)
            ak_therm = (anu0 * np.exp(-reac.gamma / self.table["Dust temperature"])).to(1 / u.s)

            tcrp = (4.36e5 * u.K ** 3 + self.table["Dust temperature"] ** 3) ** (1 / 3.)
            ak_crp = ((2.431e-2 * u.s) * self.table["Ionization rate"] * anu0 * np.exp(-reac.gamma / tcrp)).to(1 / u.s)

            rate = ak_crp + ak_photo + ak_therm
        if reac.rtype == 13:
            weight1 = mol_weight(reac.r1)
            mask_r1 = (self.reactions[:, 1] == reac.r1.ljust(10)) & (self.reactions[:, 2] == 'DESORB'.ljust(20))
            try:
                edes1 = float(self.reactions[mask_r1][0][10]) * u.K
            except IndexError:
                edes1 = 0 * u.K
            anu0 = np.sqrt(2 * site_density * c.k_B * edes1 / np.pi ** 2 / weight1)

            weight2 = mol_weight(reac.r2)
            mask_r2 = (self.reactions[:, 1] == reac.r2.ljust(10)) & (self.reactions[:, 2] == 'DESORB'.ljust(20))
            try:
                edes2 = float(self.reactions[mask_r2][0][10]) * u.K
            except IndexError:
                edes2 = 0 * u.K
            anu1 = np.sqrt(2 * site_density * c.k_B * edes2 / np.pi ** 2 / weight2)

            Rdiff0 = anu0 / (4 * np.pi * self.table["Average dust size"] ** 2 * site_density) * np.exp(-ebed * edes1 / self.table["Dust temperature"])
            Rdiff1 = anu1 / (4 * np.pi * self.table["Average dust size"] ** 2 * site_density) * np.exp(-ebed * edes2 / self.table["Dust temperature"])
            akbar = np.exp(-reac.gamma / self.table["Dust temperature"])

            if (reac.r1 == 'gH') | (reac.r1 == 'gH2') | (reac.r2 == 'gH') | (reac.r2 == 'gH2'):
                amur = weight1 * weight2 / (weight1 + weight2)
                akbar = np.exp(-4 * np.pi / c.h * (1e-8 * u.cm) * np.sqrt(2 * amur * reac.gamma * c.k_B))

            rate = akbar * (Rdiff0 + Rdiff1)
            if reac.r1 == reac.r2:
                rate *= 0.5
            if reac.p1[:1] == 'g':
                rate *= 1 - rdratio
            else:
                rate *= rdratio
        if reac.rtype==14:
            rate = 0 / u.s
        if reac.rtype==15:
            rate = 0 / u.s

        try:
            final_rate = rate.to(1 / u.s) * self.table[reac.r1] * self.table[reac.r2] * self.table[
                "n(H+2H2)"] ** 2 / self.table["Dust number density"] ** 2
        except KeyError:
            final_rate = rate.to(1 / u.s) * self.table[reac.r1] * self.table["n(H+2H2)"] / self.table["Dust number density"]
        return final_rate

    def calc_H2_CO_shielding(self):
        z = np.reshape(self.table['Height'], (len(self.table.zr_grid), len(self.table.r_grid)))
        H2_dens = np.reshape(self.table['H2'] * self.table["n(H+2H2)"],
                             (len(self.table.zr_grid), len(self.table.r_grid)))
        CO_dens = np.reshape(self.table['CO'] * self.table["n(H+2H2)"],
                             (len(self.table.zr_grid), len(self.table.r_grid)))
        H2coldens = np.trapz(H2_dens, x=z, axis=0).to(1 / u.cm ** 2).value
        COcoldens = np.trapz(CO_dens, x=z, axis=0).to(1 / u.cm ** 2).value
        ncd = H2coldens / 5e14
        f_H2 = 0.95 / (1 + ncd) ** 2 + 0.035 / np.sqrt(1 + ncd) * np.exp(-8.5e-4 * np.sqrt(1 + ncd))
        log_CO_grid = [0., 13., 14., 15., 16., 17., 18., 19.]
        log_H2_grid = [0., 19., 20., 21., 22., 23.]
        visser_teta = [[1.000, 8.080e-1, 5.250e-1, 2.434e-1, 5.467e-2, 1.362e-2, 3.378e-3, 5.240e-4],
                       [8.176e-1, 6.347e-1, 3.891e-1, 1.787e-1, 4.297e-2, 1.152e-2, 2.922e-3, 4.662e-4],
                       [7.223e-1, 5.624e-1, 3.434e-1, 1.540e-1, 3.515e-2, 9.231e-3, 2.388e-3, 3.899e-4],
                       [3.260e-1, 2.810e-1, 1.953e-1, 8.726e-2, 1.907e-2, 4.768e-3, 1.150e-3, 1.941e-4],
                       [1.108e-2, 1.081e-2, 9.033e-3, 4.441e-3, 1.102e-3, 2.644e-4, 7.329e-5, 1.437e-5],
                       [3.938e-7, 3.938e-7, 3.936e-7, 3.923e-7, 3.901e-7, 3.893e-7, 3.890e-7, 3.875e-7]]

        f_CO = interpn(
            points=(np.array(log_H2_grid), np.array(log_CO_grid)),
            values=np.array(visser_teta),
            xi=(np.log10(COcoldens), np.log10(H2coldens)),
            method='linear', fill_value=0, bounds_error=False
        )
        f_CO = 10 ** f_CO

        return [f_H2, f_CO]

    def get_top_reac(self, spec, ir, iz, limfrac=0.01):
        nrp = len(self.table.r_grid)
        mask_destreac = (self.reactions[:, 1] == spec.ljust(10)) | (self.reactions[:, 2] == spec.ljust(20))
        mask_formreac = (self.reactions[:, 3] == spec.ljust(10)) | (self.reactions[:, 4] == spec.ljust(10)) | (
                    self.reactions[:, 5] == spec.ljust(10)) | (self.reactions[:, 6] == spec.ljust(10)) | (self.reactions[:, 7] == spec.ljust(10))
        n_dest = np.sum(mask_destreac)
        n_form = np.sum(mask_formreac)
        ireacs_form = []
        ireacs_dest = []
        rates_form = []
        rates_dest = []
        for ireac in self.reactions[mask_destreac][:,0].astype(int):
            reac = Reaction(self.reactions, ireac)
            ireacs_dest.append(ireac)
            rates_dest.append(self.calc_rates(reac)[iz * nrp + ir].value)
        for ireac in self.reactions[mask_formreac][:,0].astype(int):
            reac = Reaction(self.reactions, ireac)
            ireacs_form.append(ireac)
            rates_form.append(self.calc_rates(reac)[iz * nrp + ir].value)

        rates_form = np.array(rates_form)
        ireacs_form = np.array(ireacs_form, dtype=int)

        totdestrate = np.sum(rates_dest)
        totformrate = np.sum(rates_form)

        rates_form_top = rates_form[rates_form > limfrac * totformrate]
        ireacs_form_top = ireacs_form[rates_form > limfrac * totformrate]

        rates_dest = np.array(rates_dest)
        ireacs_dest = np.array(ireacs_dest, dtype=int)
        rates_dest_top = rates_dest[rates_dest > limfrac * totdestrate]
        ireacs_dest_top = ireacs_dest[rates_dest > limfrac * totdestrate]

        form_table = QTable([ireacs_form_top, rates_form_top], names=('ireac', 'rate'))
        dest_table = QTable([ireacs_dest_top, -rates_dest_top], names=('ireac', 'rate'))
        form_table.sort('rate')
        form_table.reverse()
        dest_table.sort('rate')
        dest_table.reverse()
        final_table = vstack([form_table, dest_table])
        final_table.meta['totformrate'] = totformrate
        final_table.meta['totdestrate'] = -totdestrate

        return final_table

    def vis_top_reac(self, spec, ir, iz, limfrac=0.01, save_path: str = None,
                     figs_x: float = 12, figs_y: float = 9, show_abundance: bool = False, show_location: bool = False,
                     show_sum: bool = False):
        top_reac_table = self.get_top_reac(spec, ir, iz, limfrac=limfrac)
        formations = []
        destructions = []
        destruction_partners = []
        for ireac in top_reac_table['ireac'][top_reac_table['rate'] > 0]:
            reac = Reaction(self.reactions, ireac)
            formations.append([reac.r1, reac.r2])
        formation_rates = top_reac_table['rate'][top_reac_table['rate'] > 0].data
        for ireac in top_reac_table['ireac'][top_reac_table['rate'] < 0]:
            reac = Reaction(self.reactions, ireac)
            subd = []
            for p in [reac.p1, reac.p2, reac.p3, reac.p4, reac.p5]:
                if p != '':
                    subd.append(p)
            destructions.append(subd)
            if reac.r1 == spec:
                destruction_partners.append(reac.r2)
            else:
                destruction_partners.append(reac.r1)
        destruction_rates = top_reac_table['rate'][top_reac_table['rate'] < 0].data

        G = nx.DiGraph()

        node_colors = {}
        node_sizes = {}
        node_labels = {}
        edge_labels = {}

        G.add_node(spec)
        node_colors[spec] = 'grey'
        node_sizes[spec] = 3000
        node_labels[spec] = spec

        pos = {spec: (0,0)}

        for i, reactants in enumerate(formations):
            reac_node = f"F_reac_{i}"
            G.add_node(reac_node)
            node_colors[reac_node] = "deepskyblue"
            node_sizes[reac_node] = 100
            node_labels[reac_node] = ""

            pos[reac_node] = (-1, (i - (len(formations) - 1) / 2) * 1.5)

            G.add_edge(reac_node, spec, color="deepskyblue", weight=3)
            edge_labels[(reac_node, spec)] = f'{formation_rates[i]:.3e}'

            for j, reactant in enumerate(reactants):
                left_node_id = f"{reactant}_left"
                if left_node_id not in G:
                    G.add_node(left_node_id)
                    node_colors[left_node_id] = "dodgerblue"
                    node_sizes[left_node_id] = 2000
                    node_labels[left_node_id] = reactant

                reac_x, reac_y = pos[reac_node]
                offset_y = reac_y + (j - (len(reactants) - 1) / 2) * 0.5
                pos[left_node_id] = (-2, offset_y)

                G.add_edge(left_node_id, reac_node, color="dodgerblue", weight=1.5)

        for i, products in enumerate(destructions):
            reac_node = f"D_reac_{i}"
            G.add_node(reac_node)
            node_colors[reac_node] = "orange"
            node_sizes[reac_node] = 2500
            node_labels[reac_node] = f"+\n{destruction_partners[i]}"

            pos[reac_node] = (1, (i - (len(destructions) - 1) / 2) * 1.5)

            G.add_edge(spec, reac_node, color="orange", weight=3)
            edge_labels[(spec, reac_node)] = f'{destruction_rates[i]:.3e}'

            for j, product in enumerate(products):
                right_node_id = f"{product}_right"
                if right_node_id not in G:
                    G.add_node(right_node_id)
                    node_colors[right_node_id] = "darkorange"  # Orange
                    node_sizes[right_node_id] = 2000
                    node_labels[right_node_id] = product

                reac_x, reac_y = pos[reac_node]
                offset_y = reac_y + (j - (len(products) - 1) / 2) * 0.5
                pos[right_node_id] = (2, offset_y)

                G.add_edge(reac_node, right_node_id, color="darkorange", weight=1.5)

        fig, ax = plt.subplots(figsize=(figs_x, figs_y))

        colors = [node_colors[node] for node in G.nodes()]
        sizes = [node_sizes[node] for node in G.nodes()]

        nx.draw_networkx_nodes(
            G, pos, ax=ax, node_color=colors, node_size=sizes, alpha=0.9
        )

        edge_colors = [G[u][v]["color"] for u, v in G.edges()]
        edge_widths = [G[u][v]["weight"] for u, v in G.edges()]

        nx.draw_networkx_edges(
            G,
            pos,
            ax=ax,
            edge_color=edge_colors,
            width=edge_widths,
            arrowsize=18,
            arrowstyle="-|>",
            node_size=sizes,
        )

        nx.draw_networkx_labels(
            G, pos, labels=node_labels, ax=ax, font_size=11, font_family=rcParams['font.family'], font_weight="bold"
        )

        nx.draw_networkx_edge_labels(
            G,
            pos,
            edge_labels=edge_labels,
            ax=ax,
            font_size=10,
            font_color="black",
            rotate=False,
            font_family = rcParams['font.family'], font_weight="bold"
        )

        title_str = f"Top reactions for {spec}\n"
        nrp = len(self.table.r_grid)
        if show_location:
            title_str += f"$t$ = {self._time_grid['t[yr]'][self.index]} yr, $r$ = {self.table['Radius'][iz * nrp + ir].to(u.au).value:.2f} au, $z/r$ = {self.table['Height to radius'][iz * nrp + ir]:.2f}\n"
        if show_abundance:
            title_str += f"Abundance: {self.table[spec][iz * nrp + ir]:.3e}\n"
        title_str += f"Total formation rate = {top_reac_table.meta['totformrate']:.2e} 1/s\n"
        title_str += f"Total destruction rate = {top_reac_table.meta['totdestrate']:.2e} 1/s\n"
        if show_sum:
            mol_change = (top_reac_table.meta['totformrate']+top_reac_table.meta['totdestrate'])
            if mol_change>0:
                title_str += f"Formation {mol_change:.2e} 1/s"
            elif mol_change<0:
                title_str += f"Destruction {mol_change:.2e} 1/s"
            else:
                title_str += f"Perfect balance"
        ax.set_title(
            title_str,
            fontsize=14,
            pad=20,
        )
        plt.axis("off")
        plt.tight_layout()
        if save_path is not None:
            fig.savefig(save_path, bbox_inches='tight', dpi=300)
        else:
            plt.show()
        plt.close(fig=fig)

Ancestors

• ChemistryBase

Class variables

var age : Unit("Myr")

var folder : Union[str, os.PathLike]

var index : int

var ionization_format : str

var nonset_env : bool

var read_av : bool

var read_envelope : bool

var read_ionization : bool

var read_settling : bool

var read_uv : bool

Methods

def calc_H2_CO_shielding(self)

def calc_rates(self, reac)

def get_top_reac(self, spec, ir, iz, limfrac=0.01)

def read(self, index: int = None) ‑> CTable

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

def vis_top_reac(self, spec, ir, iz, limfrac=0.01, save_path: str = None, figs_x: float = 12, figs_y: float = 9, show_abundance: bool = False, show_location: bool = False, show_sum: bool = False)

Inherited members

• ChemistryBase:
  • add_co_isotopologues
  • calculate_column_density_towards_star
  • calculate_column_density_upwards
  • run_chemistry
  • table
  • update_hydrogen_atom_number_density

class ReadAndesHurakanData (physics: PhysicsBase = None, molar_mass: Unit("g / mol") = <Quantity 2.33 g / mol>, hydrogen_mass_fraction: float = 0.739, folder: Union[str, os.PathLike] = None, index: int = 0, read_uv: bool = False, read_ionization: bool = False, read_envelope: bool = False, ionization_format: str = 'ascii.commented_header', age: Unit("Myr") = <Quantity 1. Myr>)

Class for a special ANDES-HURAKAN pipeline

Expand source code

class ReadAndesHurakanData(ChemistryBase):
    """
    Class for a special ANDES-HURAKAN pipeline
    """
    folder: PathLike = None
    index: int = 0
    read_uv: bool = False
    read_ionization: bool = False
    read_envelope: bool = False
    ionization_format: str = "ascii.commented_header"
    age: u.Myr = 1 * u.Myr

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None, stage: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        if index is None:
            index = self.index
        reader = astropy.io.ascii.get_reader(Reader=astropy.io.ascii.CommentedHeader)
        reader.header.splitter.delimiter = '|'
        chemistry = reader.read(os.path.join(self.folder, f"Chemistry_{index:04d}"))
        # if physics[["ir", "iz"]] != chemistry[["ir", "iz"]]:
        #     raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and chemistry tables do not match!")

        self._config = self._config_read(f"setup.inp")
        self.T_star = float(self._config["Tstar[K]"]) * u.K
        self.R_star = float(self._config["Rstar[RSun]"]) * u.Rsun
        self.Tacc = float(self._config["Tacc[K]"]) * u.K

        self._flare = astropy.io.ascii.read(self._config["flare_profile_file"], names=['i', 't', 'L'])
        self.Lacc = self._flare['L'][index-1] * u.solLum

        table = CTable(
            [
                chemistry["Radius"] << u.au,
                chemistry["Relative Height"] << u.dimensionless_unscaled,
            ],
            names=["Radius", "Height to radius"]
        )

        table["Height"] = table.zr * table.r
        table["Distance to star"] = np.sqrt(table.r**2 + table.z**2)
        table["Gas density"] = chemistry["rho_gas"] << (u.g / u.cm ** 3)
        table["Dust density"] = table["Gas density"]/100
        table["Gas temperature"] = chemistry["Tgas"] << (u.K)
        table["Dust temperature"] = chemistry["Tgas"] << (u.K)
        table["n(H+2H2)"] = (chemistry["H+"] + chemistry["H"] + 2 * chemistry["H2"]) << (u.cm ** (-3))
        for species in chemistry.colnames[9:]:
            table[species] = (chemistry[species] << (u.cm ** (-3))) / table["n(H+2H2)"]

        # if self.read_uv:
        #     try:
        #         radiation = reader.read(os.path.join(self.folder, f"RadInt_StrengthTotal_{index:05d}"))
        #     except FileNotFoundError:
        #         self.logger.info("Radiation strength file %:05d not found! Trying 00000 instead", index)
        #         radiation = reader.read(os.path.join(self.folder, f"RadInt_StrengthTotal_00000"))
        #     table["G_UV"] = radiation["G_UV"] << diskchef.maps.radiation_fields.ANDES2_G0
        #
        # if self.read_envelope:
        #     try:
        #         envelope_mask_table = reader.read(os.path.join(self.folder, f"envelope_{index:05d}"))
        #     except FileNotFoundError:
        #         self.logger.info("Envelope mask file at given moment not found. Assuming static envelope")
        #         try:
        #             envelope_mask_table = reader.read(os.path.join(self.folder, f"envelope_{0:05d}"))
        #         except FileNotFoundError:
        #             self.logger.info("Envelope mask file at zero moment not found. Assuming no envelope")
        #             envelope_mask_table = np.zeros(len(table['Gas density']))
        #     try:
        #         table["Envelope mask"] = envelope_mask_table['is_envelope']
        #     except IndexError:
        #         table["Envelope mask"] = envelope_mask_table
        # else:
        #     table["Envelope mask"] = np.zeros(len(table['Gas density']))
        #
        # if self.read_ionization:
        #     try:
        #         ionization = astropy.table.Table.read(os.path.join(self.folder, f"Ionization_Rate_{index:05d}"))
        #     except FileNotFoundError:
        #         if os.path.exists(os.path.join(self.folder, f"Ionization_Rate")):
        #             self.logger.info("Ionization rate file %:05d not found! Ionization_Rate is found instead", index)
        #             ionization = astropy.table.Table.read(
        #                 os.path.join(self.folder, f"Ionization_Rate"),
        #                 format=self.ionization_format
        #             )
        #         else:
        #             self.logger.info("Ionization rate file %:05d not found! Trying 00000 instead", index)
        #             ionization = astropy.table.Table.read(
        #                 os.path.join(self.folder, f"Ionization_Rate_00000"),
        #                 format=self.ionization_format
        #             )
        #     table["X ray ionization rate"] = ionization["IR_XR[1/s]"] << u.s ** (-1)
        #     table["CR ionization rate"] = ionization["IR_CR[1/s]"] << u.s ** (-1)
        #     table["Ionization rate"] = ionization["IR_Total[1/s]"] << u.s ** (-1)

        self.physics = diskchef.physics.base.PhysicsBase(
            star_mass=float(self._config['Mstar[MSun]']) * u.Msun,
        )
        self.physics.table = table
        return table

    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()
        self.update_hydrogen_atom_number_density()

    def _config_read(self, path: PathLike) -> dict:
        out = {}
        with open(path, 'r') as configfile:
            for line in configfile.readlines():
                if (len(line.split('#')) > 1) or (line=='\n'):
                    continue
                key, value = [entry.strip() for entry in line.split(':')]
                out[key] = value
        return out

Ancestors

• ChemistryBase

Class variables

var age : Unit("Myr")

var folder : Union[str, os.PathLike]

var index : int

var ionization_format : str

var read_envelope : bool

var read_ionization : bool

var read_uv : bool

Methods

def read(self, index: int = None, stage: int = None) ‑> CTable

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

Inherited members

• ChemistryBase:
  • add_co_isotopologues
  • calculate_column_density_towards_star
  • calculate_column_density_upwards
  • run_chemistry
  • table
  • update_hydrogen_atom_number_density

class ReadAndesPhysData (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, folder: Union[str, os.PathLike] = None, index: int = 0, time_mark: Unit("s") = <Quantity 0. s>, age: Unit("Myr") = <Quantity 1. Myr>, read_envelope: bool = False, dim: Literal['2d', '3d'] = '2d', coord_type: Literal['z', 'zr'] = 'zr')

Class to read a version of ANDES (accounting for jets + no chemistry) output for following usage in diskchef.maps In testing

Expand source code

class ReadAndesPhysData(PhysicsBase):
    """
    Class to read a version of ANDES (accounting for jets + no chemistry) output for following usage in `diskchef.maps`
    In testing
    """
    folder: PathLike = None
    index: int = 0
    time_mark: u.s = 0 * u.s
    age: u.Myr = 1 * u.Myr
    read_envelope: bool = False
    dim: Literal["2d", "3d"] = '2d'
    coord_type: Literal["z", "zr"] = 'zr'

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        if index is None:
            index = self.index
        reader = astropy.io.ascii.get_reader(Reader=astropy.io.ascii.CommentedHeader)
        reader.header.splitter.delimiter = '|'
        # if physics[["ir", "iz"]] != chemistry[["ir", "iz"]]:
        #     raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and chemistry tables do not match!")
        try:
            self._config = self._config_read(os.path.join(self.folder, "../config.ini"))
        except FileNotFoundError:
            self._config = self._config_read(os.path.join(self.folder, "../config_sep.ini"))
        self.age = float(self._config['last_time_mark  [Myr]']) * u.Myr
        if self.dim == '2d':
            try:
                if index == 0:
                    physics = reader.read(os.path.join(self.folder, f"physical_structure_{index:05d}"))
                else:
                    physics = reader.read(os.path.join(self.folder, f"physical_structure_{index - 1:05d}"))
            except FileNotFoundError:
                physics = reader.read(os.path.join(self.folder, f"physical_structure"))
            table = CTable(
                [
                    physics["Radius"] << u.au,
                    physics["Height"] << u.au,
                ],
                names=["Radius", "Height"]
            )
            table["Gas temperature"] = physics["Gas temperature"] << (u.K)
            table["Dust temperature"] = physics["Dust temperature"] << (u.K)
        if self.dim == '3d':
            physics = reader.read(os.path.join(self.folder, f"all_data_{index:05d}"))
            table = CTable(
                [
                    physics["Radius"] << u.au,
                    physics["Height"] << u.au,
                    physics["Angle"] << u.rad,
                ],
                names=["Radius", "Height", "Phi"], dims='3d'
            )
            table["Gas temperature"] = physics["Temperature"] << (u.K)
            table["Dust temperature"] = physics["Temperature"] << (u.K)

        if self.coord_type == 'zr':
            table["Height to radius"] = np.nan * u.dimensionless_unscaled
            for izr in physics["iz"]:
                indices = physics["iz"] == izr
                table["Height to radius"][indices] = np.nanmedian(
                    physics["Height"][indices] / physics["Radius"][indices])
                table["Height"] = table.zr * table.r

        if self.coord_type == 'z':
            table["Height to radius"] = table["Height"]/table.r

        table["Distance to star"] = np.sqrt(table.r ** 2 + table.z ** 2)
        table["Gas density"] = physics["Gas density"] << (u.g / u.cm ** 3)
        table["Dust density"] = physics["Dust density"] << (u.g / u.cm ** 3)


        if self.read_envelope:
            try:
                if self.dim == '2d':
                    envelope_mask_table = astropy.io.ascii.read(os.path.join(self.folder, f"Envelope_test1.dat"))
                    if np.any(envelope_mask_table[["ir", "iz"]] != physics[["ir", "iz"]]):
                        raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and envelope tables do not match!")
                    table["Envelope mask"] = envelope_mask_table['is_envelope']
                elif self.dim == '3d':
                    table["Envelope mask"] = physics["env"]
            except FileNotFoundError:
                self.logger.info("Envelope mask file not found. Assuming no envelope")
                table["Envelope mask"] = np.zeros(len(table['Gas density']))
            dust_types = np.sort(np.unique(table["Envelope mask"]))
            if len(dust_types) > 1:
                for dt in dust_types:
                    dtmask = table["Envelope mask"] == dt
                    table[f"Dust density {dt}"] = table[f"Dust density"]
                    table[f"Dust temperature {dt}"] = table[f"Dust temperature"]
                    table[f"Dust density {dt}"][~dtmask] = 0
                    table[f"Dust temperature {dt}"][~dtmask] = 0
        else:
            table["Envelope mask"] = np.zeros(len(table['Gas density']))

        dust_types = np.sort(np.unique(table["Envelope mask"]))
        if len(dust_types) > 1:
            for dt in dust_types:
                dtmask = table["Envelope mask"] == dt
                table[f"Dust density {dt}"] = table[f"Dust density"]
                table[f"Dust temperature {dt}"] = table[f"Dust temperature"]
                table[f"Dust density {dt}"][~dtmask] = 0
                table[f"Dust temperature {dt}"][~dtmask] = 0

        if self.dim == '2d':
            if self.coord_type == 'zr':
                table.sort(['Height to radius', 'Radius'])
            elif self.coord_type == 'z':
                table.sort(['Height', 'Radius'])
        elif self.dim == '3d':
            if self.coord_type == 'zr':
                table.sort(['Phi', 'Height to radius', 'Radius'])
            elif self.coord_type == 'z':
                table.sort(['Phi', 'Height', 'Radius'])


        try:
            self.physics = diskchef.physics.base.PhysicsBase(
                star_mass=float(self._config['stellar mass [MSun]']) * u.solMass,
            )
        except KeyError:
            self.physics = diskchef.physics.base.PhysicsBase(
                star_mass=(float(self._config['stellar mass 1 [MSun]'])+float(self._config['stellar mass 2 [MSun]'])) * u.solMass,
            )
        self.physics.table = table

        return table

    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()

    def _config_read(self, path: PathLike) -> dict:
        out = {}
        with open(path, 'r') as configfile:
            for line in configfile.readlines():
                if len(line.split('::')) == 1:
                    continue
                value, key = [entry.strip() for entry in line.split('::')]
                out[key] = value
        return out

    def recalc_dust_dens(self):
        dust_types = np.sort(np.unique(self.table["Envelope mask"]))
        if len(dust_types) > 1:
            for dt in dust_types:
                dtmask = self.table["Envelope mask"] == dt
                self.table[f"Dust density {dt}"] = self.table[f"Dust density"]
                self.table[f"Dust temperature {dt}"] = self.table[f"Dust temperature"]
                self.table[f"Dust density {dt}"][~dtmask] = 0
                self.table[f"Dust temperature {dt}"][~dtmask] = 0

Ancestors

• PhysicsBase

Class variables

var age : Unit("Myr")

var coord_type : Literal['z', 'zr']

var dim : Literal['2d', '3d']

var folder : Union[str, os.PathLike]

var index : int

var read_envelope : bool

var time_mark : Unit("s")

Methods

def read(self, index: int = None) ‑> CTable

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

def recalc_dust_dens(self)

Inherited members

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

class ReadHurakanData (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, folder: Union[str, os.PathLike] = None, index: int = 0, time_mark: Unit("s") = <Quantity 0. s>, age: Unit("Myr") = <Quantity 1. Myr>, read_envelope: bool = False)

Class to read Hurakan output for following usage in diskchef.maps Unfinished

Expand source code

class ReadHurakanData(PhysicsBase):
    """
    Class to read Hurakan output for following usage in `diskchef.maps`
    Unfinished
    """
    folder: PathLike = None
    index: int = 0
    time_mark: u.s = 0 * u.s
    age: u.Myr = 1 * u.Myr
    read_envelope: bool = False

    @property
    def table(self) -> CTable:
        return self._table

    def read(self, index: int = None) -> CTable:
        """Return the associated diskchef.CTable with r, z, and dust and gas properties"""
        if index is None:
            index = self.index
        physics = astropy.io.ascii.read(os.path.join(self.folder, f"result-{index:04d}.dat"), header_start=3, data_start=4)#, data_end=-101)
        theta_zeros = physics["Theta"] == 0
        physics["Theta"][theta_zeros] = np.abs(np.pi - physics["Theta"][-1])
        # physics = physics[~theta_zeros]
        table = CTable(
            [
                physics["R"] << u.cm,
                physics["Theta"] << u.rad,
            ],
            names=["Distance to star", "Theta"], coord_base='sphere', coord_style='edge'
        )
        self._config = astropy.io.ascii.read(os.path.join(self.folder, f"result-{index:04d}.dat"), data_end=3)
        self.time_mark = self._config['col2'][1] * u.s
        table["Height"] = table.d * np.cos(table.theta)
        table["Radius"] = table.d * np.sin(table.theta)
        table["Height to radius"] = np.tan(np.pi/2*u.rad - table.theta)

        table["Gas density"] = physics["Rho"] << (u.g / u.cm ** 3)
        table["Dust density"] = table["Gas density"]/100
        table["Gas temperature"] = physics["T"] << (u.K)
        table["Dust temperature"] = table["Gas temperature"]

        if self.read_envelope: #to be done
            try:
                envelope_mask_table = astropy.io.ascii.read(os.path.join(self.folder, f"envelope_{0:05d}"))
                if np.any(envelope_mask_table[["ir", "iz"]] != physics[["ir", "iz"]]):
                    raise diskchef.engine.exceptions.CHEFRuntimeError("Physics and envelope tables do not match!")
                table["Envelope mask"] = envelope_mask_table['is_envelope']
            except FileNotFoundError:
                self.logger.info("Envelope mask file not found. Assuming no envelope")
                table["Envelope mask"] = np.zeros(len(table['Gas density']))
        else:
            table["Envelope mask"] = np.zeros(len(table['Gas density']))

        table.sort(['Height to radius', 'Radius'])

        self.physics = diskchef.physics.base.PhysicsBase(
            star_mass=float(self._config['col2'][2]) * u.g,
        )
        self.physics.table = table

        return table

    def __post_init__(self):
        super().__post_init__()
        self._table = self.read()

Ancestors

• PhysicsBase

Class variables

var age : Unit("Myr")

var folder : Union[str, os.PathLike]

var index : int

var read_envelope : bool

var time_mark : Unit("s")

Methods

def read(self, index: int = None) ‑> CTable

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

Inherited members

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