Module diskchef.maps.radmc_dust

Classes

class RadMCTherm (chemistry: Union[ChemistryBase, PhysicsBase] = None, line_list: List[Line] = None, executable: Union[str, os.PathLike] = 'radmc3d', folder: Union[str, os.PathLike] = 'radmc', radii_bins: Optional[None] = None, theta_bins: Optional[None] = None, outer_radius: Optional[None] = None, verbosity: int = 0, wavelengths: astropy.units.quantity.Quantity = <Quantity [9.12000000e-02, 1.00185062e-01, 1.10055336e-01, 1.20898034e-01, 1.32808958e-01, 1.45893352e-01, 1.60266826e-01, 1.76056380e-01, 1.93401528e-01, 2.12455528e-01, 2.33386735e-01, 2.56380093e-01, 2.81638767e-01, 3.09385935e-01, 3.39866766e-01, 3.73350582e-01, 4.10133237e-01, 4.50539734e-01, 4.94927096e-01, 5.43687520e-01, 5.97251841e-01, 6.56093341e-01, 7.20731931e-01, 7.91738741e-01, 8.69741171e-01, 9.55428434e-01, 1.04955764e+00, 1.15296050e+00, 1.26655065e+00, 1.39133174e+00, 1.52840633e+00, 1.67898555e+00, 1.84439991e+00, 2.02611096e+00, 2.22572425e+00, 2.44500354e+00, 2.68588630e+00, 2.95050094e+00, 3.24118552e+00, 3.56050847e+00, 3.91129125e+00, 4.29663329e+00, 4.71993939e+00, 5.18494979e+00, 5.69577320e+00, 6.25692315e+00, 6.87335783e+00, 7.55052390e+00, 8.29440466e+00, 9.11157287e+00, 1.00092488e+01, 1.09953642e+01, 1.20786320e+01, 1.32686239e+01, 1.45758543e+01, 1.60118735e+01, 1.75893699e+01, 1.93222820e+01, 2.12259213e+01, 2.33171079e+01, 2.56143191e+01, 2.81378525e+01, 3.09100054e+01, 3.39552720e+01, 3.73005596e+01, 4.09754262e+01, 4.50123423e+01, 4.94469770e+01, 5.43185138e+01, 5.96699964e+01, 6.55487093e+01, 7.20065955e+01, 7.91007152e+01, 8.68937507e+01, 9.54545592e+01, 1.04858782e+02, 1.15189513e+02, 1.26538032e+02, 1.39004611e+02, 1.52699404e+02, 1.67743413e+02, 1.84269563e+02, 2.02423878e+02, 2.22366762e+02, 2.44274428e+02, 2.68340447e+02, 2.94777460e+02, 3.23819058e+02, 3.55721846e+02, 3.90767711e+02, 4.29266308e+02, 4.71557804e+02, 5.18015875e+02, 5.69051015e+02, 6.25114158e+02, 6.86700665e+02, 7.54354701e+02, 8.28674040e+02, 9.10315352e+02, 1.00000000e+03] um>, modified_random_walk: bool = True, scattering_mode_max: int = None, nphot_therm: int = None, nphot_mono: int = None, nphot_spec: int = None, coordinate: Union[str, astropy.coordinates.sky_coordinate.SkyCoord] = None, external_source_type: Literal['Draine1978', 'WeingartnerDraine2001', 'None'] = None, external_source_multiplier: float = 1, star_radius: Optional[None] = None, star_effective_temperature: Optional[None] = None, accretion_luminosity: Optional[None] = None, accretion_temperature: Optional[None] = <Quantity 20000. K>)

