Module diskchef.physics.multidust
Module with a definition of multiple dust populations
Classes
class DustPopulation (opacity_file: Union[str, os.PathLike], table: CTable = None, total_dust_density: Unit("g / cm3") = None, dust_temperature: Unit("K") = None, mass_fraction: Unit(dimensionless) = 1.0, name: str = None, mean_mass_per_grain: Unit("g") = None, grain_particle_density: Unit("g / cm3") = <Quantity 3. g / cm3>, opacity_file_format: Literal['radmc'] = 'radmc', average_size_for_chemistry: Unit("cm") = <Quantity 1.e-05 cm>)-
Class with the data on the dust population
Args
total_dust_density- total density of all dust components,
astropy.units.Quantityin density units mass_fraction- multiplier to get the density of this dust population
average_size_for_chemistry- average grain size weighted by surface area, see Vasyunin et al. 2011 Eq. 9
(https://iopscience.iop.org/article/10.1088/0004-637X/727/2/76/pdf) Usage:
>>> # If the name is not set, use index of the created dust (unsafe!) >>> dust0 = DustPopulation("some_opacity.inp", total_dust_density=1e-15 * u.g / u.cm ** 3) >>> dust0.name 'Dust_0' >>> dust1 = DustPopulation("some_opacity.inp", total_dust_density=1e-15 * u.g / u.cm ** 3) >>> dust1.name 'Dust_1'>>> table = CTable() >>> table["Radius"] = [1, 2, 3] * u.au >>> table["Dust density"] = [100, 10, 1] * u.g / u.cm ** 3 * 1e-15 >>> table["Dust temperature"] = [100, 50, 25] * u.K >>> dust = DustPopulation("some_opacity.inp", table=table, name="Default dust") >>> dust.write_to_table() >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Dust temperature Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume AU g / cm3 K 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 ------------ ------------ ---------------- -------------------------- --------------------------- ------------------------ ------------------------------------ 1.000000e+00 1.000000e-13 1.000000e+02 1.000000e+00 7.957747e+00 1.000000e+02 1.000000e-08 2.000000e+00 1.000000e-14 5.000000e+01 1.000000e+00 7.957747e-01 5.000000e+01 1.000000e-09 3.000000e+00 1.000000e-15 2.500000e+01 1.000000e+00 7.957747e-02 2.500000e+01 1.000000e-10 >>> table.dust_list[0].name 'Default dust'>>> dust_large = DustPopulation("large_opacity.inp", table=table, name="Large dust", mass_fraction=0.3, average_size_for_chemistry=1*u.cm) >>> dust_large.write_to_table() >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Dust temperature Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume AU g / cm3 K 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ ---------------- -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- 1.000000e+00 1.000000e-13 1.000000e+02 1.000000e+00 7.957747e+00 1.000000e+02 1.000000e-08 3.000000e-01 2.387324e-15 1.000000e+02 3.000000e-14 2.000000e+00 1.000000e-14 5.000000e+01 1.000000e+00 7.957747e-01 5.000000e+01 1.000000e-09 3.000000e-01 2.387324e-16 5.000000e+01 3.000000e-15 3.000000e+00 1.000000e-15 2.500000e+01 1.000000e+00 7.957747e-02 2.500000e+01 1.000000e-10 3.000000e-01 2.387324e-17 2.500000e+01 3.000000e-16 >>> [dust.name for dust in table.dust_list] ['Default dust', 'Large dust']>>> # Now, the total dust population exceeds the total dust mass (sum of mass_fraction != 1) >>> table.dust_population_fully_set False >>> table.normalize_dust() >>> table.dust_population_fully_set True >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Dust temperature Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume AU g / cm3 K 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ ---------------- -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- 1.000000e+00 1.000000e-13 1.000000e+02 7.692308e-01 6.121344e+00 1.000000e+02 7.692308e-09 2.307692e-01 1.836403e-15 1.000000e+02 2.307692e-14 2.000000e+00 1.000000e-14 5.000000e+01 7.692308e-01 6.121344e-01 5.000000e+01 7.692308e-10 2.307692e-01 1.836403e-16 5.000000e+01 2.307692e-15 3.000000e+00 1.000000e-15 2.500000e+01 7.692308e-01 6.121344e-02 2.500000e+01 7.692308e-11 2.307692e-01 1.836403e-17 2.500000e+01 2.307692e-16>>> # Different ratio at different locations, also note that if Dust temperature was not set, it will be set as np.nan >>> table = CTable() >>> table["Radius"] = [1, 2, 3] * u.au >>> table["Dust density"] = [100, 10, 1] * u.g / u.cm ** 3 * 1e-15 >>> dust = DustPopulation("some_opacity.inp", table=table, name="Default dust", mass_fraction=[0.8, 0.7, 0.5]) >>> dust.write_to_table() >>> dust_large = DustPopulation("large_opacity.inp", table=table, name="Large dust", mass_fraction=[0.2, 0.3, 0.5], average_size_for_chemistry=1*u.cm) >>> dust_large.write_to_table() >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume AU g / cm3 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- 1.000000e+00 1.000000e-13 8.000000e-01 6.366198e+00 nan 8.000000e-09 2.000000e-01 1.591549e-15 nan 2.000000e-14 2.000000e+00 1.000000e-14 7.000000e-01 5.570423e-01 nan 7.000000e-10 3.000000e-01 2.387324e-16 nan 3.000000e-15 3.000000e+00 1.000000e-15 5.000000e-01 3.978874e-02 nan 5.000000e-11 5.000000e-01 3.978874e-17 nan 5.000000e-16 >>> table.dust_population_fully_set True>>> # Rescaling also works if mass fractions are given independently >>> dust_asteroid = DustPopulation("asteroid_opacity.inp", table=table, name="Asteroids", mass_fraction=[0.5, 0.3, 0.5], average_size_for_chemistry=1*u.km) >>> dust_asteroid.write_to_table() >>> table.dust_population_fully_set False >>> table.normalize_dust() >>> table.dust_population_fully_set True >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume Asteroids mass fraction Asteroids number density Asteroids temperature Asteroids surface area per volume AU g / cm3 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- ----------------------- ------------------------ --------------------- --------------------------------- 1.000000e+00 1.000000e-13 5.333333e-01 4.244132e+00 nan 5.333333e-09 1.333333e-01 1.061033e-15 nan 1.333333e-14 3.333333e-01 2.652582e-30 nan 3.333333e-19 2.000000e+00 1.000000e-14 5.384615e-01 4.284941e-01 nan 5.384615e-10 2.307692e-01 1.836403e-16 nan 2.307692e-15 2.307692e-01 1.836403e-31 nan 2.307692e-20 3.000000e+00 1.000000e-15 3.333333e-01 2.652582e-02 nan 3.333333e-11 3.333333e-01 2.652582e-17 nan 3.333333e-16 3.333333e-01 2.652582e-32 nan 3.333333e-21Expand source code
class DustPopulation: """ Class with the data on the dust population Args: total_dust_density: total density of all dust components, `astropy.units.Quantity` in density units mass_fraction: multiplier to get the density of this dust population average_size_for_chemistry: average grain size weighted by surface area, see Vasyunin et al. 2011 Eq. 9 (https://iopscience.iop.org/article/10.1088/0004-637X/727/2/76/pdf) Usage: >>> # If the name is not set, use index of the created dust (unsafe!) >>> dust0 = DustPopulation("some_opacity.inp", total_dust_density=1e-15 * u.g / u.cm ** 3) >>> dust0.name 'Dust_0' >>> dust1 = DustPopulation("some_opacity.inp", total_dust_density=1e-15 * u.g / u.cm ** 3) >>> dust1.name 'Dust_1' >>> table = CTable() >>> table["Radius"] = [1, 2, 3] * u.au >>> table["Dust density"] = [100, 10, 1] * u.g / u.cm ** 3 * 1e-15 >>> table["Dust temperature"] = [100, 50, 25] * u.K >>> dust = DustPopulation("some_opacity.inp", table=table, name="Default dust") >>> dust.write_to_table() >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Dust temperature Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume AU g / cm3 K 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 ------------ ------------ ---------------- -------------------------- --------------------------- ------------------------ ------------------------------------ 1.000000e+00 1.000000e-13 1.000000e+02 1.000000e+00 7.957747e+00 1.000000e+02 1.000000e-08 2.000000e+00 1.000000e-14 5.000000e+01 1.000000e+00 7.957747e-01 5.000000e+01 1.000000e-09 3.000000e+00 1.000000e-15 2.500000e+01 1.000000e+00 7.957747e-02 2.500000e+01 1.000000e-10 >>> table.dust_list[0].name 'Default dust' >>> dust_large = DustPopulation("large_opacity.inp", table=table, name="Large dust", mass_fraction=0.3, average_size_for_chemistry=1*u.cm) >>> dust_large.write_to_table() >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Dust temperature Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume AU g / cm3 K 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ ---------------- -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- 1.000000e+00 1.000000e-13 1.000000e+02 1.000000e+00 7.957747e+00 1.000000e+02 1.000000e-08 3.000000e-01 2.387324e-15 1.000000e+02 3.000000e-14 2.000000e+00 1.000000e-14 5.000000e+01 1.000000e+00 7.957747e-01 5.000000e+01 1.000000e-09 3.000000e-01 2.387324e-16 5.000000e+01 3.000000e-15 3.000000e+00 1.000000e-15 2.500000e+01 1.000000e+00 7.957747e-02 2.500000e+01 1.000000e-10 3.000000e-01 2.387324e-17 2.500000e+01 3.000000e-16 >>> [dust.name for dust in table.dust_list] ['Default dust', 'Large dust'] >>> # Now, the total dust population exceeds the total dust mass (sum of mass_fraction != 1) >>> table.dust_population_fully_set False >>> table.normalize_dust() >>> table.dust_population_fully_set True >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Dust temperature Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume AU g / cm3 K 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ ---------------- -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- 1.000000e+00 1.000000e-13 1.000000e+02 7.692308e-01 6.121344e+00 1.000000e+02 7.692308e-09 2.307692e-01 1.836403e-15 1.000000e+02 2.307692e-14 2.000000e+00 1.000000e-14 5.000000e+01 7.692308e-01 6.121344e-01 5.000000e+01 7.692308e-10 2.307692e-01 1.836403e-16 5.000000e+01 2.307692e-15 3.000000e+00 1.000000e-15 2.500000e+01 7.692308e-01 6.121344e-02 2.500000e+01 7.692308e-11 2.307692e-01 1.836403e-17 2.500000e+01 2.307692e-16 >>> # Different ratio at different locations, also note that if Dust temperature was not set, it will be set as np.nan >>> table = CTable() >>> table["Radius"] = [1, 2, 3] * u.au >>> table["Dust density"] = [100, 10, 1] * u.g / u.cm ** 3 * 1e-15 >>> dust = DustPopulation("some_opacity.inp", table=table, name="Default dust", mass_fraction=[0.8, 0.7, 0.5]) >>> dust.write_to_table() >>> dust_large = DustPopulation("large_opacity.inp", table=table, name="Large dust", mass_fraction=[0.2, 0.3, 0.5], average_size_for_chemistry=1*u.cm) >>> dust_large.write_to_table() >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume AU g / cm3 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- 1.000000e+00 1.000000e-13 8.000000e-01 6.366198e+00 nan 8.000000e-09 2.000000e-01 1.591549e-15 nan 2.000000e-14 2.000000e+00 1.000000e-14 7.000000e-01 5.570423e-01 nan 7.000000e-10 3.000000e-01 2.387324e-16 nan 3.000000e-15 3.000000e+00 1.000000e-15 5.000000e-01 3.978874e-02 nan 5.000000e-11 5.000000e-01 3.978874e-17 nan 5.000000e-16 >>> table.dust_population_fully_set True >>> # Rescaling also works if mass fractions are given independently >>> dust_asteroid = DustPopulation("asteroid_opacity.inp", table=table, name="Asteroids", mass_fraction=[0.5, 0.3, 0.5], average_size_for_chemistry=1*u.km) >>> dust_asteroid.write_to_table() >>> table.dust_population_fully_set False >>> table.normalize_dust() >>> table.dust_population_fully_set True >>> table # doctest: +NORMALIZE_WHITESPACE <CTable length=3> Radius Dust density Default dust mass fraction Default dust number density Default dust temperature Default dust surface area per volume Large dust mass fraction Large dust number density Large dust temperature Large dust surface area per volume Asteroids mass fraction Asteroids number density Asteroids temperature Asteroids surface area per volume AU g / cm3 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm 1 / cm3 K 1 / cm float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 float64 ------------ ------------ -------------------------- --------------------------- ------------------------ ------------------------------------ ------------------------ ------------------------- ---------------------- ---------------------------------- ----------------------- ------------------------ --------------------- --------------------------------- 1.000000e+00 1.000000e-13 5.333333e-01 4.244132e+00 nan 5.333333e-09 1.333333e-01 1.061033e-15 nan 1.333333e-14 3.333333e-01 2.652582e-30 nan 3.333333e-19 2.000000e+00 1.000000e-14 5.384615e-01 4.284941e-01 nan 5.384615e-10 2.307692e-01 1.836403e-16 nan 2.307692e-15 2.307692e-01 1.836403e-31 nan 2.307692e-20 3.000000e+00 1.000000e-15 3.333333e-01 2.652582e-02 nan 3.333333e-11 3.333333e-01 2.652582e-17 nan 3.333333e-16 3.333333e-01 2.652582e-32 nan 3.333333e-21 """ opacity_file: PathLike table: CTable = None total_dust_density: u.g / u.cm ** 3 = None dust_temperature: u.K = None mass_fraction: u.dimensionless_unscaled = 1. name: str = None mean_mass_per_grain: u.g = None grain_particle_density: u.g / u.cm ** 3 = 3 * u.g / u.cm ** 3 opacity_file_format: Literal['radmc'] = 'radmc' average_size_for_chemistry: u.cm = 1e-5 * u.cm _idx = 0 def __post_init__(self): if self.total_dust_density is None: if self.table is None: raise CHEFRuntimeError("table OR total_dust_density must be set for DustPopulation") self.total_dust_density = self.table["Dust density"] if self.mean_mass_per_grain is None: self.mean_mass_per_grain = 4. / 3. * np.pi * self.average_size_for_chemistry ** 3 * self.grain_particle_density if self.name is None: self.name = f"Dust_{self.__class__._idx}" self.__class__._idx += 1 @property @u.quantity_input def number_density(self) -> u.cm ** (-3): """ Number density column of given dust species """ return self.total_dust_density * self.mass_fraction / self.mean_mass_per_grain @property @u.quantity_input def surface_area_per_volume(self) -> u.cm ** (-1): """ Surface area per volume unit, important for chemistry Assumes spherical grains in default setup """ return 4 * np.pi * self.average_size_for_chemistry ** 2 * self.number_density @property @u.quantity_input def temperature(self) -> u.K: if self.dust_temperature is not None: return self.dust_temperature if self.table is not None: if "Dust temperature" in self.table.colnames: return self.table["Dust temperature"] else: return np.nan * u.K def write_to_table(self, table: CTable = None): """ Write the dust population in the table Args: `table`: if given, writes the data in this table. Else, writes in `self.table` Adds columns "{self.name} mass fraction" "{self.name} number density" "{self.name} temperature" "{self.name} surface area per volume" Also adds `self` to `table.dust_list` """ if table is None: table = self.table table[f"{self.name} mass fraction"] = self.mass_fraction table[f"{self.name} number density"] = self.number_density table[f"{self.name} temperature"] = self.temperature table[f"{self.name} surface area per volume"] = self.surface_area_per_volume try: if self not in table.dust_list: table.dust_list.append(self) except KeyError: table.dust_list = [self]Class variables
var average_size_for_chemistry : Unit("cm")var dust_temperature : Unit("K")var grain_particle_density : Unit("g / cm3")var mass_fraction : Unit(dimensionless)var mean_mass_per_grain : Unit("g")var name : strvar opacity_file : Union[str, os.PathLike]var opacity_file_format : Literal['radmc']var table : CTablevar total_dust_density : Unit("g / cm3")
Instance variables
prop number_density : Unit("1 / cm3")-
Number density column of given dust species
Expand source code
@property @u.quantity_input def number_density(self) -> u.cm ** (-3): """ Number density column of given dust species """ return self.total_dust_density * self.mass_fraction / self.mean_mass_per_grain prop surface_area_per_volume : Unit("1 / cm")-
Surface area per volume unit, important for chemistry
Assumes spherical grains in default setup
Expand source code
@property @u.quantity_input def surface_area_per_volume(self) -> u.cm ** (-1): """ Surface area per volume unit, important for chemistry Assumes spherical grains in default setup """ return 4 * np.pi * self.average_size_for_chemistry ** 2 * self.number_density prop temperature : Unit("K")-
Expand source code
@property @u.quantity_input def temperature(self) -> u.K: if self.dust_temperature is not None: return self.dust_temperature if self.table is not None: if "Dust temperature" in self.table.colnames: return self.table["Dust temperature"] else: return np.nan * u.K
Methods
def write_to_table(self, table: CTable = None)-
Write the dust population in the table
Args
table: if given, writes the data in this table. Else, writes inself.tableAdds columns "{self.name} mass fraction" "{self.name} number density" "{self.name} temperature" "{self.name} surface area per volume"Also adds
selftotable.dust_list
class DustPopulationOptool (opacity_file: str, table: CTable = None, dust_density: Unit("g / cm3") = None, dust_temperature: Unit("K") = None, name: str = None, mean_mass_per_grain: Unit("g") = None, grain_particle_density: Unit("g / cm3") = <Quantity 3. g / cm3>, opacity_file_format: Literal['radmc'] = 'radmc', size: Unit("cm") = <Quantity 1.e-05 cm>)-
Class with the data on the dust population from optool
Expand source code
class DustPopulationOptool: """ Class with the data on the dust population from optool """ opacity_file: str table: CTable = None dust_density: u.g / u.cm ** 3 = None dust_temperature: u.K = None name: str = None mean_mass_per_grain: u.g = None grain_particle_density: u.g / u.cm ** 3 = 3 * u.g / u.cm ** 3 opacity_file_format: Literal['radmc'] = 'radmc' size: u.cm = 1e-5 * u.cm _idx = 0 def __post_init__(self): # if self.dust_density is None: # if self.table is None: # raise CHEFRuntimeError("table OR total_dust_density must be set for DustPopulation") # self.dust_density = self.table[f"{self.name} density"] # if self.mean_mass_per_grain is None: # self.mean_mass_per_grain = 4. / 3. * np.pi * self.size ** 3 * self.grain_particle_density if self.name is None: self.name = f"Dust_{self.__class__._idx}" self.__class__._idx += 1 # @property # @u.quantity_input # def number_density(self) -> u.cm ** (-3): # """ # Number density column of given dust species # """ # return self.dust_density / self.mean_mass_per_grain # @property # @u.quantity_input # def surface_area_per_volume(self) -> u.cm ** (-1): # """ # Surface area per volume unit, important for chemistry # Assumes spherical grains in default setup # """ # return 4 * np.pi * self.size ** 2 * self.number_density # @property # @u.quantity_input # def temperature(self) -> u.K: # if self.dust_temperature is not None: # return self.dust_temperature # if self.table is not None: # if f"{self.name} temperature" in self.table.colnames: # return self.table[f"{self.name} temperature"] # else: # return np.nan * u.K def optool_opacity(self, materials=None, mantels=None, porosity=None, porositymantel=None, dhs=None, mmf=None, lammic_min=None, lammic_max=None, lammic_n=None, mie=False, quiet=False, verbose=False, force=False, iformat: int = 1, diana=False, dsharp=False, dsharp_no_ice=False, custom_size=None): """ Taken from DISKLAB by Dullemond and modified Compute the dust opacity for this grain model using the optool code, which is an external code by Carsten Dominik that you can download from: https://github.com/cdominik/optool So please install optool first, before you call the optool_opacity() function. It should be installed such that a simple system call os.system('optool -a 1.0') starts that code and computes (in this example) the opacity for standard settings for a grain size of 1.0 micrometer. Optool will create a file 'dustkappa.dat', which will be moved to dustkappa_1.00e+00.inp and a file dustkappa_1.00e+00.info that contains the optool command used to create it (making it easier to check later if a new opacity has to be calculated). It will then read this file. Note: The grain size is NOT a keyword, because it is a property that should be set upon the initialization of GrainModel. Arguments: materials : list of lists Each 'sublist' contains the optool name of the material and the abundance. For instance: materials = [['pyr',0.870],['c',1,80]] If not set, use optool defaults mantels : list of lists Like materials, but now for the mantel. If not set, no mantel is used porosity : float Value between 0 and 1 for the porosity porositymantel : float Value between 0 and 1 for the porosity of the mantel dhs : float Maximum volume fraction of vacuum in DHS computation mie : bool If True, then switch off DHS and use pure Mie scattering instead mmf : list of floats Use MMF with monom.sz. mmf[0] and frac.dim or fill mmf[1] lammic_min : float Minimum wavelength of wavelength grid in micron lammic_max : float Maximum wavelength of wavelength grid in micron lammic_n : float Number of wavelength points in the wavelength grid The wavelength grid will be logarithmically spaced inherit_density : bool If True, then try to read the material density from the optool opacity file header, and overwrite the self.xigrain with this value. If False, it will merely check if the value from the opacity file equals self.xigrain. quiet : bool If True, then do not display information about the progress. verbose : bool If True, then let optool print more information force : bool If True, then always recompute the opacity. Default is False, meaning it will only recompute if necessary. NOTE: Because this function calls optool from the command line with user-defined options, it may pose a web security risk when used in an online web interface. """ try: q = subprocess.check_output(['which', 'optool']) except: raise RuntimeError( 'Optool not installed. Please get optool from https://github.com/cdominik/optool and install it.') command = 'optool' argmic = self.size.to(u.micron).value if custom_size is None: command += f' -a {argmic:8.2e}' else: command += f' -amin {custom_size[0]:8.2e} -amax {custom_size[1]:8.2e} -apow {custom_size[2]:.2f}' if materials is not None: for m in materials: if m[0] == 'ol-mg100': m[0] = os.path.join(os.path.dirname(os.path.abspath(__file__)), "../dust_opacity/files/ol-mg100-Jaeger2003.lnk") command += ' -c ' + m[0].replace(' ', '') + f' {m[1]:8.2e}' else: if diana: command += ' -diana' elif dsharp: command += ' -dsharp' elif dsharp_no_ice: command += ' -dsharp-no-ice' if mantels is not None: for m in mantels: command += ' -m ' + m[0].replace(' ', '') + f' {m[1]:8.2e}' if porosity is not None: command += f' -p {porosity:8.2e}' if porositymantel is not None: command += f' {porositymantel:8.2e}' if dhs is not None: command += f' -dhs {dhs:8.2e}' if mie: command += ' -mie' if mmf is not None: if type(mmf) is float: command += f' -mmf {mmf:8.2e}' else: command += f' -mmf {mmf[0]:8.2e} {mmf[1]:8.2e}' if lammic_min is not None: command += f' -l {lammic_min:8.2e}' if lammic_max is not None: command += f' {lammic_max:8.2e}' if lammic_n is not None: command += f' {lammic_n:d}' command += f' -radmc {self.name}' if quiet: command += ' -q' if verbose: command += ' -v' filename_base = f'dustkappa_{self.name}' filename_opac = filename_base + '.inp' done = False if force: done = False if not done: if verbose: print('Running ' + command) try: os.system(command) except: raise RuntimeError('Optool failed.') # try: # shutil.move(filename_opac, self.opacity_file) # except: # raise RuntimeError(f'Moving dustkappa_{self.name}.inp failed.') data = ascii.read(filename_opac, data_start=2) data = np.array([data[f'col{i}'] for i in range(1,iformat+2)]).transpose() np.savetxt(self.opacity_file, data, header=f'{iformat}\n{len(data)}', comments='') if filename_opac != self.opacity_file: os.remove(filename_opac) else: if (verbose or not quiet): print('Opacity ' + filename_base + ' already available. Not computing new. If you want to force recompute: set force=True') def dhs_opacity(self, materials, wav_grid=None): argmic = self.size.to(u.micron).value nmat = len(materials) mass_norm = 0 if wav_grid==None: wmin = np.zeros(nmat) wmax = np.zeros(nmat) nwav = np.zeros(nmat) for im in range(len(materials)): m = materials[im] filenames = glob(os.path.join(os.path.dirname(os.path.abspath(__file__)), "../dust_opacity/files/newDHS", f"*{m[0]}*rv{argmic:.1f}.dat")) inp = np.loadtxt(filenames[0], max_rows=1) nwav[im] = inp[0] data = np.loadtxt(filenames[0], skiprows=1).transpose() wmin[im] = np.min(data[0,:]) wmax[im] = np.max(data[0,:]) if (np.max(wmin) == np.min(wmin)) & (np.max(wmax) == np.min(wmax)) & (np.max(nwav) == np.min(nwav)): wav_grid = data[0,:] else: wav_grid = np.geomspace(np.max(wmin), np.min(wmax), int(np.min(nwav))) kappa_abs = np.zeros(len(wav_grid)) for im in range(len(materials)): m = materials[im] filenames = glob(os.path.join(os.path.dirname(os.path.abspath(__file__)), "../dust_opacity/files/newDHS", f"*{m[0]}*rv{argmic:.1f}.dat")) inp = np.loadtxt(filenames[0], max_rows=1) data = np.loadtxt(filenames[0], skiprows=1).transpose() kappa = data[1,:]*3/(4*argmic*1e-4*inp[2]) * m[1] mass_norm += m[1] kappa_abs += np.interp(wav_grid, data[0,:], kappa) kappa_abs /= mass_norm np.savetxt(self.opacity_file, np.array([wav_grid, kappa_abs]).transpose(), header=f'{1}\n{len(wav_grid)}', comments='')Class variables
var dust_density : Unit("g / cm3")var dust_temperature : Unit("K")var grain_particle_density : Unit("g / cm3")var mean_mass_per_grain : Unit("g")var name : strvar opacity_file : strvar opacity_file_format : Literal['radmc']var size : Unit("cm")var table : CTable
Methods
def dhs_opacity(self, materials, wav_grid=None)def optool_opacity(self, materials=None, mantels=None, porosity=None, porositymantel=None, dhs=None, mmf=None, lammic_min=None, lammic_max=None, lammic_n=None, mie=False, quiet=False, verbose=False, force=False, iformat: int = 1, diana=False, dsharp=False, dsharp_no_ice=False, custom_size=None)-
Taken from DISKLAB by Dullemond and modified
Compute the dust opacity for this grain model using the optool code, which is an external code by Carsten Dominik that you can download from:
https://github.com/cdominik/optool
So please install optool first, before you call the optool_opacity() function. It should be installed such that a simple system call os.system('optool -a 1.0') starts that code and computes (in this example) the opacity for standard settings for a grain size of 1.0 micrometer. Optool will create a file 'dustkappa.dat', which will be moved to dustkappa_1.00e+00.inp and a file dustkappa_1.00e+00.info that contains the optool command used to create it (making it easier to check later if a new opacity has to be calculated). It will then read this file.
Note: The grain size is NOT a keyword, because it is a property that should be set upon the initialization of GrainModel.
Arguments
materials : list of lists Each 'sublist' contains the optool name of the material and the abundance. For instance: materials = [['pyr',0.870],['c',1,80]] If not set, use optool defaults
mantels : list of lists Like materials, but now for the mantel. If not set, no mantel is used
porosity : float Value between 0 and 1 for the porosity
porositymantel : float Value between 0 and 1 for the porosity of the mantel
dhs : float Maximum volume fraction of vacuum in DHS computation
mie : bool If True, then switch off DHS and use pure Mie scattering instead
mmf : list of floats Use MMF with monom.sz. mmf[0] and frac.dim or fill mmf[1]
lammic_min : float Minimum wavelength of wavelength grid in micron
lammic_max : float Maximum wavelength of wavelength grid in micron
lammic_n : float Number of wavelength points in the wavelength grid The wavelength grid will be logarithmically spaced
inherit_density : bool If True, then try to read the material density from the optool opacity file header, and overwrite the self.xigrain with this value. If False, it will merely check if the value from the opacity file equals self.xigrain.
quiet : bool If True, then do not display information about the progress.
verbose : bool If True, then let optool print more information
force : bool If True, then always recompute the opacity. Default is False, meaning it will only recompute if necessary.
NOTE: Because this function calls optool from the command line with user-defined options, it may pose a web security risk when used in an online web interface.