Read routines for a model

This module provides several routines to read the output of a ProDiMo model. All the data belonging to a ProDiMo model is put into an hierarchical data structure Data_ProDiMo.

The module provides a routine to read “all” the data of a ProDiMo model and specialized routines to read only distinct model data (e.g. only the SED of a model).

The following example reads the output files of a ProDiMo model from the current working directory. The read_prodimo() function tries to read (nearly)all ProDiMo output data which it can find, and returns the a prodimopy.read.Data_ProDiMo object, here called model.

There are also more specialized reading routines in the read module. The routines have read_ as a name prefix.

import prodimopy.read as pread

model=pread.read_prodimo()

print(model)

class prodimopy.read.DataStar(nlam, teff, r, logg, luv)[source]

Bases: object

Stellar spectrum and several stellar properties.

This object holds data from StarSpectrum.out but also from other sources such as ProDiMo.out.

It also provides some convenience functions for e.g. Lstar or the mass accretion rate derived from the UV luminosity.

Constructor for reading StarSpectrum.out. Other properties are filled later, if possible.

nlam: int: number of wavelength points

teff: float: effective temperature Unit: \(\mathrm{K}\)

r: float: radius Unit: \(\mathrm{R_\odot}\)

logg: float: surface gravity Unit: \(\mathrm{cm/s^2}\)

lam: ndarray[tuple[int, ...], dtype[float64]]: Stellar spectrum wavelength Unit: \(\mathrm{\mu m}\)

nu: ndarray[tuple[int, ...], dtype[float64]]: Stellar spectrum frequency Unit: \(\mathrm{Hz}\)

Inu: ndarray[tuple[int, ...], dtype[float64]]: Stellar Spectrum intensity Unit: \(\mathrm{erg/cm^2/s/Hz/sr}\)

InuBand: ndarray[tuple[int, ...], dtype[float64]] | None: Stellar Spectrum band averaged (from RTinterpolation.out) Unit: \(\mathrm{erg/cm^2/s/Hz/sr}\)

InuInterpol: ndarray[tuple[int, ...], dtype[float64]] | None: Stellar Spectrum interpolated (J_NU_FROM_RT, from RTinterpolation.out) Unit: \(\mathrm{erg/cm^2/s/Hz/sr}\)

wlInterpol: ndarray[tuple[int, ...], dtype[float64]] | None: wl grid for the interpolated spectrum is in micron.

mass: float | None: stellar mass Unit: \(\mathrm{M_\odot}\)

Lstar: float | None: stellar luminosity (from ProDiMo.out) Unit: \(\mathrm{L_\odot}\). Does not include accretion luminosity.

property lumAccr

The accretion luminosity Unit: \(\mathrm{erg/s}\)

Is the excess luminosity \(L_\mathrm{accr}\) due to accretion, i.e. the part that is added to the pure stellar photosphere spectrum (i.e. no accretion). We consider here only the part of the spectrum with \(\lambda \geq 0.0912\) micron.

\[L_\mathrm{accr}=L(\lambda \geq 0.0912\,\mathrm{micron}) - L_\mathrm{star}\]

Where \(L_\mathrm{star}\) is the value from Parameter.in (read from ProDiMo.out), and \(L\) is the total luminosity for \(\lambda \geq 0.0912\) micron.

property mdot

The mass accretion rate in \(\mathrm{M_\odot/yr}\) derived from the UV luminosity.

Uses:

\[L_\mathrm{accr} = G \times M_{*} \times \dot{M}_\mathrm{accr} / R\]

where we use lumAccr for \(L_\mathrm{accr}\).

property lumTot: The total stellar luminosity (full spectrum) Unit: erg/s

property lum

The “stellar” luminosity \((\lambda \geq 0.0912 \, \mathrm{\mu m})\) Unit: \(\mathrm{erg/s}\)

This includes the accretion component if present.

property lumEUV

stellar EUV luminosity Unit: erg/s

Defined as everything with energies between and 13.6 eV and 100 eV.

property lumXray

stellar X-ray luminosity Unit: erg/s

The luminosity for the part of the spectrum with energies >= 0.1 keV.

calcLum(lammin=None, lammax=None)[source]

Integrates the stellar spectrum from lammin to lammax and returns the luminosity.

Parameters:

lammin (float | None) – minimum wavelength in micron. Default: None (i.e. first wavelength point)
lammax (float | None) – maximum wavelength in micron. Default: None (i.e. last wavelength point)

Returns:

The luminosity from lammin to lammax in erg/s

Return type:

float

class prodimopy.read.DataContinuumObs(nlam=None)[source]

Bases: object

Holds the observational data for the continuum (the dust).

Holds the photometric data, spectra (experimental) and radial profiles (experimental).

Todo

proper documentation

specs: holds a list of spectra if available (wl,flux,error)

class prodimopy.read.DataContinuumImages(nlam, nx, ny, incl=None, filepath='./image.out')[source]

Bases: object

Holds the continuum images (image.out) and provides a method to read one particular image.

The coordinates x,y are the same for all images, and only stored once.

Todo

find a better solution for the multiple inclinations (Problem: don’t want to read through the whole file).

nlam: int: The number of wavelength points (== number of images)

lams: ndarray[tuple[int, ...], dtype[float64]]: The wavelengths. UNIT: micron, DIMS: (nlam)

nx: int: The number of x axis points (radial) of the image

ny: int: The number of y axis points (or theta) of the image

x: ndarray[tuple[int, ...], dtype[float64]]: x coordinates. UNIT: au, DIMS: (nx,ny)

y: ndarray[tuple[int, ...], dtype[float64]]: y coordinates. UNIT: au, DIMS: (nx,ny)

property inclination

getImage(wl, iinc=0)[source]

Reads the intensities at a certain wavelength (image) from the image.out file.

The image with the closest wavelength to the given wavelength (wl) will be returned.

Parameters:

wl (float) – the wavelength in micron of the requested image
iinc (int) – the inclination index in case of multiple inclinations. Default iinc=0

Returns:

x, y coordinates DIM(nx,ny), intensities DIM(nx,ny), UNIT: erg/cm^2/s/Hz/sr, and the exact wavelength Unit: micron

Return type:

tuple

class prodimopy.read.DataSEDAna(nlam, nx)[source]

Bases: object

Holds the analysis data for the Spectral Energy Distribution (SEDana.out).

Initializes the SED analysis data.

Parameters:

nlam (int) – The number of wavelength points.
nx (int) – The number of radial points.

lams: Wavelength points shape=(nlam)