RadMCTherm(chemistry: Union[diskchef.chemistry.base.ChemistryBase, diskchef.physics.base.PhysicsBase] = None, line_list: List[diskchef.lamda.line.Line] = None, executable: Union[str, os.PathLike] = 'radmc3d', folder: Union[str, os.PathLike] = 'radmc', radii_bins: Optional[int] = None, theta_bins: Optional[int] = None, outer_radius: Optional[astropy.units.quantity.Quantity] = None, verbosity: int = 0, wavelengths: astropy.units.quantity.Quantity = , modified_random_walk: bool = True, scattering_mode_max: int = None, nphot_therm: int = None, nphot_mono: int = None, nphot_spec: int = None, coordinate: Union[str, astropy.coordinates.sky_coordinate.SkyCoord] = None, external_source_type: Literal['Draine1978', 'WeingartnerDraine2001', 'None'] = None, external_source_multiplier: float = 1, star_radius: Optional[astropy.units.quantity.Quantity] = None, star_effective_temperature: Optional[astropy.units.quantity.Quantity] = None, accretion_luminosity: Optional[astropy.units.quantity.Quantity] = None, accretion_temperature: Optional[astropy.units.quantity.Quantity] = )

Expand source code

@dataclass
class RadMCTherm(RadMCBase):
    star_radius: Union[None, u.Quantity] = None
    star_effective_temperature: Union[None, u.Quantity] = None
    accretion_luminosity: Union[None, u.Quantity] = None
    accretion_temperature: Union[None, u.Quantity] = 2e4 * u.K

    def __post_init__(self):
        super().__post_init__()
        self.dust_species = self.chemistry.table.dust_list
        yb08 = diskchef.physics.yorke_bodenheimer.YorkeBodenheimer2008()
        if (self.star_radius is None) and (self.star_effective_temperature is not None):
            self.star_radius = yb08.radius(self.chemistry.physics.star_mass)
            self.logger.warn(
                "Only star effective temperature was given, not star radius! Taking radius from trktime0210")
        elif (self.star_radius is not None) and (self.star_effective_temperature is None):
            self.star_effective_temperature = yb08.effective_temperature(self.chemistry.physics.star_mass)
            self.logger.warn(
                "Only  star radius was given, not effective temperature! Taking effective temperature from trktime0210")
        elif (self.star_radius is None) and (self.star_effective_temperature is None):
            self.star_radius = yb08.radius(self.chemistry.physics.star_mass)
            self.star_effective_temperature = yb08.effective_temperature(self.chemistry.physics.star_mass)
            self.logger.info("Taking effective temperature and radius from trktime0210")

    @property
    def mode(self) -> str:
        """Mode of RadMC: `mctherm`, `image`, `spectrum`, `sed`"""
        return "mctherm"

    def create_files(self) -> None:
        super().create_files()
        self.dustopac()
        self.dust_density()
        self.stars()
        self.external_source()

    def stars(self, out_file: PathLike = None) -> None:
        """Writes the `stars.inp` file"""
        if out_file is None:
            out_file = os.path.join(self.folder, 'stars.inp')
        with open(out_file, 'w') as file:
            if self.accretion_luminosity is None:
                print('2', file=file)  # Typically 2 at present
                print(f'1 {len(self.wavelengths)}', file=file)  # number of stars, wavelengths
                print(f'{self.star_radius.to(u.cm).value} '
                      f'{self.chemistry.physics.star_mass.to(u.g).value} '
                      f'0 0 0',
                      file=file)
                print('\n'.join(f"{entry:.7e}" for entry in self.wavelengths.to_value(u.um)), file=file)
                print(-self.star_effective_temperature.to(u.K).value, file=file)
            else:
                print('2', file=file)  # Typically 2 at present
                print(f'2 {len(self.wavelengths)}', file=file)
                print(f'{self.star_radius.to(u.cm).value} '
                      f'{self.chemistry.physics.star_mass.to(u.g).value} '
                      f'0 0 0',
                      file=file)
                print(f'{self.accretion_effective_radius.to_value(u.cm)} '
                      f'0 '
                      f'0 0 0',
                      file=file)
                print('\n'.join(f"{entry:.7e}" for entry in self.wavelengths.to_value(u.um)), file=file)
                print(-self.star_effective_temperature.to(u.K).value, file=file)
                print(-self.accretion_temperature.to(u.K).value, file=file)

    @property
    def accretion_effective_radius(self):
        return np.sqrt(self.accretion_luminosity / (
                4 * np.pi * c.sigma_sb * self.accretion_temperature ** 4)).cgs

    def dustopac(self, out_file: PathLike = None) -> None:
        """Writes the `dustopac.inp` file"""
        if out_file is None:
            out_file = os.path.join(self.folder, 'dustopac.inp')
        with open(out_file, 'w') as file:
            print('2', file=file)  # Typically 2 at present
            print(len(self.dust_species), file=file)  # Number of species
            for dust_spice in self.dust_species:
                print('-----------------------------', file=file)
                type = "dustkapscatmat_" if os.path.basename(dust_spice.opacity_file).startswith(
                    "dustkapscatmat_") else "dustkappa_"
                print('10' if type == "dustkapscatmat_" else '1', file=file)
                print('0', file=file)
                name_without_whitespaces = re.sub(r"\s*", "", dust_spice.name)
                print(name_without_whitespaces, file=file)
                shutil.copyfile(dust_spice.opacity_file, os.path.join(self.folder, f"{type}{name_without_whitespaces}.inp"))
            print('-----------------------------', file=file)

    def dust_density(self, out_file: PathLike = None) -> None:
        """Writes the dust density file"""
        if out_file is None:
            out_file = os.path.join(self.folder, 'dust_density.inp')
        self.interpolate("Dust density")
        with open(out_file, 'w') as file:
            print('1', file=file)  # Typically 1 at present
            print(self.nrcells, file=file)
            print(len(self.dust_species), file=file)  # Number of species
            for dust_spice in self.dust_species:
                self.interpolate(f"{dust_spice.name} mass fraction")
                print(
                    '\n'.join(
                        f"{entry:.7e}" for entry
                        in (
                                self.polar_table["Dust density"]
                                * self.polar_table[f"{dust_spice.name} mass fraction"]
                        ).to(u.g * u.cm ** (-3)).value),
                    file=file
                )

    @u.quantity_input
    def run(
            self,
            inclination: u.deg = 0 * u.deg, position_angle: u.deg = 0 * u.deg,
            distance: u.pc = 140 * u.pc, velocity_offset: u.km / u.s = 0 * u.km / u.s,
            threads: int = 1, npix: int = 100,
    ) -> None:
        self.logger.info("Running radmc3d")
        start = time.time()
        command = (f"{self.executable} {self.mode} "
                   f"setthreads {threads} "
                   )
        self.logger.info("Running radmc3d for dust temperature calculation: %s", command)
        proc = subprocess.run(
            command,
            cwd=self.folder,
            text=True,
            capture_output=True,
            shell=True
        )
        self.logger.info("radmc3d finished after %s", timedelta(seconds=time.time() - start))
        self.catch_radmc_messages(proc)

    def read_dust_temperature(self) -> None:
        """Read RadMC3D dust temperature output into `table`

        Current limitations:

        * Only one dust population
        """
        self.polar_table["RadMC Dust temperature"] = np.loadtxt(
            os.path.join(self.folder, "dust_temperature.dat")
        )[3:] << u.K
        number_of_zeroes = np.sum(self.polar_table["RadMC Dust temperature"] == 0)
        if number_of_zeroes != 0:
            self.logger.error("RadMC never achieved %d points. "
                              "Recalculate with higher nphot_therm "
                              "and/or remove high-density regions "
                              "out of the model grid "
                              "(highly likely, you can forfeit radii ~< 1 au)."
                              "Values are set to 2.7 K", number_of_zeroes)
            self.polar_table["RadMC Dust temperature"][self.polar_table["RadMC Dust temperature"] == 0] = 2.7 * u.K
        self.interpolate_back("RadMC Dust temperature")
        self.table["Original Dust temperature"] = self.table["Dust temperature"]
        self.table["Dust temperature"] = self.table["RadMC Dust temperature"]
        self.table.check_zeros("Dust temperature")

Ancestors

• RadMCBase
• MapBase

Subclasses

• RadMCThermMono

Class variables

var accretion_luminosity : Optional[None]

var accretion_temperature : Optional[None]

var star_effective_temperature : Optional[None]

var star_radius : Optional[None]

Instance variables

var accretion_effective_radius

Expand source code

@property
def accretion_effective_radius(self):
    return np.sqrt(self.accretion_luminosity / (
            4 * np.pi * c.sigma_sb * self.accretion_temperature ** 4)).cgs

Methods

def dust_density(self, out_file: Union[str, os.PathLike] = None) ‑> None

Writes the dust density file

Expand source code

def dust_density(self, out_file: PathLike = None) -> None:
    """Writes the dust density file"""
    if out_file is None:
        out_file = os.path.join(self.folder, 'dust_density.inp')
    self.interpolate("Dust density")
    with open(out_file, 'w') as file:
        print('1', file=file)  # Typically 1 at present
        print(self.nrcells, file=file)
        print(len(self.dust_species), file=file)  # Number of species
        for dust_spice in self.dust_species:
            self.interpolate(f"{dust_spice.name} mass fraction")
            print(
                '\n'.join(
                    f"{entry:.7e}" for entry
                    in (
                            self.polar_table["Dust density"]
                            * self.polar_table[f"{dust_spice.name} mass fraction"]
                    ).to(u.g * u.cm ** (-3)).value),
                file=file
            )