r15: The r15 radius (when 15% of the emission are reached. shape=(nlam)

r85: The r85 radius (when 85% of the emission are reached. shape=(nlam)

z15: The z15 height for each radial grid point (when 15% of the emission are reached from the top). shape=(nlam,nx)

z85: The z85 height for each radial grid point (when 85% of the emission are reached from the top). shape=(nlam,nx)

class prodimopy.read.DataSED(nlam, distance)[source]

Bases: object

Holds the data for the Spectral Energy Distribution (SED).

In case of multiple inclinations the data for all inclinations is stored. One can select a SED for a given inclination via sedinc()

Parameters:

nlam (int) – The number of wavelength points.
distance (float) – The distance of the target in pc.

nlam: int: number of wavelength points

distance: float: distance for which the SED was calculated. UNIT: pc

lam: ndarray[tuple[int, ...], dtype[float64]]: the wavelength points in micron DIM: nlam

nu: ndarray[tuple[int, ...], dtype[float64]]: the frequency points in Hz DIM: nlam

sedAna: DataSEDAna | None: the SED analysis data. Holds the data for the emission origin

fnuFaceonErg: ndarray[tuple[int, ...], dtype[float64]]: the flux density in erg/cm^2/Hz/s of the face-on SED (i.e. not from raytracing)

property fnuErg: float64: The flux density F_nu in erg/s/cm^2/Hz, without extinction.

property nuFnuW: float64: The flux density nuF_nu in W/m^2; without extinction.

property fnuJy: float64: The flux density in Jy; without extinction.

property inclination: float64: The current inclination. See prodimopy.read.sedinc().

property Lbol: float64: The bolometric luminosity in L_sun. See setLbolTbol().

property Tbol: float64: The bolometric temperature in K. See setLbolTbol().

setLbolTbol()[source]

Calculates the bolometric temperature and luminosity for the given values of the SED. Quite approximate at the moment.

Todo

check the correctness of the calculation
provide a proper reference for the procedure we follow.
allow for more flexibility

extinction(sedObs, extmodel='F99')[source]

Return the extinction for the given model.

Note

Please cite https://joss.theoj.org/papers/10.21105/joss.07023 if you use this function in your work.

Parameters:

sedObs (DataContinuumObs) – The observed data for the continuum. Needs the R_V and E_BV values.
extmodel (str) – The model to use for the extinction. Default: “F99” For supported models see astropy dust_extinction

Return type:

The extinction values for the given model. shape=(nlam) E.g. do flux*extinction(model.sedObs).

class prodimopy.read.DataLineEstimateRInfo(ix, iz, Fcolumn, tauLine, tauDust, z15, z85, ztauD1, Ncol)[source]

Bases: object

Data container for the extra radial information for a single line estimate. This object corresponds to one line of the radial information in FlineEstimates.out, i.e. it represents the vertical column at grid point ix.

ix: int: The x-coordinate index.

iz: int: The z-coordinate index.

Fcolumn: float: the line flux for the column (check again, i guess the iz coordinate is irrelevant for that); UNIT: W/m^2

tauLine: float: the total vertical optical depth of the line measured from tauDust=1 upwards f(r)

tauDust: float: the total vertical optical depth of the dust at the wl of the line as f(r)

z15: float: z-level where the line flux reaches 15% of the total flux as f(r) (integrated from top to bottom of the disk); UNIT: au

z85: float: z-level where the line flux reaches 85% as f(r) (integrated from top to bottom of the disk); UNIT: au

ztauD1: float: z-level where taudust_ver(lambda_line)=1; UNIT: au Is None for FlineEstimates version < 3

Ncol: float: The column density of the species as f(r) measured from ztauD1 upwards; UNIT: cm^-2 It considers both halves of the disks (i.e. for the case of tauDust <1 ) so the values can be different to cdnmol where only one half of the disk is considered. Is None for FlineEstimates version < 3

class prodimopy.read.DataLineEstimate(ident, wl, jup, jlow, flux)[source]

Bases: object

Data container for a single line estimate.

This includes also the radial information (e.g. line emitting region). The radial information is only read and put into memory when accessed.

Parameters:

ident (str)
wl (float)
jup (int)
jlow (int)
flux (float)

ident: str: The line identification (species name)

wl: float: The wavelength of the line in micron

jup: int: Upper level as defined in ProDiMo (not necessarily the real one!)

jlow: int: Lower level as defined in ProDiMo (not necessarily the real one!).

flux: float: The integrated line flux in W/m^2

property frequency: float: Frequency of the line in GHz.

property flux_Jy: float: The flux of the line in [Jansky km s^-1].

property rInfo: list[DataLineEstimateRInfo]

The radial information for the line estimate (e.g. used for line origin plots). It used lazy initialisation, e.g. if the rInfo object is accessed the firs time the data is read from FlineEstimates.out.

Returns:: A list with nx entries.
Return type:: list[DataLineEstimateRInfo]

property rInfo_Fcolumn: ndarray: The line flux per column as a numpy array, one entry per radial grid point (see Fcolumn). DIM: (nx)

property rInfo_tauLine: ndarray: The total vertical optical depth of the line per radial grid point as a numpy array (see tauLine). DIM: (nx)

property rInfo_tauDust: ndarray: The total vertical optical depth of the dust at the line wavelength per radial grid point as a numpy array (see tauDust). DIM: (nx)

property rInfo_Ncol: ndarray: The column density of the species per radial grid point as a numpy array (see Ncol). DIM: (nx)

luminosityErgs(distance)[source]

Returns the luminosity in erg/s for the given distance

simply calls luminosityWatt() and converts the result to erg/s.

Parameters:: distance (float)
Return type:: float

luminosityWatt(distance)[source]

Parameters:: distance (float) – The distance in pc, to calculate the line luminosity.
Returns:: Returns the line luminosity in Watt
Return type:: float

class prodimopy.read.DataLineProfile(nvelo, restfrequency=None)[source]

Bases: object

Data container for a spectral line profile for a single spectral line.

nvelo: int: number of velocity points of the profile.

restfrequency: float: The rest frequency of the line. Useful for conversions. Optional. UNIT: GHz

velo: ndarray[tuple[int, ...], dtype[float64]]: The velocity grid of the line profile. UNIT: kms/s, DIMS: (nvelo)

property flux: ndarray[tuple[int, ...], dtype[float64]]: The flux at each velocity point. UNIT: given by flux_unit, DIMS: (nvelo)

property flux_dust: ndarray[tuple[int, ...], dtype[float64]]: The flux of the dust (continuum) at each velocity point. UNIT: given by flux_unit, DIMS: (nvelo)

property flux_conv: ndarray[tuple[int, ...], dtype[float64]]: The flux at each velocity point for the convolved spectrum. UNIT: given by flux_unit, DIMS: (nvelo)

property frequency: float: The frequencies according to the velocity grid. Is relative to the rest frequency. UNIT: GHz FIXME: restfreq can be None …

property fluxErgAng: The fluxes in units of erg/s/cm^2/Angstrom.

Deprecated since version 2.4.0: Use flux_unit instead.

property flux_unit: str: The flux unit for the line profile (e.g. for plotting). Possible choices: Jy (default), mJy and ErgAng.

convolve(R)[source]

Convolves the given profile to the spectral resolution R. The profile is convolved with a Gaussian where R determines the FWHM of that Gaussian. The results are stored in flux_conv

Parameters:: R – The spectral resolution.
Returns:: The FWHM of the Gaussian in units of km/s i.e. the spectral resolution.
Return type:: float

class prodimopy.read.DataLine[source]

Bases: object

Data container for a single spectral line read from line_flux.out.

All lines can then be accessed via model.lines, where model is a Data_ProDiMo object. To print a decent looking list of all lines (in the notebook or console) one can use the following code:

print(*model.lines,sep="\n")

This works also with a list of lines selected via selectLines():. But of course a simple for loop will also do the trick.

wl: float: Rest-wavelength of the line. UNIT: micron

frequency: float: Rest-frequency of the line. UNIT: GHz

prodimoInf: str: Some information String from ProDiMo (i.e. transition)

species: str: The chemical species of the line. Derived from ident.

ident: str: The identification of the line (i.e. as given in LineTransferList.in)

distance: float | None: The distance used for the line calculations. UNIT: pc

property flux: float: The integrated line flux for the current inclination. UNIT: W/m^2

property fcont: float: The continuum flux at the line position for the current inclination. UNIT: Jy

property profile: list[DataLineProfile]: The line profile(s) for this particular line for each inclination.

property inclination: float: The current inclination. Unit: degrees.

property flux_Jy: float: The flux of the line in Jansky km s^-1.

property fcontErgAng: float: The flux of the continuum in erg/s/cm^2/Angstrom

property luminosityErgs: Returns the luminosity in erg/s

property luminosityWatt: Returns the line luminosity in Watt

property FWHM: float

Calculate the Full Width at Half Maximum (FWHM) of a flux profile.

This method computes the FWHM by finding the velocity range where the flux is above half of the maximum value. The calculation searches for the leftmost and rightmost points where the flux crosses the half-maximum threshold.

Parameters:: convolved (bool) – If True, the FWHM is calculated from the convolved flux profile.
Returns:: The FWHM in km/s, calculated as the difference between the right and left velocity values at half maximum (vright - vleft).
Return type:: float

Notes

The half maximum is calculated as the average of the maximum and minimum flux values: 0.5 * (max(flux) + min(flux)) , so the continuum is taken care off
The left boundary search starts from index 0 and stops when flux exceeds half maximum or reaches the end of the array
The right boundary search starts from the last index and stops when flux exceeds half maximum, or reaches index 0
The method assumes that vright > vleft in the velocity profile

property FWHM_convolved: float: Calculate the Full Width at Half Maximum (FWHM) of a convolved flux profile. For details see FWHM()

convolve(R)[source]

Convolves all the profiles (each inclination) for this line. simply calls DataLineProfile.convolve()

Parameters:: R (float) – The desired spectral resolution.
Returns:: The FWHM of the Gaussian in units of km/s i.e. the spectral resolution.
Return type:: float

class prodimopy.read.ModelParameters(filename='Parameter.out', directory='.')[source]

Bases: MutableMapping

Access to the Parameters of a model. Reads Parameter.out and provides some utility functions fo access Parameters.

Opens the file and fills a dictionary with all Parameters. The f90nml package is used, so type conversion is properly done. This dictionary is readonly.

The get routine has some special treatment for certain parameters (e.g. return a list with the correct len). However, special treatment for some fields might be missing.

Parameters:

filename (string) – The Filename of the Parameter file. Default: Parameter.out
directory (string) – The directory that contains the ProDiMo model output. Default: “.”
..todo:: –
- make some utility routine, mainly unit conversion, for fields like Mdisk etc.

mapping: nml.NameList

class prodimopy.read.DataFLiTsSpec[source]

Bases: object

Data container for a FLiTs spectrum. Currently provides only the wavelength grid and the flux in Jy, as produced by FLiTs.

wl: ndarray[tuple[int], float64]: Wavelength grid in micron.

flux: ndarray[tuple[int], float64]: Flux in Jy.

flux_cont: ndarray[tuple[int], float64]: Continuum flux of the spectrum in Jy.

conv_wl: ndarray[tuple[int], float64] | None: Convolved Wavelength grid in micron.

conv_flux: ndarray[tuple[int], float64] | None: Convolved Flux in Jy.

conv_flux_cont: ndarray[tuple[int], float64] | None: Convolved Continuum flux of the spectrum in Jy.

conv_R: ndarray[tuple[int, ...], dtype[float64]] | None: Convolved spectrum. Unit: Jy.

convolve(specR, sample=1, contReturn=False, inplace=False)[source]

Returns a convolved version of the Spectrum. Does not change the original FLiTs spectrum.

Parameters:: specR (int) – The desired spectral resolution.
Returns:: tuple containing: wls(array_like): array of wavelength points in micron, flux(array_like): flux values for each wavelength point in Jansky.
Return type:: (tuple)

Todo

allow for different units.

class prodimopy.read.DataGas(nlam)[source]

Bases: object

Holds the data for the gas (mainly from dust_opac.out).

Todo

currently only the init cs is read, not the x,y data

class prodimopy.read.DataDust(amin, amax, apow, nsize, nlam)[source]

Bases: object

Holds the data for the dust (mainly from dust_opac.out)

Todo

dust composition is not yet read.

amin: float: The minimum grain size of the size distribution. UNIT: micron

amax: float: The maximum grain size of the size distribution. UNIT: micron

apow: float: The power-law exponent of the grain size distribution.

nsize: int: The number for grain size bins

lam: ndarray[tuple[int, ...], dtype[float64]]: The wavelength grid for the dust opacities. UNIT: micron DIM: (nlam)

energy: ndarray[tuple[int, ...], dtype[float64]]: The energy grid for the dust opacities (for convenience). UNIT: eV DIM: (nlam)

kext: ndarray[tuple[int, ...], dtype[float64]]: The dust extinction coefficient for each wavelength per g of dust. UNIT: cm² g^-1 DIM: (nlam)

kabs: ndarray[tuple[int, ...], dtype[float64]]: The dust absorption coefficient for each wavelength per g of dust. UNIT: cm² g^-1 DIM: (nlam)

ksca: ndarray[tuple[int, ...], dtype[float64]]: The dust scattering coefficient for each wavelength per g of dust. UNIT: cm² g^-1 DIM: (nlam)

ksca_an: ndarray[tuple[int, ...], dtype[float64]]: The dust anisotropic (approximation) scattering coefficient for each wavelength per g of dust. UNIT: cm² g^-1 DIM: (nlam)

asize: ndarray[tuple[int, ...], dtype[float64]] | None: The grain size for each size bin. UNIT: micron DIM: (nsize)

sigmaa: ndarray[tuple[int, ...], dtype[float64]] | None: The radial surface density profiles for each individual grain. size UNIT: g cm^-2 DIM: (nsize,nr)

fsize: ndarray[tuple[int, ...], dtype[float64]] | None: The normalized dust size distribution function. DIM: (nsize,nr,nz)

fsize_rho: ndarray[tuple[int, ...], dtype[float64]] | None: The dust size distribution function. Density for each grain size bin. UNIT: UNIT: g cm^-3 DIM: (nsize,nr,nz)

class prodimopy.read.DataElements[source]

Bases: object

Data for the Element abundances (the input). Holds the data from Elements.out. The data is stored as OrderedDict (same order as in Elements.out) with the element names as keys.

abun: OrderedDict[str, float]: Ordered dictionary holding the element abundances relative to hydrogen.

abun12: OrderedDict[str, float]: abundances in the +12 notation

massRatio: OrderedDict[str, float]: mass ratios

amu: OrderedDict[str, float]: atomic mass unit

muHamu: float | None: rho = muH*n<H> with muH/amu from Elements.out

class prodimopy.read.DataSpecies[source]

Bases: object

Data for the Species (the input).

Holds the data from Species.out.

The data is stored as OrderedDict with the element names as keys.

mass: OrderedDict[str, float]: The mass of the species UNIT: g.

charge: OrderedDict[str, int]: The charge of the species UNIT: e.

chemPot: OrderedDict[str, float]: The chemical potential as determined by ProDiMo.

get_spdata(name)[source]

Returns all data of the species with name.

Parameters:: name (str or int) – The name of the species. if int then use it as index
Return type:: (mass,charge,chemPot)

class prodimopy.read.DataBgSpec(nlam)[source]

Bases: object

Background field input spectrum.

Todo

is not yet included in Data_ProDiMo

class prodimopy.read.Chemistry(name)[source]

Bases: object

Data container for chemistry analysis. Holds the information for one particular molecule.

Also provides some convenience functions for further analysis.

Parameters:: name (string) – The name of the model (can be empty).

name: str: The name of the model (can be empty)

totfrate: ndarray[tuple[int, ...], dtype[_ScalarType_co]]: Total formation rate at each spatial grid point.

totdrate: ndarray[tuple[int, ...], dtype[_ScalarType_co]]: Total destruction rate at each spatial grid point.

gridf: ndarray[tuple[int, ...], dtype[_ScalarType_co]]: Contains for each individual grid point a sorted list of the formation reactions including the index (0) and the rate (1)

gridd: ndarray[tuple[int, ...], dtype[_ScalarType_co]]: Contains for each individual grid point a sorted list of the destruction reactions including the index (0) and the rate (1)

fridxs: list: List of all formation reaction indices for this species.

dridxs: list: List of all destruction reaction indices for this species.

species: str: name of species for which the chamanalysis should be done

model: Data_ProDiMo: A reference to the model, for which we do the analysis.

reac_rates_ix_iz(ix=None, iz=None, locau=None, lowRatesFrac=0.001)[source]

Function that analyse the Chemistry manually via a point-by-point analysis for a given species. Shows the most important formation and destruction rates for the given point ix,iz (or in au via Parameter locau).

Parameters:

ix (x-index corresponding to desired radial location in grid, starting at 0)
iz (z-index corresponding to desired radial location in grid, starting at 0)
locau (tuple[float, float]) – the desired coordinates in au (x,z). The routine then finds the closest grid point for those coordinates.
lowRatesFrac (float) – Only rates with rate/total_rate > lowRatesFrac are printed. Default: 1.e-3 This is useful to avoid printing all the reactions that are not important.

reac_to_str(reac, idx, rate=None)[source]

Converts a Reaction to the output format we want for the chemanalysis.

Parameters:: reac (prodimopy.chemistry.network.Reaction)

get_reac_grid(level, rtype)[source]

Gets the counting index of the x important Reaction. If it does not exist it is set to 0.

Parameters:

level (int) – The level of importance 1 means most important, 2 second most important etc.
type (str) – Formation (pass f) or destruction Reaction (pass d)

Returns:

And array of shape (model.nx,model.nz) where the first one contains the index of the Reactions for the list of reactions from chemanalysis (the output file) and the second one the rate itself (for this Reaction at each point).

Return type:

tuple

class prodimopy.read.DataLineObs(flux, flux_err, fwhm, fwhm_err, flag)[source]

Bases: DataLine

Holds the observational data for one line.

Differently to DataLine DataLineObs does not care about multiple inclinations. So simply some setter and getters are used to fill the list quantities (first entry).

property flux: The integrated line flux for the current inclination. UNIT: W/m^2

property profile: The line profile(s) for this particular line for each inclination.

class prodimopy.read.Data_ProDiMo(name)[source]

Bases: object

Data container for most of the output produced by ProDiMo.

The class also includes some convenience functions and derives/calculates some additional quantities not directly included in the ProDiMo output.

Todo

maybe put all the (real) observational data into one object (e.g. also the lines) e.g.DataObservations
maybe also have something like DataObservables (hold, SED, imaged, lines etc.)

Parameters:: name (string) – The name of the model (can be empty).

name: str: The name of the model (can be empty).

directory: str | None: The directory from which the model was read. Is e.g. set by read_prodimo(). Can be a relative path.

params: ModelParameters: Dictionary that allows to access the models Parameters from Parameter.out.

nx: int: The number of spatial grid points in the x (radial) direction

nz: int: The number of spatial grid points in the z (vertical) direction

nspec: int | None: The number of chemical species included in the model.

nheat: int: The number of heating processes included in the model.

ncool: int: The number of cooling processes included in the model.

p_dust_to_gas: float: The global dust to gas mass ratio (single value, given Parameter).

p_v_turb: float: The global turbulent velocity (single value). UNIT: km s^-1

p_rho_grain: float: The global grain mass density (the density of one dust grain). UNIT: g cm^-3.

p_mdisk: float: The disk mass in solar units. Is taken from ProDiMo.out. Represents the parameter value. Please note if also an envelope is included or the structure comes from a 1D or 2D interface this value is not the actually disk mass.

td_fileIdx: int | None: The file index (prefix) in case of a time-dependent chemistry model (Starts at 1).

age: float: The age of the model in years. Is taken from Parameter.out. np.inf for steady-state model. UNIT: years.

x: ndarray[tuple[int, ...], dtype[float64]]: The x coordinates (radial direction). UNIT: au, DIMS: (nx,nz)

z: ndarray[tuple[int, ...], dtype[float64]]: The z coordinates (vertical direction). UNIT: au, DIMS: (nx,nz)

rhog: ndarray[tuple[int, int], float64]: The gas density. UNIT: g cm^-3, DIMS: (nx,nz)

rhod: ndarray[tuple[int, int], float64]: The dust density. UNIT: g cm^-3, DIMS: (nx,nz)

d2g: ndarray[tuple[int, int], float64]: The dust to gas mass ratio from ProDiMo UNIT: , DIMS: (nx,nz)

nHtot: ndarray[tuple[int, int], float64]: The total hydrogen number density. UNIT: cm^-3, DIMS: (nx,nz)

muH: float: The conversion constant from nHtot to rhog. It is assumed that it is constant throughout the disk. It is given by rhog/nHtot. UNIT: g.

NHver: ndarray[tuple[int, int], float64]: Vertical total hydrogen column density. nHtot is integrated from the disk surface to the midplane at each radial grid point. The intermediate results are stored at each grid point. For example NHver[:,0] gives the total column density as a function of radius. UNIT: cm^-2, DIMS: (nx,nz).

NHrad: ndarray[tuple[int, int], float64]: Radial total hydrogen column density. Integrated along radial rays, starting from the star. Otherwise same behaviour as NHver. UNIT: cm^-2, DIMS: (nx,nz)

nd: ndarray[tuple[int, int], float64]: The dust number density. UNIT: cm^-3, DIMS: (nx,nz).

tg: ndarray[tuple[int, int], float64]: The gas temperature. UNIT: K, DIMS: (nx,nz).

td: ndarray[tuple[int, int], float64]: The dust temperature. UNIT: K, DIMS: (nx,nz).

pressure: ndarray[tuple[int, int], float64]: The gas pressure. UNIT: erg cm^-3, DIMS: (nx,nz).

soundspeed: ndarray[tuple[int, int], float64]: The isothermal sound speed. UNIT: km s^-1, DIMS: (nx,nz).

velocity: ndarray[tuple[int, int], float64]: The velocity field (vector) given as vx,vy,vz. UNIT: km s^-1, DIMS: (nx,nz,3).

damean: ndarray[tuple[int, int], float64]: The mean dust particle radius. Is defined as <a3>**(1/3). UNIT: micron, DIMS: (nx,nz).

da2mean: ndarray[tuple[int, int], float64]: The surface weighted mean dust particle radius. UNIT: micron, DIMS: (nx,nz)

dNlayers: ndarray[tuple[int, int], float64]: The number of ice layers on the dust grains. DIMS: (nx,nz).

Hx: ndarray[tuple[int, int], float64]: The X-ray energy deposition rate per hydrogen nuclei. UNIT: erg <H>^-1, DIMS: (nx,nz).

zetaX: ndarray[tuple[int, int], float64]: X-ray ionisation rate per hydrogen nuclei. UNIT: s^-1, DIMS: (nx,nz).

zetaCR: ndarray[tuple[int, int], float64]: Cosmic-ray ionisation rate per molecular hydrogen (H2) UNIT: s^-1, DIMS: (nx,nz).

zetaSTCR: ndarray[tuple[int, int], float64]: Stellar energetic particle ionisation rate per H2. UNIT: s^-1, DIMS: (nx,nz).

tauX1: ndarray[tuple[int, int], float64]: Radial optical depth at 1 keV (for X-rays). DIMS: (nx,nz).

tauX5: ndarray[tuple[int, int], float64]: Radial optical depth at 5 keV (for X-rays). DIMS: (nx,nz).

tauX10: ndarray[tuple[int, int], float64]: Radial optical depth at 10 keV (for X-rays). DIMS: (nx,nz).

AVrad: ndarray[tuple[int, int], float64]: Radial visual extinction (measured from the star outwards). DIMS: (nx,nz).

AVver: ndarray[tuple[int, int], float64]: Vertical visual extinction (measured from the disk surface to the mid-plane). UNIT:, DIMS: (nx,nz).

AV: ndarray[tuple[int, int], float64]: Given by min([AVver[ix,iz],AVrad[ix,iz],AVrad[nx-1,iz]-AVrad[ix,iz]]), hence the lowest visual extinction at a certain point. It is assumed radiation can escape either vertically upwards, radially inwards or radially outwards. DIMS: (nx,nz).

nlam: int: The number of wavelength bands used in the continuum radiative transfer.

lams: ndarray[tuple[int, ...], dtype[float64]]: The band wavelengths (centers) considered in the radiative transfer. UNIT: microns, DIMS: (nlam).

lambs: ndarray[tuple[int, ...], dtype[float64]] | None: The band wavelengths (edges) considered in the radiative transfer. UNIT: microns, DIMS: (nlam+1). This can be None as the edges wavelengths are only available with XrayRT or gasRT mode.

chi: ndarray[tuple[int, ...], dtype[float64]]: Geometrical UV radiation field in units of the Drain field. UNIT: Draine field, DIMS: (nx,nz).

chi1_shielded: ndarray[tuple[int, ...], dtype[float64]]: The shielded UV radiation field (i.e. self-shielding is taken into account). UNIT: Draine field 91.2-205 nm, DIMS: (nx,nz).

chi1_unshielded: ndarray[tuple[int, ...], dtype[float64]]: The unshielded UV radiation field (no self-shielding). UNIT: Draine field 91.2-205 nm, DIMS: (nx,nz).

chi2_shielded: ndarray[tuple[int, ...], dtype[float64]]: The shielded UV radiation field (i.e. self-shielding is taken into account). UNIT: Habing field 91.2-111nm, DIMS: (nx,nz).

chi2_unshielded: ndarray[tuple[int, ...], dtype[float64]]: The unshielded UV radiation field (no self-shielding). UNIT: Habing field 91.2-111 nm, DIMS: (nx,nz).

chiRT: ndarray[tuple[int, ...], dtype[float64]]: UV radiation field as properly calculated in the radiative transfer, in units of the Drain field. UNIT: Draine field, DIMS: (nx,nz).

kappaRos: ndarray[tuple[int, int], float64] | None: Rosseland mean opacity. In case of gas radiative transfer for the dust plus the gas. UNIT: cm^-1, DIMS: (nx,nz).

tauchem: ndarray[tuple[int, int], float64]: Chemical timescale (steady-state) UNIT: yr, DIMS: (nx,nz).

taucool: ndarray[tuple[int, int], float64]: Cooling timescale. UNIT: yr, DIMS: (nx,nz).

taudiff: ndarray[tuple[int, int], float64]: Vertical radiative diffusion timescale (using the Rosseland mean opacities). UNIT: yr, DIMS: (nx,nz).

spnames: Dict[str, int]: Dictionary providing the index of a particular species (e.g. spnames[“CO”]). This index can be used for arrays having an species dimension (like nmol). The electron is included. DIMS: (nspec).

solved_chem: ndarray[tuple[int, ...], dtype[int32]]: Flag for chemistry solver (values, 0,1,2) (see ProDiMo code). 1 means everything okay, 0 failure, 2 time-dependent step needed. DIMS: (nx,nz).

rateH2form: ndarray[tuple[int, int], float64]: The H2 formation ratio. UNIT: s^-1, DIMS: (nx,nz).

rateH2diss: ndarray[tuple[int, int, int], float64]: The different H2 dissociation rates. UNIT: s^-1, DIMS: (nx,nz,3).

isoratio_12CO13CO: ndarray[tuple[int, int], float64] | None: Isotope ratio of 12CO to 13CO, if it was estimated in the model. DIMS: (nx,nz)

heat_names: list[str]: All the names of the heating processes.

cool_names: list[str]: All the names of the cooling processes.

heat_mainidx: ndarray[tuple[int, ...], dtype[int32]]: Index of the main heating process at the given grid point. UNIT: , DIMS: (nx,nz).

cool_mainidx: ndarray[tuple[int, ...], dtype[int32]]: Index of the main cooling process at the given grid point. UNIT: , DIMS: (nx,nz).

lineEstimates: list[DataLineEstimate] | None: All the line estimates from FlineEstimates.out. Each spectral line in FlineEstimates corresponds to one prodimopy.read.DataLineEstimate object.

lines: list[DataLine] | None: All the spectral lines from line_flux.out (proper Line transfer). Each spectral line in line_flux.out is represented by one prodimopy.read.DataLine object. This also includes multiple inclinations.

contImages: DataContinuumImages | None: Holds the continuum images data (image.out) if available. The full images are only read if requested for a particular wavelength.

star: DataStar | None: The stellar properies including the (unattenuated) input spectrum.

gas: DataGas | None: Holds various properties of the gas component (e.g. opacities).

dust: DataDust | None: Holds various properties of the dust component (e.g. opacities).

env_dust: DataDust | None: Holds various properties of the dust component (e.g. opacities) of the envelope. Only relevant if ProDiMo is used in the envelope mode.

dust_charge_mean: ndarray[tuple[int, ...], dtype[float64]] | None: The mean dust charge, if it was calculated. DIMS: (nx,nz).

dust_charge_min: ndarray[tuple[int, ...], dtype[float64]] | None: The minimum dust charge, if it was calculated. DIMS: (nx,nz).

dust_charge_max: ndarray[tuple[int, ...], dtype[float64]] | None: The max dust charge, if it was calculated. DIMS: (nx,nz).

elements: DataElements: Holds the element abundances (Elements.in)

species: DataSpecies: Holds the initial species data (Species.in).

sedObs: DataContinuumObs | None: Holds the provided SED observations (photometry, spectra, extinction etc.).

lineObs: list[DataLineObs] | None: Holds the provide line observations (e.g. LINEObs.dat and line profiles).

FLiTsSpec: DataFLiTsSpec | None: Holds the FLiTs spectrum if it exists.

property mstar: float: The stellar mass in solar units. Is taken from ProDiMo.out Deprecated: use mass instead.

property starSpec: DataStar: The stellar properties including the (unattenuated) input spectrum. Deprecated: use star instead.

property sed: DataSED | None: Not sure if this is the best solution, but need a getter. It returns the SED for the current inclination.

property nmol: ndarray[tuple[int, ...], dtype[float64]]: Number densities of all chemical species (mainly molecules). UNIT: cm^-3, DIMS: (nx,nz,nspec).

property cdnmol: ndarray[tuple[int, ...], dtype[float64]]: Vertical column number densities for each chemical species at each point in the disk. Integrated from the surface to the midplane at each radial grid point. UNIT: cm^-2, DIMS: (nx,nz,nspec).

property rcdnmol: ndarray[tuple[int, ...], dtype[float64]]: Radial column number densities for each species at each point in the disk. Integrated from the star outwards along fixed radial rays given by the vertical grid. UNIT: cm^-2, DIMS: (nx,nz,nspec)

property sdg: ndarray[tuple[int, ...], dtype[float64]]: The gas vertical surface density. This is for only one half of the disk (i.e. from z=0 to z+), like in ProDiMo. For the total surface density multiply by two. UNIT: g cm^-2, DIMS: (nx,nz)

property sdd: ndarray[tuple[int, ...], dtype[float64]]: The dust vertical surface density. This is for only one half of the disk (i.e. from z=0 to z+), like in ProDiMo. For the total surface density multiply by two. UNIT: g cm^-2, DIMS: (nx,nz)

property mgasdisk: float

The total disk gas mass (integrated rhog).

Returns:: Disk gas mass in g.
Return type:: float

property mdustdisk: float

The total disk dust mass (integrated rhod).

Returns:: Disk dust mass in g.
Return type:: float

property cool: ndarray[tuple[int, ...], dtype[float64]]: Cooling rates for the various cooling processes. UNIT: erg cm^-3 s^-1 DIMS: (nx,nz,ncool).

property heat: ndarray[tuple[int, ...], dtype[float64]]: Heating rates for the various heating processes. UNIT: erg cm^-3 s^-1 DIMS: (nx,nz,nheat).

property radFields: ndarray[tuple[int, ...], dtype[float64]]: Radiation field (mean intensity) for each wavelength band. UNIT: erg s^-1 cm^-2 sr^-1 Hz^-1, DIMS: (nx,nz,nlam).

property vol: ndarray[tuple[int, ...], dtype[float64]]: The volume for each grid point. UNIT: cm³, DIMS: (nx,nz).

property chemnet: ReactionNetworkPout: The chemical reactions of the model. Holds what is in Reactions.out. Is only read if needed.

property velocity_xz: ndarray[tuple[int, ...], dtype[float64]]: : The absolute xz-component of the velocity field. UNIT: km s^-1, DIMS: (nx,nz).

sedinc(incidx=0)[source]

Get the SED for for the given inclination index iinc.

Parameters:: incidx (int) – The inclination index for which the SED should be retrieved.
Returns:: The SED for the specified inclination index.
Return type:: DataSED

getSEDAnaMask(lam)[source]

Returns a numpy mask where only the grid points outside the emission origin area for the given wavelength (continuum) are marked as invalid.

This mask can be used in e.g. avg_quantity()

Warning

in case of optically thin emission only the points of one half of the disk are considered. In that case an average quantity of this box(mask) will not be completely correct.

Parameters:: lam (float) – the wavelength for which the emission origin region should be determined: UNITS: micron
Returns:: Numpy mask with DIMS: (nx,nz)
Return type:: array_like(float,ndim=2)

getLine(wl, ident=None, incidx=0)[source]

Returns the spectral line (from the line RT) closest to the given wavelength.

Parameters:

wl (float) – The wavelength which is used for the search. UNIT: micron.
ident (str) –
A line identification string which should also be considered for the search.

Please note that it is not the chemical species name. For example if you have to lines for the same CO transition but one for NLTE (ident=CO) and one for LTE (ident=CO_lte) the species will be the same in both. To properly select those lines you need to provided the ident as the wl is identical.
incidx (int) – select the inclination for the line. (default: 0) . 0 means the first one. If only one inclination exists always this one will be used (i.e. the value of incidx is ignored).

Returns:

Returns None if no lines are included in the model.

Return type:

DataLine

selectLines(ident)[source]

Returns a list of all line fluxes included in the line transfer, for the given line ident.

Parameters:: ident (str) – The line identification (species name) as defined in ProDiMo. The ident is not necessarily equal to the underlying chemical species name (e.g. isotopologues, ortho-para, or cases like N2H+ and HN2+)
Returns:: All lines for the given ident, or an empty list if nothing was found. None if no lines are available at all.
Return type:: list[DataLine]

gen_specFromLineEstimates(wlrange=[10, 15], ident=None, specR=3000, unit='W', contOnly=False, noCont=False, faceon=False)[source]

Generates a “Spectrum” from the line estimates results of ProDiMo and convolves it to the given spectral resolution. If the SED was also calculated the line fluxes will be added to the continuum, if not a continuum of zero flux is assumed.

This routine can become very slow, depending on the number of lines within a given wavelength range. It can be more efficient to produce smaller chunks with not so many lines.

Parameters:

wlrange (array_like) – Generate the spectrum in the wavelength range [start,end] Default: [10,15] Units: micron
ident (str | list[str] | None) – only the lines with this ident (like given in the line estimates) are considered. Default: None (all lines in the given wlrange are considered) if ident is a list of strings, all those molecules will be considered
specR (int) – the desired spectral resolution. If None a spectrum with simple the line fluxes added (i.e. as “delta function”) is returned. Default: 3000
unit (str) – desired unit for the output. Current choices W (W/m^2/Hz) or ‘Jy’ (Jansky). W is the default option.
contOnly (boolean) – only do it for the continuum (do not add any lines). Default: False
noCont (boolean) – assume zero continuum Default: False
faceon (bool) – whether to use the face-on SED for the continuum calculation. Default: False

Returns:

tuple containing: wls(array_like): array of wavelength points in micron, flux(array_like): flux values for each wavelength point in W m^-2 Hz^-1 or Jy (depending on the unit parameter)

Return type:

(tuple)

getLineEstimate(ident, wl)[source]

Finds and returns the line estimate for the given ident which is closest to the given wavelength. If the ident is unknown None is returned. If the ident exist the routine always returns an line estimate object.

Parameters:

ident (string) – The line identification string. Note that the ident is not necessarily equal to the corresponding chemical species.
wl (float) – the wavelength which is used for the search. UNIT: micron.

Returns:

Returns None if the line estimate is not found.

Return type:

prodimopy.read.DataLineEstimate

Notes

The parameters of this method are not consistent with getLine(). This might be confusing. However, there is a reason for this. The number of included line estimates in a ProDiMo model can be huge and just searching with the wavelength might become slow. However, this probably can be improved.

To speed up the reading of FlineEstimate.out the extra radial info (see DataLineEstimateRInfo of all the line estimates is not read. However, in this routine it is read for the single line estimate found. One should use this routine to access the radial info for a single line estimate.

selectLineEstimates(ident, wlrange=None)[source]

Returns a list of all line estimates (i.e. all transitions) for the given line ident and/or in the given wavelength range.

Parameters:

ident (str | None) – The line identification (species name) as defined in ProDiMo. The ident is not necessarily equal to the underlying chemical species name (e.g. isotopologues, ortho-para, or cases like N2H+ and HN2+). If wlrange is set ident can be None in that case all lines in the given wavelength range are returned.
wlrange (All lines in the given wavelength range [start,end] and the given ident) – are returned. if the ident is None all lines are returned. Default: None Units: micron

Returns:

List of prodimopy.read.DataLineEstimate objects, or empty list if nothing was found.

Return type:

list(prodimopy.read.DataLineEstimate)

Notes

In this routine the additional radial information of the line estimate (see DataLineEstimateRInfo) is not read.

Examples

# select all line estimates for CO
lests=model.selectLineEstimates("CO")

# select only the 1300 micron CO line
lests=model.selectLineEstimates("CO",wlrange=[1300,1310])

# select all linest in the given wavelength range
lests=model.selectLineEstimates(None,wlrange=[1200,1400])

# print all of them to see what we got
for lest in lests:
  print(lest)

getLineOriginMask(lineEstimate)[source]

Returns a numpy mask where the grid points outside the emission origin area for the given lineEstimate are marked as invalid.

This mask can be used in e.g. avg_quantity()

Parameters:: lineEstimate (DataLineEstimate) – A line estimate object for which the mask should be created.
Returns:: Numpy mask with DIMS: (nx,nz)
Return type:: NDArray[np.float64]

getAbun(spname)[source]

Returns the abundance for a given species.

Parameters:: spname (string) – The name of the chemical species as define in ProDiMo.
Returns:: an array of dimension [nx,nz] with species abundance or None if the species was not found.
Return type:: array_like(float,ndim=2)

Notes

The abundance is given relative to the total hydrogen nuclei density (i.e. nmol(spname)/nHtot)

A warning is printed in case the species was not found in spnames.

get_TotSpeciesMass(spname=None)[source]

Returns the total mass (integrated over the whole disk) for all or one given species.

The total mass is calculated by summing up the product of the number density and the volume at each grid point and multiplying it with the mass of one molecule. A warning is printed in case the species was not found in spnames.

Parameters:: spname (string) – The name of the chemical species as defined in ProDiMo. If None an array with the species masses of all elements will be returned.
Returns:: The total mass of the species in grams or an array with the masses of all species. None if the species was not found.
Return type:: float or array_like(float)

get_toomreQ(mstar=None)[source]

Returns the Toomre Q parameter as a function of radius. (for the midplane).

Q is given by

Q(r)=k(r)*cs(r)/(pi * G * sigma(r))

for k we use the keplerian frequency, cs is the soundspeed (from the mdoel) and sigma is the total gas surface density (both half’s of the disk).

Parameters:: mstar (float) – The stellar mass in solar units (optional). If None the value from the model is taken.
Returns:: an array of dimension [nx] with the Q param
Return type:: array_like(float,ndim=1)

get_VKep(mstar=None, type=2)[source]

Returns the Keplerian rotation velocity [km/s] as a function of radius and height (same as in ProDiMo)

Type 1: vkep=sqrt((G* mstar) / r) Type 2: vkep=sqrt((G* mstar * x**2) / r**3) Type 3: with pressure gradient (see Rosenfeld+ 2013)

Parameters:: mstar (float) – The stellar mass in solar units (optional). If None the value from the model is taken.
Returns:: Rotation velocity [km/s]
Return type:: array_like(float,ndim=1)

get_KeplerOmega(mstar=None)[source]

Returns the Keplerian orbital frequency [1/s] (for the midplane).

omega=sqrt(G*mstar/r**3)

Parameters:: mstar (float) – The stellar mass in solar units (optional). If None the value from the model is taken.
Returns:: Keplerian orbital frequency as function of radius [1/s]
Return type:: array_like(float,ndim=1)

int_ring(quantity, r1=None, r2=None)[source]

Integrates the given quantity (needs to be a function of r) from rin to rout.

No interpolation is done. Simply the nearest neighbour for rin and rout are taken.

The integration is done in cgs units and also the return value is in cgs units.

For example if the total gas surface density is passed this routine gives the disk mass. (if r1=rin and r2=rout).

Parameters:

quantity (array_like(float,ndim=1)) – the quantity to average. DIMS: (nx).
r1 (float) – inner radius [au]. optional DEFAULT: inner disk radius
r2 (float) – inner radius [au]. optional DEFAULT: outer disk radius

avg_quantity(quantity, weight='gmass', weightArray=None, mask=None)[source]

Calculates the weighted mean value for the given field (e.g. td). Different weights can be used (see weight parameter)

Parameters:

quantity (array_like(float,ndim=2)) – the quantity to average. DIMS: (nx,nz).
weight (str) – option for the weight. Values: gmass (gas mass) dmass (dust mass) or vol (Volume). DEFAULT: gmass
weightArray (array_like(float,ndim=2)) – Same dimension as quantity. If weightArray is not None it is used as the weight for calculating the mean value.
mask (array_like(boolean,ndim=2)) – A mask with the same dimension as quantity. All elements with True are masked, i.e. not considered in the average (see numpy masked_array). For example one can pass the result of getLineOriginMask() to calculate average quantities over the emitting region of a line.

Examples

Here is an example for making a mass weighted average for the gas temperature in the main line emitting region for the CO 1.3 mm line. model is the model data read with read_prodimo().

avgtg=model.avg_quantity(model.tg,weight="gmass",
                  mask=model.getLineOriginMask(model.getLineEstimate("CO",1300.)))

analyse_chemistry(species, to_txt=True, filenameChemistry='chemanalysis.out', screenout=False)[source]

Routine to create a Chemistry analysis object for the given species, if the chemanalysis data is available.

Parameters:

species (str) – The species to analyse.
to_txt (bool) – Write info about formation and destruction reactions for the selected molecule to a txt file.
filenameChemistry (str) – The name of the file that holds the reaction rates. Just in case it has a different name, usually one does not need to change that.
screenout (bool) – If True all the form./dest. reactions are printed on the screen.

Returns:

Object that holds all the required information and can be use for the plotting routines or for reac_rates_ix_iz().

Return type:

Chemistry

freezeout_timescale(species, alpha=1.0)[source]

Calculates the freezeout timescale for the given species, on the full grid using Eq. 12 from Rab+ (2017) .

Parameters:

species (string) – The chemical species for which the freezeout timescale should be calculated.
alpha (float) – The sticking coefficient. Default: 1.0

Returns:

array_like(float,ndim=2) – The freezeout timescale for the given species in years. DIMS: (nx,nz)
.. todo:: –
- get sticking coefficient from parameter file. However, that might not be possible if an equation is used.

prodimopy.read.read_prodimo(directory='.', name=None, readlineEstimates=True, readObs=True, readImages=True, filename='ProDiMo.out', filenameLineEstimates='FlineEstimates.out', filenameLineFlux='line_flux.out', filenameFLiTs='specFLiTs.out', td_fileIdx=None)[source]

Reads in all (not all yet) the output of a ProDiMo model from the given model directory.

Parameters:

directory (str) – the directory of the model (optional). One can use ~ for the home directory
name (str) – an optional name for the model (e.g. can be used in the plotting routines)
readlineEstimates (boolean) – read the line estimates file (can be slow, so it is possible to deactivate it)
readObs (boolean) – try to read observational data (e.g. SEDobs.dat …)
filename (str) – the filename of the main prodimo output
filenameLineEstimates (str) – the filename of the line estimates output
filenameLineFlux (str) – the filename of the line flux output
filenameFLiTs (str) – the filename of the flits output
td_fileIdx (str | int) – in case of time-dependent model the index of a particular output age can be provided (e.g. “03”) if it is an int try to expand the index properly to a string (currently 1 becomes “0001”)

Returns:

the ProDiMo model data or None in case of errors.

Return type:

prodimopy.read.Data_ProDiMo

prodimopy.read.read_elements(directory, filename='Elements.out')[source]

Reads the Elements.out file. Also adds the electron “e-” to the Elements.

Parameters:

directory (str) – the directory of the model
filename (str) – an alternative Filename

Returns:

the Elements model data or None in case of errors.

Return type:

DataElements

prodimopy.read.read_species(directory, filename='Species.out')[source]

Reads the Species.out file. Also adds the electron “e-” if necessary.

Also sets data.nspec as only here I can know if e- was a real species in ProDiMo or not. If it was a real species I don’t need to add that.

Parameters:

directory (str) – the directory of the model
pdata (prodimopy.read.Data_ProDiMo) – the ProDiMo Model data structure, where the data is stored
filename (str) – an alternative Filename

Returns:

DataSpecies – the Species model data or None in case of errors.
..todo –
- stochiometry is still missing

Return type:

DataSpecies | None

prodimopy.read.read_lineEstimates(directory, pdata, filename='FlineEstimates.out')[source]

Read FlineEstimates.out Can only be done after ProDiMo.out is read.

Parameters:

directory (str) – the directory of the model
pdata (prodimopy.read.Data_ProDiMo) – the ProDiMo Model data structure, where the data is stored
filename (str) – an alternative Filename

prodimopy.read.read_linefluxes(directory, filename='line_flux.out')[source]

Reads the line fluxes output. Can deal with multiple inclinations.

Parameters:

directory (str) – the model directory.
filename (str) – the filename of the output file (optional).

Returns:

List of lines.

Return type:

list(prodimopy.read.DataLine)

prodimopy.read.read_lineObs(directory, nlines, filename='LINEobs.dat')[source]: Reads the lineobs Data. the number of lines have to be known.

prodimopy.read.read_lineObsProfile(filename, directory='.')[source]: reads a line profile file which can be used for ProDiMo

prodimopy.read.read_gas(directory, filename='gas_cs.out')[source]

Reads gas_cs.out

Return type:: DataGas

prodimopy.read.read_continuumImages(directory, filename='image.out')[source]

Reads the image.out file.

To avoid unnecessary memory usage, the Intensities are not read in this routine. For this use get_image()

prodimopy.read.read_continuumObs(directory, filename='SEDobs.dat')[source]

Read observational continuum data (SED).

Looks for and reads the following files:

SEDobs.dat photometric data points.
PREFIXspec.dat where prefix can be Spitzer, ISO, PACS, SPIRE, JWST
extinct.dat extinction data used to redden the simulated SED.

All files are optional. The information is then used to e.g. overplot the observational data in plot_sed()).

prodimopy.read.read_FLiTs(directory, filename='specFLiTs.out')[source]

Reads the FLiTs spectrum.

Returns:: The Spectrum.
Return type:: prodimopy.read.DataFLiTsSpec

prodimopy.read.read_sed(directory, filename='SED.out', filenameAna='SEDana.out')[source]: Reads the ProDiMo SED output including the analysis data.

prodimopy.read.read_star(directory, filenameSpec='StarSpectrum.out', filenameRTinterpol='RTinterpolation.out')[source]

Reads StarSpectrum.out and RTinterpolation.out if it exists.

Calls read_starSpec() for the StarSpectrum.out

Parameters:

directory (str) – The model directory.
filenameSpec (str) – The filename for the star spectrum. Default: “StarSpectrum.out”
filenameRTinterpol (str) – The filename for the RT interpolation data. Default: “RTinterpolation.out”

prodimopy.read.read_starSpec(directory, filename='StarSpectrum.out')[source]

Reads StarSpectrum.out

Parameters:

directory (str) – The model directory.
filename (str) – The filename. Default: “StarSpectrum.out”

prodimopy.read.read_bgSpec(directory, filename='BgSpectrum.out')[source]

Reads the BgSpectrum.out file.

Returns:

prodimopy.read.DataBgSpec
the background spectra or None if not found.

prodimopy.read.read_dust(directory, filename='dust_opac.out')[source]

Reads dust_opac.out and if it exists dust_sigmaa.out. Does not read the dust composition yet.

Return type:: DataDust

prodimopy.read.read_dust_sigmaa(directory, filename='dust_sigmaa.out')[source]

Reads the radial surface density profiles for each grain size.

Parameters:

directory (str) – The directory where the file is located (the ProDiMo model directory).
filename (str) – The filename. Default: “dust_sigmaa.out”

Return type:

tuple[ndarray[tuple[int, …], dtype[_ScalarType_co]], ndarray[tuple[int, …], dtype[_ScalarType_co]]]

prodimopy.read.read_dust_fsize(directory, pdata, filename='dust_fsize.out')[source]

Reads the dust size distribution function if available (either dust_fisize.out or dust_fsize_rho.out).

Parameters:

directory (str) – The directory where the file is located (the ProDiMo model directory).
pdata (Data_ProDiMo) – The model object.
filename (str) – The filename. Default: “dust_fsize.out”

Returns:

Dust size distribution funciton DIM: (na,nx,nz) na is the number of grain sizes.

Return type:

NDArray

prodimopy.read.read_dust_charge(directory, filename='dust_charge.out')[source]

Reads some dust charge information, if dust charge chemistry was used.

Parameters:

directory (str) – The directory where the file is located (the ProDiMo model directory).
filename (str) – The filename. Default: “dust_charge.out”

Returns:

mean dustcharc, minimum dust charge and maximum dustcharge

Return type:

tuple of NDArray

prodimopy.read.calc_NHrad_oi(data)[source]

Calculates the radial column density from out to in (border of model space to star/center)

Todo

is this really needed, should rather be a method within prodimopy.read.Data_ProDiMo

prodimopy.read.calc_surfd(data)[source]

Caluclates the gas and dust vertical surface densities at every point in the model.

Only one half of the disk is considered. If one needs the total surface density simply multiply the value from the midplane (zidx=0) by two.

prodimopy.read.analyse_chemistry(species, model, to_txt=True, td_fileIdx=None, filenameChemistry='chemanalysis.out', screenout=True)[source]

Deprecated since version 2.5: Use prodimopy.read.Data_ProDiMo.analyse_chemistry() instead.

Function that analyses the chemistry in a similar way to chemanalysis.pro.

Parameters:

species (str) – The species name one wants to analyze.
model (Data_ProDiMo) – the ProDiMo model data.
to_txt (boolean) – Write info about formation and destruction reactions for the selecte molecule to a txt file. Default: True
td_fileIdx (str or int) – For time-dependent chemanalysis. Provide here the idx for the timestep. E.g. “0001” for the first one.
filenameChemistry (str) – The name of the file that holds the reaction rates. Default: chemanalysis.out. Just in case it has a different name, usually one does not need to change that.
screenout (boolean) – If True all the form./dest. reactions are printed on the screen. Default: True

Returns:

Object that holds all the required information and can be use for the plotting routines, and to analyse the chemistry point by point func:~pread.Chemistry.reac_rates_ix_iz.

Return type:

prodimopy.read.Chemistry

prodimopy.read.reac_rates_ix_iz(model, chemana, ix=None, iz=None, locau=None, lowRatesFrac=0.001)[source]

Deprecated since version 2.5: Use prodimopy.read.Chemistry.reac_rates_ix_iz() instead.

Function that analyses the chemana manually via a point-by-point analysis for a given species. Shows the most important formation and destruction rates for the given point ix,iz (or in au via Parameter locau).

Parameters:

model (Data_ProDiMo) – the model data
chemana (Chemistry) – data resulting from analyse_chemistry() on a single species
ix (int) – ix corresponding to desired radial location in grid, starting at 0
iz (int) – iz corresponding to desired radial location in grid, starting at 0
locau (array_like) – the desired coordinates in au [x,z]. The routine then finds the closest grid point for those coordinates.
lowRatesFrac (float) – Only rates with rate/total_rate > lowRatesFrac are printed. Default: 1.e-3 This is useful to avoid printing all the reactions that are not important.