def dustopac(self, out_file: Union[str, os.PathLike] = None) ‑> None

Writes the dustopac.inp file

Expand source code

def dustopac(self, out_file: PathLike = None) -> None:
    """Writes the `dustopac.inp` file"""
    if out_file is None:
        out_file = os.path.join(self.folder, 'dustopac.inp')
    with open(out_file, 'w') as file:
        print('2', file=file)  # Typically 2 at present
        print(len(self.dust_species), file=file)  # Number of species
        for dust_spice in self.dust_species:
            print('-----------------------------', file=file)
            type = "dustkapscatmat_" if os.path.basename(dust_spice.opacity_file).startswith(
                "dustkapscatmat_") else "dustkappa_"
            print('10' if type == "dustkapscatmat_" else '1', file=file)
            print('0', file=file)
            name_without_whitespaces = re.sub(r"\s*", "", dust_spice.name)
            print(name_without_whitespaces, file=file)
            shutil.copyfile(dust_spice.opacity_file, os.path.join(self.folder, f"{type}{name_without_whitespaces}.inp"))
        print('-----------------------------', file=file)

def read_dust_temperature(self) ‑> None

Read RadMC3D dust temperature output into table

Current limitations:

• Only one dust population

Expand source code

def read_dust_temperature(self) -> None:
    """Read RadMC3D dust temperature output into `table`

    Current limitations:

    * Only one dust population
    """
    self.polar_table["RadMC Dust temperature"] = np.loadtxt(
        os.path.join(self.folder, "dust_temperature.dat")
    )[3:] << u.K
    number_of_zeroes = np.sum(self.polar_table["RadMC Dust temperature"] == 0)
    if number_of_zeroes != 0:
        self.logger.error("RadMC never achieved %d points. "
                          "Recalculate with higher nphot_therm "
                          "and/or remove high-density regions "
                          "out of the model grid "
                          "(highly likely, you can forfeit radii ~< 1 au)."
                          "Values are set to 2.7 K", number_of_zeroes)
        self.polar_table["RadMC Dust temperature"][self.polar_table["RadMC Dust temperature"] == 0] = 2.7 * u.K
    self.interpolate_back("RadMC Dust temperature")
    self.table["Original Dust temperature"] = self.table["Dust temperature"]
    self.table["Dust temperature"] = self.table["RadMC Dust temperature"]
    self.table.check_zeros("Dust temperature")

def stars(self, out_file: Union[str, os.PathLike] = None) ‑> None

Writes the stars.inp file

Expand source code

def stars(self, out_file: PathLike = None) -> None:
    """Writes the `stars.inp` file"""
    if out_file is None:
        out_file = os.path.join(self.folder, 'stars.inp')
    with open(out_file, 'w') as file:
        if self.accretion_luminosity is None:
            print('2', file=file)  # Typically 2 at present
            print(f'1 {len(self.wavelengths)}', file=file)  # number of stars, wavelengths
            print(f'{self.star_radius.to(u.cm).value} '
                  f'{self.chemistry.physics.star_mass.to(u.g).value} '
                  f'0 0 0',
                  file=file)
            print('\n'.join(f"{entry:.7e}" for entry in self.wavelengths.to_value(u.um)), file=file)
            print(-self.star_effective_temperature.to(u.K).value, file=file)
        else:
            print('2', file=file)  # Typically 2 at present
            print(f'2 {len(self.wavelengths)}', file=file)
            print(f'{self.star_radius.to(u.cm).value} '
                  f'{self.chemistry.physics.star_mass.to(u.g).value} '
                  f'0 0 0',
                  file=file)
            print(f'{self.accretion_effective_radius.to_value(u.cm)} '
                  f'0 '
                  f'0 0 0',
                  file=file)
            print('\n'.join(f"{entry:.7e}" for entry in self.wavelengths.to_value(u.um)), file=file)
            print(-self.star_effective_temperature.to(u.K).value, file=file)
            print(-self.accretion_temperature.to(u.K).value, file=file)

Inherited members

• RadMCBase:
  • amr_grid
  • catch_radmc_messages
  • create_files
  • external_source
  • interpolate
  • interpolate_back
  • mode
  • radmc3d
  • run
  • wavelength_micron

class RadMCThermMono (chemistry: Union[ChemistryBase, PhysicsBase] = None, line_list: List[Line] = None, executable: Union[str, os.PathLike] = 'radmc3d', folder: Union[str, os.PathLike] = 'radmc', radii_bins: Optional[None] = None, theta_bins: Optional[None] = None, outer_radius: Optional[None] = None, verbosity: int = 0, wavelengths: astropy.units.quantity.Quantity = <Quantity [1.00000000e-02, 1.09729293e-02, 1.20405177e-02, 1.32119750e-02, 1.44974067e-02, 1.59079019e-02, 1.74556282e-02, 1.91539374e-02, 2.10174801e-02, 2.30623323e-02, 2.53061342e-02, 2.77682421e-02, 3.04698957e-02, 3.34344011e-02, 3.66873319e-02, 4.02567499e-02, 4.41734470e-02, 4.84712111e-02, 5.31871172e-02, 5.83618476e-02, 6.40400427e-02, 7.02706860e-02, 7.71075269e-02, 8.46095441e-02, 9.28414545e-02, 1.01874271e-01, 1.11785918e-01, 1.22661897e-01, 1.34596032e-01, 1.47691275e-01, 1.62060591e-01, 1.77827941e-01, 1.95129342e-01, 2.14114048e-01, 2.34945830e-01, 2.57804398e-01, 2.82886943e-01, 3.10409843e-01, 3.40610526e-01, 3.73749521e-01, 4.10112707e-01, 4.50013774e-01, 4.93796932e-01, 5.41839882e-01, 5.94557071e-01, 6.52403270e-01, 7.15877495e-01, 7.85527313e-01, 8.61953566e-01, 9.45815554e-01, 1.03783672e+00, 1.13881089e+00, 1.24960914e+00, 1.37118727e+00, 1.50459410e+00, 1.65098047e+00, 1.81160919e+00, 1.98786596e+00, 2.18127126e+00, 2.39349353e+00, 2.62636353e+00, 2.88189013e+00, 3.16227766e+00, 3.46994492e+00, 3.80754602e+00, 4.17799333e+00, 4.58448253e+00, 5.03052027e+00, 5.51995432e+00, 6.05700685e+00, 6.64631078e+00, 7.29294983e+00, 8.00250228e+00, 8.78108917e+00, 9.63542705e+00, 1.05728860e+01, 1.16015530e+01, 1.27303021e+01, 1.39688705e+01, 1.53279428e+01, 1.68192432e+01, 1.84556367e+01, 2.02512396e+01, 2.22215421e+01, 2.43835410e+01, 2.67558871e+01, 2.93590457e+01, 3.22154733e+01, 3.53498111e+01, 3.87890977e+01, 4.25630026e+01, 4.67040818e+01, 5.12480588e+01, 5.62341325e+01, 6.17053160e+01, 6.77088069e+01, 7.42963951e+01, 8.15249090e+01, 8.94567062e+01, 9.81602111e+01, 1.07710506e+02, 1.18189976e+02, 1.29689025e+02, 1.42306850e+02, 1.56152301e+02, 1.71344815e+02, 1.88015454e+02, 2.06308029e+02, 2.26380341e+02, 2.48405547e+02, 2.72573651e+02, 2.99093140e+02, 3.28192787e+02, 3.60123625e+02, 3.95161107e+02, 4.33607489e+02, 4.75794431e+02, 5.22085865e+02, 5.72881128e+02, 6.28618411e+02, 6.89778538e+02, 7.56889112e+02, 8.30529071e+02, 9.11333677e+02, 1.00000000e+03] um>, modified_random_walk: bool = True, scattering_mode_max: int = None, nphot_therm: int = None, nphot_mono: int = None, nphot_spec: int = None, coordinate: Union[str, astropy.coordinates.sky_coordinate.SkyCoord] = None, external_source_type: Literal['Draine1978', 'WeingartnerDraine2001', 'None'] = None, external_source_multiplier: float = 1, star_radius: Optional[None] = None, star_effective_temperature: Optional[None] = None, accretion_luminosity: Optional[None] = None, accretion_temperature: Optional[None] = <Quantity 20000. K>, mcmono_wavelengths: astropy.units.quantity.Quantity = <Quantity [0.0912 , 0.09504824, 0.09905885, 0.10323869, 0.10759491, 0.11213493, 0.11686653, 0.12179778, 0.1269371 , 0.13229329, 0.13787548, 0.14369321, 0.14975643, 0.15607548, 0.16266117, 0.16952475, 0.17667795, 0.18413297, 0.19190256, 0.2 ] um>)

Class to run radmc3d mctherm and radmc3d mcmono one after another

Expand source code

@dataclass
class RadMCThermMono(RadMCTherm):
    """Class to run `radmc3d mctherm` and `radmc3d mcmono` one after another"""
    wavelengths: u.Quantity = field(default=np.geomspace(0.01, 1000, 125) * u.um)
    mcmono_wavelengths: u.Quantity = field(default=np.geomspace(0.0912, 0.2, 20) * u.um)

    @u.quantity_input
    def run(
            self,
            inclination: u.deg = 0 * u.deg, position_angle: u.deg = 0 * u.deg,
            distance: u.pc = 140 * u.pc, velocity_offset: u.km / u.s = 0 * u.km / u.s,
            threads: int = 1, npix: int = 100,
    ) -> None:
        super().run(
            inclination=inclination, position_angle=position_angle,
            distance=distance, velocity_offset=velocity_offset,
            threads=threads, npix=npix
        )

        self.logger.info("Running radmc3d")
        start = time.time()
        command = (f"{self.executable} mcmono "
                   f"setthreads {threads} "
                   )
        self.logger.info("Running radmc3d for radiation field: %s", command)
        proc = subprocess.run(
            command,
            cwd=self.folder,
            text=True,
            capture_output=True,
            shell=True
        )
        self.logger.info("radmc3d finished after %s", timedelta(seconds=time.time() - start))
        self.catch_radmc_messages(proc)

    def mcmono_wavelength_micron(self, out_file: PathLike = None) -> None:
        """Creates a `mcmono_wavelength_micron.inp` file"""

        if out_file is None:
            out_file = os.path.join(self.folder, 'mcmono_wavelength_micron.inp')

        with open(out_file, 'w') as file:
            print(len(self.mcmono_wavelengths), file=file)
            print('\n'.join(f"{entry.to(u.um).value:.7e}" for entry in self.mcmono_wavelengths), file=file)

    def create_files(self) -> None:
        super().create_files()
        self.mcmono_wavelength_micron()

    def read_radiation_strength(self) -> None:
        """Read RadMC3D mcmono radiation strength output into `table`"""
        radiation = np.loadtxt(
            os.path.join(self.folder, "mean_intensity.out"), skiprows=4
        ).reshape(-1, len(self.mcmono_wavelengths), order="F") << (u.erg / u.s / u.cm ** 2 / u.Hz / u.sr)

        self.polar_table["Radiation strength"] = np.abs(
            4 * np.pi * u.sr * np.trapz(
                radiation, self.mcmono_wavelengths.to(u.Hz, equivalencies=u.spectral()),
                axis=1)
        ).to(u.erg / u.s / u.cm ** 2)

        number_of_zeroes = np.sum(radiation == 0)
        if number_of_zeroes != 0:
            self.logger.warning("RadMC never achieved %d points. "
                                "Recalculate with higher nphot_therm "
                                "and/or remove high-density regions "
                                "out of the model grid "
                                "(highly likely, you can forfeit radii ~< 1 au). ", number_of_zeroes)

        self.interpolate_back("Radiation strength")
        self.table["Radiation strength"][self.table["Radiation strength"] < 0] = 0
        self.table.check_zeros("Radiation strength")

Ancestors

• RadMCTherm
• RadMCBase
• MapBase

Class variables

var mcmono_wavelengths : astropy.units.quantity.Quantity

var wavelengths : astropy.units.quantity.Quantity

Methods

def mcmono_wavelength_micron(self, out_file: Union[str, os.PathLike] = None) ‑> None

Creates a mcmono_wavelength_micron.inp file

Expand source code

def mcmono_wavelength_micron(self, out_file: PathLike = None) -> None:
    """Creates a `mcmono_wavelength_micron.inp` file"""

    if out_file is None:
        out_file = os.path.join(self.folder, 'mcmono_wavelength_micron.inp')

    with open(out_file, 'w') as file:
        print(len(self.mcmono_wavelengths), file=file)
        print('\n'.join(f"{entry.to(u.um).value:.7e}" for entry in self.mcmono_wavelengths), file=file)

def read_radiation_strength(self) ‑> None

Read RadMC3D mcmono radiation strength output into table

Expand source code

def read_radiation_strength(self) -> None:
    """Read RadMC3D mcmono radiation strength output into `table`"""
    radiation = np.loadtxt(
        os.path.join(self.folder, "mean_intensity.out"), skiprows=4
    ).reshape(-1, len(self.mcmono_wavelengths), order="F") << (u.erg / u.s / u.cm ** 2 / u.Hz / u.sr)

    self.polar_table["Radiation strength"] = np.abs(
        4 * np.pi * u.sr * np.trapz(
            radiation, self.mcmono_wavelengths.to(u.Hz, equivalencies=u.spectral()),
            axis=1)
    ).to(u.erg / u.s / u.cm ** 2)

    number_of_zeroes = np.sum(radiation == 0)
    if number_of_zeroes != 0:
        self.logger.warning("RadMC never achieved %d points. "
                            "Recalculate with higher nphot_therm "
                            "and/or remove high-density regions "
                            "out of the model grid "
                            "(highly likely, you can forfeit radii ~< 1 au). ", number_of_zeroes)

    self.interpolate_back("Radiation strength")
    self.table["Radiation strength"][self.table["Radiation strength"] < 0] = 0
    self.table.check_zeros("Radiation strength")

Inherited members

• RadMCTherm:
  • amr_grid
  • catch_radmc_messages
  • create_files
  • dust_density
  • dustopac
  • external_source
  • interpolate
  • interpolate_back
  • mode
  • radmc3d
  • read_dust_temperature
  • stars
  • wavelength_micron
• RadMCBase:
  • run