temet.cosmo package

Submodules

temet.cosmo.auxcatalog module

Cosmological simulations - auxiliary catalog for additional derived galaxy/halo properties. The functions here are rarely called directly. Instead they are typically invoked from within a particular auxCat request.

fofRadialSumType(sP, pSplit, ptProperty, rad, method='B', ptType='all')

Compute total/sum of a particle property (e.g. mass) for those particles enclosed within one of the SO radii already computed and available in the group catalog. Use one of four methods:

Method A: do individual halo loads per halo, one loop over all halos.
Method B: do a full snapshot load per type, then halo loop and slice per FoF, to cut down on I/O ops.
Method C: per type: full snapshot load, spherical aperture search per FoF (brute-force global).
Method D: per type: full snapshot load, construct octtree, spherical aperture search per FoF (global).

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
ptProperty (str) – particle/cell quantity to apply reduction operation to.
rad (str) – a radius available in the group catalog, e.g. ‘Group_R_Crit200’ [code units].
method (str) – see above. Note! Methods A and B restrict this calculation to FoF particles only, whereas method C does a full particle search over the entire box in order to compute the total/sum for each FoF halo.
ptType (str) – if ‘all’, then sum over all types (dm, gas, and stars), otherwise just for the single specified type.

Returns:

a 2-tuple composed of

result (ndarray): 1d array, result for each subhalo.
attrs (dict): metadata.

Warning

This was an early example of a catalog generating function, and is left mostly for reference as a particularly simple example. In practice, its functionality can be superseded by subhaloRadialReduction().

pSplitBounds(sP, pSplit, minStellarMass=None, minHaloMass=None, indivStarMags=False, partType=None, cenSatSelect=None, equalSubSplit=True)

For a given pSplit, determine an efficient work split and return the required processing for this task, in the form of the list of subhaloIDs to process and the global snapshot index range required in load to cover these subhalos.

Parameters:

pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
minStellarMass (float) – apply lower limit on mstar_30pkpc_log.
minHaloMass (float) – apply lower limit on mhalo_200_log.
indivStarMags (bool) – make sure return covers the full PartType4 size.
partType (str or None) – if not None, use this to decide particle-based split, otherwise use ‘gas’.
cenSatSelect (str or None) – if not None, restrict to ‘cen’, ‘sat’, or ‘all’.
equalSubSplit (bool) – subdivide a pSplit based on equal numbers of subhalos, rather than particles.

Returns:

a 3-tuple composed of

subhaloIDsTodo (list[int]): the list of subhalo IDs to process by this task.
indRange (dict): the index range for the particle load required to cover these subhalos.
nSubsSelection (int): the number of subhalos to be processed.

findHalfLightRadius(rad, vals, frac=0.5, mags=True)

Linearly interpolate in rr (squared radii) to find the half light radius, given a list of values[i] corresponding to each particle at rad[i].

Parameters:

rad (ndarray) – list of squared radii.
vals (ndarray) – list of values.
frac (float) – if 0.5, then half-light radius, otherwise e.g. 0.2 for r20.
mags (bool) – input vals are magnitudes, i.e. conversion to linear luminosity needed.

Returns:

half-light radius in 3D or 2D (if rad is input 3D or 2D).

Return type:

float

subhaloRadialReduction(sP, pSplit, ptType, ptProperty, op, rad, ptRestrictions=None, weighting=None, scope='subfind', minStellarMass=None, minHaloMass=None, cenSatSelect=None)

Compute a reduction operation (total/sum, weighted mean, etc) of a particle property (e.g. mass) for those particles of a given type, optionally enclosed within a given radius.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
ptType (str) – particle type e.g. ‘gas’, ‘stars’, ‘dm’, ‘bhs’.
ptProperty (str) – particle/cell quantity to apply reduction operation to.
op (str) – reduction operation to apply. ‘sum’, ‘mean’, ‘max’ or custom user-defined string.
rad (float or str) – if a scalar, then [physical kpc], otherwise a string label for a given radial restriction specification e.g. ‘rvir’ or ‘2rhalfstars’ (see _radialRestriction()).
ptRestrictions (dict) – apply cuts to which particles/cells are included. Each key,val pair in the dict specifies a particle/cell field string in key, and a [min,max] pair in value, where e.g. np.inf can be used as a maximum to enforce a minimum threshold only.
weighting (str) – if not None, then use this additional particle/cell property as the weight.
scope (str) – Calculation is restricted to subhalo particles only if scope=='subfind' (default), or FoF particles if scope=='fof'. If scope=='global', currently a full non-chunked snapshot load and brute-force distance computations to all particles for each subhalo (can change to tree method).
minStellarMass (str or float) – minimum stellar mass of subhalo to compute in log msun (optional).
minHaloMass (str or float) – minimum halo mass to compute, in log msun (optional).
cenSatSelect (str) – exclusively process ‘cen’, ‘sat’, or ‘all’.

Returns:

a 2-tuple composed of

result (ndarray): 1d array, value for each subhalo.
attrs (dict): metadata.

subhaloStellarPhot(sP, pSplit, iso=None, imf=None, dust=None, Nside=1, rad=None, modelH=True, bands=None, sizes=False, indivStarMags=False, fullSubhaloSpectra=False, redshifted=False, emlines=False, seeing=None, minStellarMass=None, minHaloMass=None)

Compute the total band-magnitudes (or half-light radii if sizes==True), per subhalo, under a given assumption of an iso(chrone) model, imf model, dust model, and radial restrction. If using a dust model, can include multiple projection directions per subhalo.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
iso (str) – isochrone library, as in cosmo.stellarPop.sps initialization.
imf (str) – stellar IMF, as in cosmo.stellarPop.sps initialization.
dust (str or None) – dust model, as in cosmo.stellarPop.sps initialization.
Nside (int or str or None) – if None, then no 2D projections are done, should be used if e.g. if no viewing angle dependent dust is requested. If an integer, then a Healpix specification, which defines the multiple viewing angles per subhalo. If str, one of ‘efr2d’ (several edge-on and face-on viewing angles), or ‘z-axis’ (single projection per subhalo along z).
rad (float or str) – if a scalar, then [physical kpc], otherwise a string label for a given radial restriction specification e.g. ‘rvir’ or ‘2rhalfstars’ (see _radialRestriction()).
modelH (bool) – use our model for neutral hydrogen masses (for extinction), instead of snapshot values.
bands (list[str] or None) – over-ride default list of broadbands to compute.
sizes (bool) – instead of band-magnitudes, save half-light radii.
indivStarMags (bool) – save the magnitudes for every PT4 (wind->NaN) in all subhalos.
fullSubhaloSpectra (bool) – save a full spectrum vs wavelength for every subhalo.
redshifted (bool) – all the stellar spectra/magnitudes are computed at sP.redshift and the band filters are then applied, resulting in apparent magnitudes. If False (default), stars are assumed to be at z=0, spectra are rest-frame and magnitudes are absolute (if dust is unresolved) or apparent (if dust is resolved).
emlines (bool) – include nebular emission lines.
seeing (float or None) – if not None, then instead of a binary inclusion/exclusion of each star particle based on the rad aperture, include all stars weighted by the fraction of their light which enters the rad aperture, assuming it is spread by atmospheric seeing into a Gaussian with a sigma of seeing [units of arcseconds at sP.redshift].
minStellarMass (str or float) – minimum stellar mass of subhalo to compute in log msun (optional).
minHaloMass (str or float) – minimum halo mass to compute, in log msun (optional).

Returns:

a 2-tuple composed of

result (ndarray): 1d array, value for each subhalo.
attrs (dict): metadata.

mergerTreeQuant(sP, pSplit, treeName, quant, smoothing=None)

For every subhalo, compute an assembly/related quantity using a merger tree.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
treeName (str) – specify merger tree.
quant (str) – specify subhalo quantity (available in merger tree) or custom user-defined quantity.
smoothing (None or list) – if not None, then smooth the quantity along the time dimension according to the tuple specification [method,windowSize,windowVal,order] where method should be mm (moving median of windowSize), ma (moving average of windowSize), or poly (poly fit of order N). The window size is given in units of windowVal which can be only snap.

Returns:

a 2-tuple composed of

result (ndarray): 1d array, value for each subhalo.
attrs (dict): metadata.

tracerTracksQuant(sP, pSplit, quant, op, time, norm=None)

For every subhalo, compute a assembly/accretion/related quantity using the pre-computed tracker track catalogs (through time).

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
quant (str) – specify tracer tracks catalog quantity.
op (str) – statistical operation to apply across all the child tracers of the subhalo.
time (str or None) – sample quantity at what time? e.g. ‘acc_time_1rvir’.
norm (str or None) – normalize quantity by a second quantity, e.g. ‘tvir_tacc’.

Returns:

a 2-tuple composed of

result (ndarray): 1d array, value for each subhalo.
attrs (dict): metadata.

subhaloCatNeighborQuant(sP, pSplit, quant, op, rad=None, proj2D=None, subRestrictions=None, subRestrictionsRel=None, minStellarMass=None, minHaloMass=None, cenSatSelect=None)

For every subhalo search for spatially nearby neighbors, using a tree, and compute some reduction operation over them. In this case, the search radius is globally constant and/or set per subhalo. Alternatively, perform an adaptive search such that we find at least N>=1 neighbor, and similarly compute a reduction operation over their properties.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
quant (str) – subhalo quantity to apply reduction operation to.
op (str) – reduction operation to apply. ‘min’, ‘max’, ‘mean’, ‘median’, and ‘sum’ compute over the requested quant for all nearby subhalos within the search, excluding this subhalo. ‘closest_rad’ returns the distance of the closest neighbor satisfying the requested restrictions. ‘d[3,5,10]_rad’ returns the distance to the 3rd, 5th, or 10th nearest neighbor, respectively. ‘closest_quant’ returns the quant of this closest neighbor. ‘count’ returns the number of identified neighbors.
rad (str, float, or list) – physical kpc if float, or a string as recognized by _radialRestriction(), or if a list/tuple, then we compute radial profiles for each subhalo instead of a single value, in which case the tuple provides the profile parameters {radMin,radMax,radNumBins,radBinsLog,radRvirUnits} as in subhaloRadialProfile().
proj2D (list or None) – if not None, do 3D profiles, otherwise 2-tuple specifying (i) integer coordinate axis in [0,1,2] to project along or ‘face-on’ or ‘edge-on’, and (ii) depth in code units (None for full box).
subRestrictions (list) – apply cuts to which subhalos are searched over. Each item in the list is a 3-tuple consisting of {field name, min value, max value}, where e.g. np.inf can be used as a maximum to enforce a minimum threshold only. This is the only option which modifies the search target sample, as opposite to the search origin sample.
subRestrictionsRel (dict) – as above, but every field is understood to be relative to the current subhalo value, which is the normalization, e.g. {‘gt’:1.0} requires that neighbors have a strictly larger value, while {‘lt’:0.5} requires neighbors have a value half as large or smaller.
minStellarMass (str or float) – minimum stellar mass of subhalo to compute in log msun (optional).
minHaloMass (str or float) – minimum halo mass to compute, in log msun (optional).
cenSatSelect (str) – exclusively process ‘cen’, ‘sat’, or ‘all’.

Returns:

a 2-tuple composed of

result (ndarray): 1d or 2d array, containing result(s) for each processed subhalo.
attrs (dict): metadata.

wholeBoxColDensGrid(sP, pSplit, species, gridSize=None, onlySFR=False, allSFR=False)

Compute a 2D grid of gas column densities [cm^-2] covering the entire simulation box. For example to derive the neutral hydrogen CDDF. Strategy is a chunked load of the snapshot files, for each using SPH-kernel deposition to distribute the mass of the requested species (e.g. HI, CIV) in all gas cells onto the grid.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
species (str) – the gas species/sub-component to grid.
gridSize (int or None) – if specified, override the default grid cell size [code units].
onlySFR (bool) – if True, only include SFR > 0 gas cells.
allSFR (bool) – if True, assume that all gas with SFR > 0 gas a unity fraction of the given species. E.g. for H2, assume star formation gas cells are entirely made of molecular hydrogen.

Returns:

a 2-tuple composed of

result (ndarray): 2d array, gridded column densities.
attrs (dict): metadata.

wholeBoxCDDF(sP, pSplit, species, gridSize=None, omega=False)

Compute the column density distribution function (CDDF, i.e. histogram) of column densities given a full box column density grid.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
species (str) – the gas species/sub-component to grid.
gridSize (int or None) – if specified, override the default grid cell size [code units].
omega (bool) – if True, then instead compute the single number Omega_species = rho_species / rho_crit,0.

Returns:

a 2-tuple composed of

result (ndarray): 1d array, the CDDF.
attrs (dict): metadata.

Note

There is unfortunate duplication/lack of generality between this function and wholeBoxColDensGrid() (e.g. in the projection depth specification) which is always called. To define a new catalog for a CDDF, it must be specified twice: the actual grid, and the CDDF.

subhaloRadialProfile(sP, pSplit, ptType, ptProperty, op, scope, weighting=None, proj2D=None, ptRestrictions=None, subhaloIDsTodo=None, radMin=-1.0, radMax=3.7, radBinsLog=True, radNumBins=100, radRvirUnits=False, Nside=None, Nngb=None, minHaloMass=None, minStellarMass=None, cenSatSelect='cen')

Compute subhalo radial profiles (e.g. total/sum, weighted mean, and so on) of a particle property (e.g. mass) for those particles of a given type.

Parameters:

sP (simParams) – simulation instance.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total].
ptType (str) – particle type, e.g. ‘gas’, ‘stars’, ‘dm’, ‘bhs’.
ptProperty (str) – particle/cell quantity to apply reduction operation to.
op (str) – reduction operation to apply. ‘min’, ‘max’, ‘mean’, ‘median’, and ‘sum’ for example.
scope (str) – which particles/cells are included for each profile. If scope=='global', then all snapshot particles are used, and we do the accumulation in a chunked snapshot load. Self/other halo terms are decided based on subfind membership, unless scope==’global_fof’, then on group membership. If scope=='fof' or ‘subfind’ then restrict to FoF/subhalo particles only, respectively, and do a restricted load according to pSplit. In this case, only the self-halo term is computed. If scope=='subfind_global' then only the other-halo term is computed, approximating the particle distribution using an already computed subhalo-based accumlation auxCat, e.g. ‘Subhalo_Mass_OVI’. If scope=='global_spatial', then use pSplit to decompose the work via a spatial subset of the box, such that we are guaranteed to have access to all the particles/cells within radMax of all halos being processed, to enable more complex global scope operations. If scope=='global_tngcluster', then use the file structure of the reconstructed TNG-Cluster simulation to achieve global scope profiles.
weighting (str) – if not None, then use this additional particle/cell property as the weight.
proj2D (list or None) – if not None, do 3D profiles, otherwise 2-tuple specifying (i) integer coordinate axis in [0,1,2] to project along or ‘face-on’ or ‘edge-on’, and (ii) depth in code units (None for full box).
ptRestrictions (dict) – apply cuts to which particles/cells are included. Each key,val pair in the dict specifies a particle/cell field string in key, and a [min,max] pair in value, where e.g. np.inf can be used as a maximum to enforce a minimum threshold only.
subhaloIDsTodo (list) – if not None, then process this explicit list of subhalos.
radMin (int) – minimum radius for profiles, should be [log] if radBinsLog == True, else [linear]. should be [code units] if radRvirUnits == False, else [dimensionless rvir units].
radMax (int) – maximum radius for profiles, as above.
radBinsLog (bool) – if True, input radMin and radMax are in log, otherwise linear.
radNumBins (int) – number of radial bins for profiles.
radRvirUnits (bool) – if True, radMin and radMax are in units of rvir of each halo. Note that in this case binning is always linear, whereas if False, binning is always logarithmic.
Nside (int or None) – if not None, should be a healpix parameter (2,4,8,etc). In this case, we do not compute a spherically averaged radial profile per halo, but instead a spherical healpix sampled set of shells/radial profiles, where the quantity sampling at each point uses a SPH-kernel with Nngb.
Nngb (int or None) – must be specified, if and only if Nside is also specified. The neighbor number.
minHaloMass (str or float) – minimum halo mass to compute, in log msun (optional).
minStellarMass (str or float) – minimum stellar mass of subhalo to compute in log msun (optional).
cenSatSelect (str) – exclusively process ‘cen’, ‘sat’, or ‘all’.

Returns:

a 2-tuple composed of

result (ndarray): 2d array, radial profile for each subhalo.
attrs (dict): metadata.

Note

For scopes global and global_fof, four profiles are saved: [all, self-halo, other-halo, diffuse], otherwise only a single ‘all’ profile is computed.

temet.cosmo.cloudy module

Use previously computed CLOUDY tables of photo-ionization models to derive ionic abundances and emissivities.

class cloudyIon(sP, el=None, res='lg_c17', redshiftInterp=False, order=3)

Bases: object

Use pre-computed Cloudy table to derive ionic abundances for simulation gas cells.

property atomicSymbols: List of atomic symbols, e.g. ‘Mg’, we have stored.

property atomicNames: List of atomic names, e.g. ‘Magnesium’, we have stored.

static atomicMasses(self): Dict of all atomic name:atomic weight pairs we have stored.

static ionList(): List of all elements + ion numbers we track.

atomicMass(elements): Return atomic mass (A_r in AMUs) for the given element(s).

numToRoman(ionNums): Map numeric ion numbers to roman numeral strings.

formatWithSpace(str, name=False): Convert a string of the type e.g. ‘Mg2’ or ‘O6’ to ‘Mg II’ or ‘O VI’.

slice(element, ionNum, redshift=None, dens=None, metal=None, temp=None)

Return a 1D slice of the table specified by a value in all other dimensions (only one input can remain None).

Parameters:

element (str) – name or symbol.
ionNum (str or int) – roman numeral (e.g. IV) or numeral starting at 1 (e.g. CII is ionNum=2).
redshift (float or None) – redshift.
dens (float or None) – hydrogen number density [log cm^-3].
metal (float or None) – metallicity [log solar].
temp (float or None) – temperature [log K].

Returns:

1d array of ionization fraction [log] as a function of the input which is None.

Return type:

ndarray

frac(element, ionNum, dens, metal, temp, redshift=None)

Interpolate the ion abundance table, return log(ionization fraction). Input gas properties can be scalar or np.array(), in which case they must have the same size. e.g. x = ion.frac('O','VI',-3.0,0.0,6.5) or x = ion.frac('O',6,dens=-3.0,metal=0.0,temp=6.5,redshift=2.2)

Parameters:

element (str) – name or symbol
ionNum (str or int) – roman numeral (e.g. IV) or numeral starting at 1 (e.g. CII is ionNum=2) where I = neutral (e.g. HeI = He), II = single ionized (e.g. HeII = He+) (e.g. HII region = fully ionized hydrogen, HeIII = fully ionized helium)
dens (ndarray) – hydrogen number density [log cm^-3]
temp (ndarrray) – temperature [log K]
metal (ndarray) – metallicity [log solar]

Returns:

ionization fraction per cell [log].

Return type:

ndarray

calcGasMetalAbundances(sP, element, ionNum, indRange=None, assumeSolarAbunds=False, assumeSolarMetallicity=False, sfGasTemp='cold')

Compute abundance mass fraction of the given metal ion for gas particles in the whole snapshot, optionally restricted to an indRange.

Parameters:

sP (simParams) – simulation instance.
element (str) – name or symbol.
ionNum (str or int) – roman numeral (e.g. IV) or numeral starting at 1 (e.g. CII is ionNum=2).
indRange (2-tuple) – if not None, the usual particle/cell-level index range to load.
assumeSolarAbunds (bool) – assume solar abundances (metal ratios), thereby ignoring GFM_Metals field.
assumeSolarMetallicity (bool) – assume solar metallicity, thereby ignoring GFM_Metallicity field.
sfGasTemp (str) – must be ‘cold’ (i.e. 1000 K), ‘hot’ (i.e. 5.73e7 K), ‘effective’ (snapshot value), or ‘both’, in which case abundances from both cold and hot phases are combined, each given fractional densities of the total gas density based on the two-phase model.

Returns:

mass fraction of the requested ion, relative to the total cell gas mass [linear].

Return type:

ndarray

class cloudyEmission(sP, line=None, res='lg_c23', redshiftInterp=False, order=3)

Bases: object

Use pre-computed Cloudy table to derive line emissivities for simulation gas cells.

lineWavelength(lines): Return the [rest, vacuum] wavelength of a line (or multiple lines) given its name, in [Ang].

slice(line, redshift=None, dens=None, metal=None, temp=None)

Return a 1D slice of the table specified by a value in all other dimensions (only one input can remain None).

Parameters:

line (str) – name of line (species, ion number, and wavelength triplet) (or abbreviation).
redshift (float or None) – redshift.
dens (float or None) – hydrogen number density [log cm^-3].
metal (float or None) – metallicity [log solar].
temp (float or None) – temperature [log K].

Returns:

1d array of volume emissivity [log erg/cm^3/s] as a function of the input which is None.

Return type:

ndarray

emis(line, dens, metal, temp, redshift=None)

Interpolate the line emissivity table for gas cell(s) with the given properties. Input gas properties can be scalar or np.array(), in which case they must have the same size.

Parameters:

line (str) – name of line (species, ion number, and wavelength triplet) (or abbreviation).
dens (float or None) – hydrogen number density [log cm^-3].
temp (float or None) – temperature [log K].
metal (float or None) – metallicity [log solar].
redshift (float or None) – if input, then interpolate also in redshift.

Returns:

1d array of volume emissivity, per cell [log erg/cm^3/s].

Return type:

ndarray

calcGasLineLuminosity(sP, line, indRange=None, dustDepletion=False, assumeSolarAbunds=False, assumeSolarMetallicity=False, sfGasTemp='cold')

Compute luminosity of line emission in for the given ‘line’, for gas particles in the whole snapshot, optionally restricted to an indRange.

Parameters:

sP (simParams) – simulation instance.
line (str) – name of line (species, ion number, and wavelength triplet) (or abbreviation).
indRange (2-tuple) – if not None, the usual particle/cell-level index range to load.
dustDepletion (bool) – apply a dust-depletion model for a given species.
assumeSolarAbunds (bool) – assume solar abundances (metal ratios), thereby ignoring GFM_Metals field.
assumeSolarMetallicity (bool) – assume solar metallicity, thereby ignoring GFM_Metallicity field.
sfGasTemp (str) – must be ‘cold’ (i.e. 1000 K), ‘hot’ (i.e. 5.73e7 K), ‘effective’ (snapshot value), or ‘both’, in which case abundances from both cold and hot phases are combined, each given fractional densities of the total gas density based on the two-phase model.

Returns:

luminosity of line emission, per cell [linear erg/s].

Return type:

ndarray

temet.cosmo.cloudyGrid module

Run grids of CLOUDY photo-ionization models for ion abundances, emissivities, or cooling rates.

getEmissionLines(): Return the list of emission lines (lineList above) that we save from CLOUDY runs.

loadFG11UVB(redshifts=None): Load the Faucher-Giguerre (2011) UVB at one or more redshifts and convert to CLOUDY units.

loadFG20UVB(redshifts=None): Load the Faucher-Giguere (2020) UVB at redshifts (or all available) and convert to CLOUDY units.

loadUVB(uvb='FG11', redshifts=None): Load the UVB at one or more redshifts.

loadUVBRates(uvb='FG11'): Load the photoionization [1/s] and photoheating [erg/s] rates for a given UVB.

cloudyUVBInput(gv): Generate the cloudy input string for a given UVB, by default with the self-shielding attenuation (at >= 13.6 eV) using the Rahmati+ (2013) fitting formula.

makeCloudyConfigFile(gridVals): Generate a CLOUDY input config file for a single run.

runCloudySim(gv, temp): Create a config file and execute a single CLOUDY run (e.g. within a thread).

runGrid(redshiftInd, nThreads=71, res='lg', uvb='FG11'): Run a sequence of CLOUDY models over a parameter grid at a redshift (one redshift per job).

collectOutputs(res='lg', uvb='FG11'): Combine all CLOUDY outputs for a grid into our master HDF5 table used for post-processing.

collectEmissivityOutputs(res='lg', uvb='FG11'): Combine all CLOUDY (line emissivity) outputs for a grid into a master HDF5 table.

collectCoolingOutputs(res='grackle', uvb='FG11'): Combine all CLOUDY (cooling function) outputs for a grid into a master HDF5 table.

temet.cosmo.clustering module

Calculations for TNG clustering paper.

twoPointAutoCorrelationPeriodicCube(sP, cenSatSelect='all', minRad=10.0, numRadBins=20, colorBin=None, cType=None, mstarBin=None, mType=None, jackKnifeNumSub=4)

Calculate the two-point auto-correlation function (of galaxies/subhalos) in a periodic cube geometry.

Parameters:

sP (simParams) – simulation instance.
cenSatSelect (str) – restrict to one of ‘cen’, ‘sat’, or ‘all’ (default).
minRad (float) – minimum radius [code units].
numRadBins (int) – number of radial bins.
colorBin (2-tuple or None) – if not None, a [min,max] tuple giving the minimum and maximum value of galaxy color to be included in the calculation. In this case we also require the corresponding cType to be specified.
cType (2-tuple or None) – if colorBin is not None, then should be [bands,simColorsModel] specifying the two bands to derive the color, and the model (iso,imf,dust).
mstarBin (2-tuple or None) – if not None, a [min,max]] tuple giving the minimum and maximum value of galaxy stellar mass to be included in the calculation. In this case we also require the corresponding mType to be specified.
mType (str or None) – if mstarBin is not None, then should be a string defining the stellar mass field (of plot.quantities.simSubhaloQuantity()) and so also units.
jackKnifeNumSub (int or None) – if not None, gives the number \(N\) of linear subdivisions for jackknife error estimation, requiring \(N^3\) additional tpcf calculations.

Returns:

a 4-tuple composed of

rad (ndarray): distance bin mid-points [code units].
xi (ndarray): the two-point auto correlation function \(\chi(r) = \rm{DD}/\rm{RR} - 1\) where xi[i] is computed between rad[i]:rad[i+1].
xi_err (ndarray): one sigma errors derived from the covariance (None if jackKnifeNumSub is None).
covar (ndarray): 2d covariance matrix (None if jackKnifeNumSub is None).

Note

Automatically caching. Result is saved to sP.derivPath/clustering/.

twoPointAutoCorrelationParticle(sP, partType, partField, pSplit=None): Calculate the [weighted] two-point auto-correlation function in a periodic cube geometry. Instead of per-subhalo/galaxy, this is a per-particle based calculation. Given the prohibitive expense, we do this with a Monte Carlo sampling of the first term of the particle-particle tpcf, saving intermediate results which are then accumulated as available. If pSplit==None and no partial files exist, full/single compute and return. If pSplit!=None and the requested partial file exists, concatenate all and return intermediate result. if pSplit!=None and the requested partial file does not exist, then run a partial computation and save.

isolationCriterion3D(sP, rad_pkpc, cenSatSelect='all', mstar30kpc_min=9.0): For every subhalo, record the maximum nearby subhalo mass (a few types) which falls within a 3d spherical aperture of rad_pkpc. Look for only cenSatSelect types in this search.

conformityRedFrac(sP, radRange=[0.0, 20.0], numRadBins=40, isolationRadPKpc=500.0, colorBin=None, cType=None, mstarBin=None, mType=None, jackKnifeNumSub=4, cenSatSelectSec='all', colorSplitSec=None): Calculate the conformity signal for -isolated- primaries subject to: {colorBin, cType, mstarBin, mType}. Specifications for these 4 are the same as for twoPointAutoCorrelationPeriodicCube(). A series of numRadBins are linearly spaced from radRange[0] to radRange[1] in physical Mpc. cenSatSelect is applied to the secondary sample, as are {colorSplitSec,cType} which define the red/blue split for secondaries.

lightcone3DtoSkyCoords(pos, vel, sP, velType): Transform pos [Nx3] and vel[Nx3] in code units, representing particle or galaxy positions in the periodic cube, into (ra,dec,redshift). The observer is assumed to be placed at the origin of the cube, (0,0,0), and the view direction is right now hardcoded.

temet.cosmo.color module

Calculations for optical stellar light of galaxies and galaxy colors.

loadColors(sP, quantName): Wrap the below function including some logic to take only a single lowercase string ‘quantName’ input. Loads either a color (difference of two band magntiudes), or just a band magnitude, for every subhalo.

loadSimGalColors(sP, simColorsModel, colorData=None, bands=None, projs=None, rad=''): Load band-magnitudes either from snapshot photometrics or from auxCat SPS modeling, and return band magnitudes (if bands is a single band), convert to a color (if bands contains two bands), or return complete loaded data (if bands is None). If loaded data is passed in with bands, do then magnitude computation/slicing without re-loading.

stellarPhotToSDSSColor(photVector, bands): Convert the GFM_StellarPhotometrics[] or SubhaloStellarPhotometrics[] vector into a specified color, by choosing the right elements and handling any necessary conversions.

calcSDSSColors(bands, redshiftRange=None, eCorrect=False, kCorrect=False, petro=False): Load the SDSS data files and compute a requested color, optionally restricting to a given galaxy redshift range, correcting for extinction, and/or doing a K-correction.

calcMstarColor2dKDE(bands, gal_Mstar, gal_color, Mstar_range, mag_range, sP=None, simColorsModel=None, kCorrected=True): Quick caching of (slow) 2D KDE calculation of (Mstar,color) plane for SDSS z<0.1 points if kCorrected==False, then tag filename (assume this is handled prior in calcSDSSColors) if sP is None, otherwise for simulation (Mstar,color) points if sP is specified.

compareOldNewMags(): Compare stellar_photometrics and my new sdss subhalo mags, and BuserUconverted vs sdss_u.

plotDifferentUPassbands(): Buser’s U filter from BC03 vs. Johnson UX filter from Bessel+ 98.

temet.cosmo.hydrogen module

Modeling of the UVB and hydrogen states following Rahmati. Credit to Simeon Bird for many ideas herein. Also full box or halo based analysis of hydrogen/metal content.

photoCrossSec(freq, ion='H I')

Find photoionisation cross-section (for a given ion) in cm^2 as a function of frequency. This is zero for energies less than nuthr (13.6 eV = 1 Ryd, for HI), and then falls off like E^-3. From Verner+ (1996), the Opacity Project, values are from Table 1 of astro-ph/9601009.

Parameters:

freq (ndarray[float]) – frequency in eV (must be numpy array).
ion (str) – specify coefficients.

Returns:

cross-section [cm^2].

Return type:

ndarray[float]

photoCrossSecGray(freq, J_nu, ion)

Compute a gray i.e. frequency/spectrum-averaged cross section. :param freq: frequency i.e. energy [Rydberg]. :type freq: ndarray[float] :param J_nu: uvb intensity [linear erg/s/cm^2/Hz/sr]. :type J_nu: ndarray[float] :param ion: specify coefficients. :type ion: str

Returns:: cross-section [cm^2].
Return type:: float

photoRate(freq, J_nu, ion)

Compute the photoionization rate of (H, HeI, HeII) given the UVB intensity. :param freq: frequency i.e. energy [Rydberg]. :type freq: ndarray[float] :param J_nu: uvb intensity [linear erg/s/cm^2/Hz/sr]. :type J_nu: ndarray[float] :param ion: specify coefficients. :type ion: str

Returns:: rate [1/s].
Return type:: float

uvbEnergyDensity(freq, J_nu, eV_min=6.0, eV_max=13.6)

Compute the energy density of the UVB between eV_min and eV_max. :param freq: frequency i.e. energy [Rydberg]. :type freq: ndarray[float] :param J_nu: uvb intensity [linear erg/s/cm^2/Hz/sr]. :type J_nu: ndarray[float] :param eV_min: minimum energy [eV]. :type eV_min: float :param eV_max: maximum energy [eV]. :type eV_max: float

Returns:: energy density U [erg/cm^3].
Return type:: float

uvbPhotoionAtten(log_hDens, log_temp, redshift)

Compute the reduction in the photoionisation rate at an energy of 13.6 eV at a given density [log cm^-3] and temperature [log K], using the Rahmati+ (2012) fitting formula. Note the Rahmati formula is based on the FG09 UVB; if you use a different UVB, the self-shielding critical density will change somewhat.

For z < 5 the UVB is probably known well enough that not much will change, but for z > 5 the UVB is highly uncertain; any conclusions about cold gas absorbers at these redshifts need to marginalise over the UVB amplitude here.

At energies above 13.6eV the HI cross-section reduces like \(\nu^{-3}\). Account for this by noting that self-shielding happens when tau=1, i.e \(\tau = n*\sigma*L = 1\). Thus a lower cross-section requires higher densities. Assume then that HI self-shielding is really a function of tau, and thus at a frequency \(\nu\), the self-shielding factor can be computed by working out the optical depth for the equivalent density at 13.6 eV. ie, for \(\Gamma(n, T)\), account for frequency dependence with:

\(\Gamma( n / (\sigma(13.6) / \sigma(\nu) ), T)\).

So that a lower x-section leads to a lower effective density. Note Rydberg ~ 1/wavelength, and 1 Rydberg is the energy of a photon at the Lyman limit, ie, with wavelength 911.8 Angstrom.

neutral_fraction(nH, sP, temp=10000.0, redshift=None): The neutral fraction from Rahmati+ (2012) Eqn. A8.

get_H2_frac(nH)

Get the molecular fraction for neutral gas from the ISM pressure: only meaningful when nH > 0.1. From Bird+ (2014) Eqn 4, e.g. the pressure-based model of Blitz & Rosolowsky (2006).

Parameters:: nHI (ndarray[float]) – neutral hydrogen number density [cm^-3].

neutralHydrogenFraction(gas, sP, atomicOnly=True, molecularModel=None): Get the total neutral hydrogen fraction, by default for the atomic component only. Note that given the SH03 model, none of the hot phase is going to be neutral hydrogen, so in fact we should remove the hot phase mass from the gas cell mass. But this is subdominant and should be less than 10%. If molecularModel is not None, then return instead the H2 fraction itself, using molecularModel as a string for the particular H2 formulation. Note that in all cases these are ratios relative to the total hydrogen mass of the gas cell.

hydrogenMass(gas, sP, total=False, totalNeutral=False, totalNeutralSnap=False, atomic=False, molecular=False, indRange=None): Calculate the (total, total neutral, atomic, or molecular) hydrogen mass per cell. Here we use the calculations of Rahmati+ (2012) for the neutral fractions as a function of density. Return still in code units, e.g. [10^10 Msun/h].

calculateCDDF(N_GridVals, binMin, binMax, binSize, sP, depthFrac=1.0)

Calculate the CDDF (column density distribution function) f(N) given an input array of HI or metal column densities values [cm^-2], from a grid of sightlines covering an entire box.

Parameters:

N_GridVals – column density values in [log cm^-2].
binMin (float) – column densities in [log cm^-2].
binMax (float) – column densities in [log cm^-2].
binSize (float) – in [log cm^-2].
depthFrac (float) – is the fraction of sP.boxSize over which the projection was done (for dX).

temet.cosmo.kCorr module

Compute K-corrections following Chilingarian+ (2010).

kCorrections(filter_name, redshifts, color_name, color_values)

Obtain k-correction for a given filter, redshift, and color.

Parameters:

filter_name – name of the filter to calculate K-correction for, e.g. ‘u’, ‘g’, ‘J2’.
redshift – redshift of a galaxy, should be between 0.0 and 0.5.
color_name – name of the color, e.g. ‘u - g’, ‘g - r’, ‘ui’.
color_value – value of galaxy color, as specified in color_name.

Returns:

K-correction in specified filter for given redshift and color.

test(): Benchmark.

temet.cosmo.lightcone module

Cosmological simulations - generating and working with lightcones. (w/ Andres Aramburo-Garcia)

get_cone(sP, group, config, snap_index): Load coordinates of a particle type, or halos/subhalos, transform into the lightcone geometry, subset, and optionally load additional fields for lightcone mmembers. Return transformed positions, velocities, additional fields, and indices back into the original periodic snapshots.

lightcone_coordinates(sP, group, pos, vel, config, snap_index): Compute the ra, dec, and redshift given the position and velocity within the lightcone.

generate_lightcone(index_todo=None): Generate a lightcone from a set of saved snapshots of a cubic periodic cosmological volume. Our technique is to apply the volume remapping to reshape the cubic box into an elongated cuboid domain with significant extent along the line-of-sight distance, taken to be the x-axis, while (ra,dec) are mapped from the (y,z) coordinates. If index_todo is not None, then process only this single snapshot.

finalize_lightcone(): Write total counts in analogy to normal snapshots to facilitate loading.

plot_z_distributions(): Debugging. Load lightcone files and plot redshift distributions.

temet.cosmo.mergertree module

Cosmological simulations - working with merger trees (SubLink, LHaloTree, C-Trees).

loadMPB(sP, id, fields=None, treeName='SubLink', fieldNamesOnly=False): Load fields of main-progenitor-branch (MPB) of subhalo id from the given tree.

loadMDB(sP, id, fields=None, treeName='SubLink', fieldNamesOnly=False): Load fields of main-descendant-branch (MDB) of subhalo id from the given tree.

loadMPBs(sP, ids, fields=None, treeName='SubLink')

Load multiple MPBs at once (e.g. all of them), optimized for speed, with a full tree load (high mem). Basically a rewrite of illustris_python/sublink.py under specific conditions (hopefully temporary).

Parameters:

sP (simParams) – simulation instance.
ids (list[int]) – list of subhalo IDs to load.
fields (list[str]) – list of field names to load, or None for all (not recommended).
treeName (str) – which merger tree to use? ‘SubLink’ or ‘SubLink_gal’.

Returns:

Dictionary of MPBs where keys are subhalo IDs, and the contents of each dict value is another dictionary of identical stucture to the return of loadMPB().

Return type:

dict

loadTreeFieldnames(sP, treeName='SubLink'): Load names of fields available in a mergertree.

insertMPBGhost(mpb, snapToInsert=None): Insert a ghost entry into a MPB dict by interpolating over neighboring snapshot information. Appropriate if e.g. a group catalog if corrupt but snapshot files are ok. Could also be used to selectively wipe out outlier points in the MPB.

mpbPositionComplete(sP, id, extraFields=None): Load a particular MPB of subhalo id, and return it along with a filled version of SubhaloPos which interpolates for any skipped intermediate snapshots as well as back beyond the end of the tree to the beginning of the simulation. The return indexed by snapshot number.

mpbSmoothedProperties(sP, id, fillSkippedEntries=True, extraFields=None): Load a particular subset of MPB properties of subhalo id, and smooth them in time. These are currently: position, mass (m200_crit), virial radius (r200_crit), virial temperature (derived), velocity (subhalo), and are inside [‘sm’] for smoothed versions. Also attach time with snap/redshift. Note: With the Sublink* trees, the group properties are always present for all subhalos and are identical for all subhalos in the same group.

quantMPB(sim, mpb, quant): Return a particular quantity from the MPB, where quant is a string. A simplified version of e.g. simSubhaloQuantity(). Can be generalized in the future. Returned units should be consistent with simSubhaloQuantity() of the same quantity name.

debugPlot(): Testing MPB loading and smoothing.

plot_tree(sP, subhaloID, saveFilename, treeName='SubLink', dpi=100, ctName='inferno', output_fmt='png'): Visualize a full merger tree of a given subhalo. saveFilename can either be a string, BytesIO buffer, or None. If a buffer, output_fmt should specify the required format (e.g. ‘pdf’, ‘png’, ‘jpg’). If None, then the plot is rendered herein and the final image array (uint8) is returned, for e.g. image manipulation procedures.

test_plot_tree()

test_plot_tree_mem(haloID=19)

temet.cosmo.perf module

Performance, scaling and timings analysis.

getCpuTxtLastTimestep(filePath): Parse cpu.txt for last timestep number and number of CPUs/tasks and total CPU hours.

loadCpuTxt(basePath, keys=None, hatbMin=0, skipWrite=False): Load and parse Arepo cpu.txt, save into hdf5 format. If hatbMin>0, then save only timesteps with active time bin above this value.

loadTimebinsTxt(basePath): Load and parse Arepo timebins.txt, save into hdf5 format.

loadMemoryTxt(basePath): Parse the memory.txt file for a run (prelimiany).

plotCpuTimes(): Plot code time usage fractions from cpu.txt. Note that this function is being automatically run and the resultant plot uploaded to http://www.illustris-project.org/w/images/c/ce/cpu_tng.pdf as of May 2016 and modifications should be made with caution.

plotTimebins(): Plot analysis of timebins throughout the course of a run.

plotTimebinsFrame(pStyle='white', conf=0, timesteps=None): Plot analysis of timebins at one timestep.

scalingPlots(): Construct strong (fixed problem size) and weak (Npart scales w/ Ncores) scaling plots based on the Hazel Hen (2016) and Hawk (2021) tests.

temet.cosmo.skirt module

Prepare simulation outputs for SKIRT9 radiative transfer runs, execute, and vis/analyze.

createParticleInputFiles(sP, subhaloID, params, fofScope=False, star_model='fsps', dust_model='const04'): Create the .txt files representing the stellar populations, and gas/dust, to be read in by SKIRT.

driver_single(): Prepare and execute a single SKIRT run of one subhalo. Note that all input and output files are created in the working directory, and left there after for visualization. By default, three projections are made simultaneously: edge-on, face-on, and one random orientation.

driver_sample(sP, subhaloIDs): Prepare and execute a large number of jobs.

convolve_cube_with_filter(data, wavelengths, filterName, raw=False): Return convolution of IFU (x,y,wave) cube with a given broadband. If raw is True, then do not do any unit conversions (e.g. for stats0 analysis).

vis(subhaloID=343503): Load and visualize SKIRT output. This routine simply looks for files in the current working directory corresponding to the input subhalo ID.

plot_sed(subhaloID=343503): Plot a spectral energy distribution .dat output file. As above, simply looks for existing output files in the current working directory.

temet.cosmo.spectrum module

Synthetic absorption spectra: generation.

lsf_matrix(instrument)

Create a (wavelength-dependent) kernel, for the line spread function (LSF) of the given instrument.

Parameters:

instrument (str) – string specifying the instrumental setup.

Returns:

a 3-tuple composed of

lsf_mode (int): integer flag specifying the type of LSF.
lsf (ndarray): 2d array, with dimensions corresponding to wavelength of the instrumental grid, and discrete/pixel kernel size, respectively. Each entry is normalized to unity.
lsf_dlambda (ndarray): 1d array, the LSF FWHM at the same wavelengths.

create_wavelength_grid(line=None, instrument=None): Create a wavelength grid (i.e. x-axis of a spectrum) to receieve absorption line depositions. Must specify one, but not both, of either ‘line’ or ‘instrument’. In the first case, a local spectrum is made around its rest-frame central wavelength. In the second case, a global spectrum is made corresponding to the instrumental properties.

varconvolve(arr, kernel): Convolution (non-fft) with variable kernel.

deposit_single_line(wave_edges_master, tau_master, f, gamma, wave0, N, b, z_eff, debug=False)

Add the absorption profile of a single transition, from a single cell, to a spectrum. Global method, where the original master grid is assumed to be very high resolution, such that no sub-sampling is necessary (re-sampling onto an instrument grid done later).

Parameters:

wave_edges_master (array[float]) – bin edges for master spectrum array [ang].
tau_master (array[float]) – master optical depth array.
N (float) – column density in [1/cm^2].
b (float) – doppler parameter in [km/s].
f (float) – oscillator strength of the transition
gamma (float) – sum of transition probabilities (Einstein A coefficients) [1/s]
wave0 (float) – central wavelength, rest-frame [ang].
z_eff (float) – effective redshift, i.e. including both cosmological and peculiar components.
debug (bool) – if True, return local grid info and do checks.

Returns:

None.

create_spectra_from_traced_rays(sP, line, instrument, rays_off, rays_len, rays_cell_dl, rays_cell_inds, cell_dens, cell_temp, cell_vellos, nThreads=54)

Given many completed rays traced through a volume, in the form of a composite list of intersected cell pathlengths and indices, extract the physical properties needed (dens, temp, vellos) and create the final absorption spectrum, depositing a Voigt absorption profile for each cell.

Parameters:

sP (simParams) – simulation instance.
line (str) – string specifying the line transition.
instrument (str) – string specifying the instrumental setup.
rays_off (array[int]) – first entry from tuple return of util.voronoiRay.rayTrace().
rays_len (array[int]) – second entry from tuple return of util.voronoiRay.rayTrace().
rays_cell_dl (array[float]) – third entry from tuple return of util.voronoiRay.rayTrace().
rays_cell_inds (array[int]) – fourth entry from tuple return of util.voronoiRay.rayTrace().
cell_dens (array[float]) – gas per-cell densities of a given species [linear ions/cm^3]
cell_temp (array[float]) – gas per-cell temperatures [linear K]
cell_vellos (array[float]) – gas per-cell line of sight velocities [linear km/s]
z_lengths (array[float]) – the comoving distance to each z_vals relative to sP.redshift [pMpc]
z_vals (array[float]) – a sampling of redshifts, starting at sP.redshift
nThreads (int) – parallelize calculation using this threads (serial computation if one)

generate_rays_voronoi_fullbox(sP, projAxis=2, nRaysPerDim=2000, raysType='voronoi_rndfullbox', subhaloIDs=None, pSplit=None, integrateQuant=None, search=False)

Generate a large grid of (fullbox) rays by ray-tracing through the Voronoi mesh.

Parameters:

sP (simParams) – simulation instance.
projAxis (int) – either 0, 1, or 2. only axis-aligned allowed for now.
nRaysPerDim (int) – number of rays per linear dimension (total is this value squared).
raysType (str) – either ‘voronoi_fullbox’ (equally spaced), ‘voronoi_rndfullbox’ (random), or ‘sample_localized’ (distributed around a given set of subhalos).
subhaloIDs (list) – if raysType is ‘sample_localized’ (only), then a list of subhalo IDs.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total]. Note that we follow a spatial subdivision, so the total job number should be an integer squared.
integrateQuant (str) – if None, save rays for future use. otherwise, directly perform and save the integral of the specified gas quantity along each ray.
search (bool) – if True, return existing data only, do not calculate new files.

integrate_along_saved_rays(sP, field, pSplit=None)

Integrate a physical (gas) property along the line of sight, based on already computed and saved rays. The result has units of [pc] * [field] where [field] is the original units of the physical field as loaded, unless field is a number density, in which case the result (column density) is in [cm^-2].

Parameters:

sP (simParams) – simulation instance.
field (str) – any available gas field.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total]. Note that we follow a spatial subdivision, so the total job number should be an integer squared.

concat_integrals(sP, field)

Combine split files for line-of-sight quantity integrals into a single file.

Parameters:

sP (simParams) – simulation instance.
field (str) – any available gas field.

generate_spectra_from_saved_rays(sP, ion='Si II', instrument='4MOST-HRS', pSplit=None, solar=False)

Generate a large number of spectra, based on already computed and saved rays.

Parameters:

sP (simParams) – simulation instance.
ion (str) – space separated species name and ionic number e.g. ‘Mg II’.
instrument (str) – specify wavelength range and resolution, must be known in instruments dict.
pSplit (list[int]) – standard parallelization 2-tuple of [cur_job_number, num_jobs_total]. Note that we follow a spatial subdivision, so the total job number should be an integer squared.
solar (bool) – if True, do not use simulation-tracked metal abundances, but instead use the (constant) solar value.

concat_spectra(sP, ion='Fe II', instrument='4MOST-HRS', solar=False)

Combine split files for spectra into a single file.

Parameters:

sP (simParams) – simulation instance.
ion (str) – space separated species name and ionic number e.g. ‘Mg II’.
instrument (str) – specify wavelength range and resolution, must be known in instruments dict.
solar (bool) – if True, do not use simulation-tracked metal abundances, but instead use the (constant) solar value.

test_conv(): Debug check behavior and benchmark variable convolution.

test_conv_master(): Debug check convolving on master grid vs inst grid.

temet.cosmo.spectrum_analysis module

Synthetic absorption spectra: analysis and derived quantities.

load_spectra_subset(sim, ion, instrument, solar, mode, num=None, EW_minmax=None, dv=0.0)

Load a subset of spectra from a given simulation and ion.

Parameters:

sim (simParams) – simulation instance.
ion (str) – space separated species name and ionic number e.g. ‘Mg II’.
instrument (str) – specify wavelength range and resolution, must be known in instruments dict.
solar (bool) – if True, do not use simulation-tracked metal abundances, but instead use the (constant) solar value.
num (int) – how many individual spectra to show.
mode (str) – either ‘all’, ‘random’, ‘evenly’, or ‘inds’.
inds (list[int]) – if mode is ‘inds’, then the list of specific spectra indices to plot. num is ignored.
EW_minmax (list[float]) – minimum and maximum EW to plot [Ang].
dv (float) – if not zero, then take as a velocity window (+/-), convert the wavelength axis to velocity, and subset spectra to only this vel range.

load_absorber_spectra(sim, line, instrument, solar, EW_minmax=None, dwave=0.0, dv=0.0)

Load the (local) spectra for each absorber, from a given simulation and ion. Note that the absorber catalog is available for each line separately, so one must be specified, and this line is used for the EW_minmax restriction (and returned EW values). However, the flux is the observable i.e. the full spectrum including all lines of this ion.

Parameters:

sim (simParams) – simulation instance.
line (str) – transition name e.g. ‘Mg II 2796’, since each line has its own absorber catalog.
instrument (str) – specify wavelength range and resolution, must be known in instruments dict.
solar (bool) – if True, do not use simulation-tracked metal abundances, but instead use the (constant) solar value.
EW_minmax (list[float]) – minimum and maximum EW to plot [Ang].
dwave (float) – if not zero, then take as a wavelength window (+/-) around each absorber.
dv (float) – if not zero, then take as a velocity window (+/-), convert the wavelength axis to velocity, and subset spectra to only this vel range.

absorber_catalog(sP, ion, instrument, solar=False)

Detect and chatacterize absorbers, handling the possibility of one or possibly multiple per sightline, defined as separated but contiguous regions of absorption. Create an absorber catalog, i.e. counts and offsets per sightline, and compute their EWs.

Parameters:

sP (simParams) – simulation instance.
ion (str) – space separated species name and ionic number e.g. ‘Mg II’.
instrument (str) – specify wavelength range and resolution, must be known in instruments dict.
solar (bool) – if True, do not use simulation-tracked metal abundances, but instead use the (constant) solar value.

cell_to_absorber_map(sP, ion, instrument, solar=False)

For each absorber, defined as a contiguous set of pixels in a spectrum, identify the gas cells that contribute to this absorber. Return ‘spec_index’ is an ordered list of global gas cell indices, stored with a length (counts_abs) and offset (offset_abs) approach, where the length and offset are per absorber. Also derived and returned: the column density per absorber, and the fraction of this (total) column density that is recovered by the cells included here as ‘significantly’ contributing to the absorber.

Parameters:

sP (simParams) – simulation instance.
ion (str) – space separated species name and ionic number e.g. ‘Mg II’.
instrument (str) – specify wavelength range and resolution, must be known in instruments dict.
solar (bool) – if True, do not use simulation-tracked metal abundances, but instead use the (constant) solar value.

calc_statistics_from_saved_rays(sP, ion)

Calculate useful statistics based on already computed and saved rays. Results depend on ion, independent of actual transition. Results depend on the physical properties along each sightline, not on the absorption spectra.

Parameters:

sP (simParams) – simulation instance.
ion (str) – space separated species name and ionic number e.g. ‘Mg II’.

test_abs_coldens(sim, ion, instrument): Check if the column densities of all absorbers of a given spectrum sum to the separately integrated total column density of that spectrum.

temet.cosmo.spectrum_test module

Synthetic absorption spectra generation (testing/non-production).

create_spectrum_from_traced_ray(sP, line, instrument, cell_dens, cell_dx, cell_temp, cell_vellos): Given a completed (single) ray traced through a volume, and the properties of all the intersected cells (dens, dx, temp, vellos), create the final absorption spectrum, depositing a Voigt absorption profile for each cell.

deposit_single_line_local(wave_edges_master, tau_master, f, gamma, wave0, N, b, z_eff, debug=False)

Add the absorption profile of a single transition, from a single cell, to a spectrum. Local method, where master grid is sub-sampled and optical depth is deposited back onto master.

Parameters:

wave_edges_master (array[float]) – bin edges for master spectrum array [ang].
tau_master (array[float]) – master optical depth array.
N (float) – column density in [1/cm^2].
b (float) – doppler parameter in [km/s].
f (float) – oscillator strength of the transition
gamma (float) – sum of transition probabilities (Einstein A coefficients) [1/s]
wave0 (float) – central wavelength, rest-frame [ang].
z_eff (float) – effective redshift, i.e. including both cosmological and peculiar components.
debug (bool) – if True, return local grid info and do checks.

Returns:

None.

generate_spectrum_uniform_grid(): Generate an absorption spectrum by ray-tracing through a uniform grid (deposit using sphMap).

generate_spectrum_voronoi(use_precomputed_mesh=True, compare=False, debug=1, verify=True)

Generate a single absorption spectrum by ray-tracing through the Voronoi mesh.

Parameters:

use_precomputed_mesh (bool) – if True, use pre-computed Voronoi mesh connectivity from VPPP, otherwise use tree-based, connectivity-free method.
compare (bool) – if True, run both methods and compare results.
debug (int) – verbosity level for diagnostic outputs: 0 (silent), 1, 2, or 3 (most verbose).
verify (bool) – if True, brute-force distance calculation verify parent cell at each step.

generate_spectra_voronoi_halo(): Generate a large grid of (halocentric) absorption spectra by ray-tracing through the Voronoi mesh.

benchmark_line(): Deposit many random lines.

temet.cosmo.stellarPop module

Stellar population synthesis, evolution, photometrics.

class sps(sP, iso='padova07', imf='chabrier', dustModel='cf00_res_conv', order=3, redshifted=False, emlines=False)

Bases: object

Use pre-computed FSPS stellar photometrics tables to derive magnitudes for simulation stars.

basePath = '/vera/u/dnelson/temet/temet//tables/fsps/'

imfTypes = {'chabrier': 1, 'kroupa': 2, 'salpeter': 0}

isoTracks = ['mist', 'padova07', 'parsec', 'basti', 'geneva']

stellarLib = ['miles', 'basel', 'csk']

dustModels = ['none', 'cf00', 'cf00_res_eff', 'cf00b_res_conv', 'cf00_res_conv', 'cf00_res3_conv']

bands = ['wfpc2_f439w', 'isaac_ks', 'iras_12', 'jwst_f356w', 'scuba_450wb', 'wfc3_uvis_f218w', 'pacs_100', 'buser_b', 'i2300', 'wfc3_ir_f140w', 'wfc_acs_f814w', 'wfc3_uvis_f814w', 'nicmos_f160w', 'wfcam_k', 'steidel_un', 'cousins_r', 'iras_25', 'wfcam_z', 'wfcam_y', 'wfc_acs_f625w', 'wfc3_ir_f098m', 'cousins_i', 'newfirm_j1', 'newfirm_j3', 'i1500', 'newfirm_j2', 'wise_w1', 'ps1_i', 'wise_w3', 'wise_w2', 'wise_w4', 'bessell_l', 'bessell_m', 'ps1_r', 'ps1_y', 'ps1_z', 'jwst_f277w', 'sdss_u', 'des_r', 'des_y', 'des_z', 'bessell_lp', 'wfc3_uvis_f438w', '2mass_ks', 'steidel_rs', 'des_g', 'des_i', 'uvot_m2', 'wfc3_ir_f110w', 'stromgren_u', 'stromgren_v', 'wfc_acs_f850lp', 'wfc3_uvis_f390w', 'stromgren_y', 'stromgren_b', 'wfpc2_f850lp', 'wfc3_uvis_f336w', 'spire_500', 'wfpc2_f555w', 'fors_v', 'uvot_w2', 'fors_r', 'wfc3_uvis_f606w', 'iras_100', 'newfirm_h2', 'newfirm_h1', 'wfc_acs_f475w', 'wfpc2_f814w', 'jwst_f070w', 'mips_24', 'wfc3_uvis_f775w', 'suprimecam_r', 'suprimecam_v', 'suprimecam_i', 'wfc_acs_f435w', 'suprimecam_b', 'mips_70', '2mass_h', '2mass_j', 'sdss_g', 'wfpc2_f300w', 'sdss_i', 'wfc3_uvis_f475w', 'sdss_r', 'wfc3_uvis_f850lp', 'newfirm_k', 'sdss_z', 'wfc3_ir_f160w', 'wfc_acs_f555w', 'wfc3_ir_f125w', 'i2800', 'vista_y', 'megacam_u', 'megacam_r', 'suprimecam_z', 'wfpc2_f450w', 'wfc3_uvis_f225w', 'wfc_acs_f606w', 'galex_nuv', 'vista_h', 'vista_k', 'vista_j', 'wfpc2_f336w', 'megacam_i', 'b', 'wfc3_ir_f105w', 'galex_fuv', 'mips_160', 'v', 'uvot_w1', 'jwst_f444w', 'wfc3_uvis_f555w', 'wfc3_uvis_f275w', 'wfcam_j', 'spire_350', 'wfcam_h', 'jwst_f200w', 'pacs_160', 'jwst_f150w', 'wfpc2_f255w', 'wfc_acs_f775w', 'jwst_f115w', 'spire_250', 'scuba_850wb', 'wfpc2_f606w', 'ps1_g', 'iras_60', 'nicmos_f110w', 'pacs_70', 'megacam_z', 'irac_4', 'irac_2', 'irac_3', 'irac_1', 'megacam_g', 'cfht_i', 'suprimecam_g', 'cfht_b', 'u', 'steidel_i', 'steidel_g', 'jwst_f090w', 'cfht_r']

computePhotTable(iso, imf, saveFilename): Compute a new photometrics table for the given (iso,imf,self.dust) using fsps.

filters(select=None): Return name of available filters.

has_filter(filterName): Return True or False if the pre-computed table contains the specified filter/band.

prep_dust_models(): Do possibly expensive pre-calculations for (resolved) dust model.

prep_filters(): Extract filter properties in case we want them later.

convertSpecToSDSSUnitsAndAttenuate(spec, output_wave=None): Convert a spectrum from FSPS in [Lsun/Hz] to SDSS units [10^-17 erg/cm^2/s/Ang] and, if self.sP.redshift > 0, attenuate the spectrum by the luminosity distance. If self.sP.redshift == 0, the rest-frame spectrum is returned at an assumed distance of 10pc (i.e. absolute magnitudes). If output_wave is not None, should be in Angstroms.

redshiftSpectrum(spec, output_wave=None): If self.sP.redshift > 0, attenuate the spectrum by the luminosity distance and redshift the wavelength. If self.sP.redshift == 0, the rest-frame spectrum is returned at an assumed distance of 10pc (i.e. absolute magnitudes). Note that the spectrum is sampled back onto the original output_wave points, which therefore become observer-frame instead of rest-frame.

mags(band, ages_logGyr, metals_log, masses_logMsun): Interpolate table to compute magnitudes [AB absolute] in requested band for input stars.

mags_code_units(sP, band, gfm_sftime, gfm_metallicity, masses_code, retFullSize=False): Do unit conversions (and wind particle filtering) on inputs, and return mags() results. If retFullSize is True, return same size as inputs with wind set to nan, otherwise filter out wind/nan values and compress return size.

calcStellarLuminosities(sP, band, indRange=None, rotMatrix=None, rotCenter=None): Compute (linear) luminosities in the given band, using either snapshot-stored values or on the fly sps calculation, optionally restricted to indRange. Note that wind is here returned as NaN luminosity, assuming it is filtered out elsewhere, e.g. in gridBox(). Return are linear luminosities in units of [Lsun/Hz].

dust_tau_model_mags(bands, N_H, Z_g, ages_logGyr, metals_log, masses_msun, ret_indiv=False, ret_full_spectrum=False, output_wave=None, rel_vel=None): For a set of stars characterized by their (age,Z,M) values as well as (N_H,Z_g) calculated from the resolved gas distribution, do the Model (C) attenuation on the full spectra, sum together, and convolve the resulting total L(lambda) with the transmission function of multiple bands, returning a dict of magnitudes, one for each band. If ret_indiv==True, then the individual magnitudes for every member star are instead returned separately. If ret_full_spectrum==True, the full aggregate spectrum, summed over all member stars, as a function of wavelength, is instead returned. If output_wave is not None, then this output spectrum is interpolated to the requested wavelength grid (in Angstroms, should be rest or observed frame depending on self.redshifted). If rel_vel is not None, then add LoS peculiar velocity shifts (physical km/s) of each star. Note: Will return apparent magnitudes if self.sP.redshift > 0, or absolute magnitudes if self.sP.redshift == 0. Will return redshifted (and attenuated) spectra if self.redshifted and self.sP.redshift is nonzero, otherwise rest-frame spectra.

resolved_dust_mapping(pos_in, hsml, mass_nh, quant_z, pos_stars_in, projCen, projVec=None, rotMatrix=None, pxSize=1.0): Compute line of sight quantities per star for a resolved dust attenuation calculation. Gas (pos,hsml,mass_nh,quant_z) and stars (pos_stars) are used for the gridding of the gas and the target (star) list. projVec is a [3]-vector, and the particles are rotated about projCen [3] such that it aligns with the projection direction. pxSize is in physical kpc.

debug_check_redshifting(redshift=0.8): Verify we understand what is going on with redshifting and apparent vs. absolute magnitudes of the band magnitudes (from FSPS) and the band magnitudes derived from our convolving our spectra with the bandpass filters manually.

debug_dust_plots(): Plot intermediate aspects of the resolved dust calculation.

debug_plot_spectra(): Check mock spectra.

debug_check_rawspec(): Check spectral tables.

temet.cosmo.util module

Helper functions related to cosmo box simulations.

redshiftToSnapNum(redshifts=None, times=None, sP=None, recalculate=False): Convert one or more input redshifts to closest matching snapshot numbers for a given sP.

validSnapList(sP, maxNum=None, minRedshift=None, maxRedshift=None, onlyFull=False, reqTr=False, reqFluidQuants=False): Return a list of all snapshot numbers which exist.

multiRunMatchedSnapList(runList, method='expand', **kwargs): For an input runList of dictionaries containing a sP key corresponding to a simParams for each run, produce a ‘matched’/unified set of snapshot numbers, one set per run, with all the same length, e.g. for comparative analysis at matched redshifts, or for rendering movie frames comparing runs at the same redshift. If method is ‘expand’, inflate the snapshot lists of all runs to the size of the maximal (duplicates are then guaranteed). If method is ‘condense’, shrink the snapshot lists of all runs to the size of the minimal (skips are then guaranteed).

snapNumToRedshift(sP, snap=None, time=False, all=False): Convert snapshot number(s) to redshift or time (scale factor or non-cosmological sim time).

snapNumToAgeFlat(sP, snap=None): Convert snapshot number to approximate age of the universe at that time.

crossMatchSubhalosBetweenRuns(sP_from, sP_to, subhaloInds_from_search, method='LHaloTree'): Given a set of subhaloInds_from_search in sP_from, find matched subhalos in sP_to. Can implement many methods. For now, uses external (pre-generated) postprocessing/SubhaloMatching/ for TNG_method runs, or postprocessing/SubhaloMatchingToDark/ for Illustris/TNG to DMO runs, or postprocessing/SubhaloMatchingToIllustris/ for TNG->Illustris runs. Return is an int32 array of the same size as input, where -1 indicates no match.

correctPeriodicDistVecs(vecs, sP): Enforce periodic B.C. for distance vectors (effectively component by component).

correctPeriodicPosVecs(vecs, sP): Enforce periodic B.C. for positions (add boxSize to any negative points, subtract boxSize from any points outside box).

correctPeriodicPosBoxWrap(vecs, sP): For an array of positions [N,3], determine if they span a periodic boundary (e.g. half are near x=0 and half are near x=BoxSize). If so, wrap the high coordinate value points by a BoxSize, making them negative. Suitable for plotting particle positions in global coordinates. Return indices of shifted coordinates so they can be shifted back, in the form of dict with an entry for each shifted dimension and key equal to the dimensional index.

periodicDists(pt, vecs, sP, chebyshev=False, xyzMin=False)

Calculate distances correctly taking into account periodic boundary conditions.

Parameters:

pt (list[3] or [N,3]) – if pt is one point, distance from pt to all vecs. if pt is several points, distance from each pt to each vec (must have same number of points as vecs).
vecs (list[3,N]) – position array in periodic 3D space.
sP (simParams) – simulation instance.
chebyshev (bool) – use Chebyshev distance metric (greatest difference in positions along any one axis)

periodicDists2D(pt, vecs, sP, chebyshev=False, xyzMin=False)

Calculate 2D distances correctly taking into account periodic boundary conditions.

Parameters:

pt (list[2] or [N,2]) – if pt is one point, distance from pt to all vecs. if pt is several points, distance from each pt to each vec (must have same number of points as vecs).
vecs (list[2,N]) – position array in periodic 3D space.
sP (simParams) – simulation instance.
chebyshev (bool) – use Chebyshev distance metric (greatest difference in positions along any one axis)

periodicDistsSq(pt, vecs, sP): As cosmo.util.periodicDists() but specialized, without error checking, and no sqrt. Works either for 2D or 3D, where the dimensions of pt and vecs should correspond.

periodicPairwiseDists(pts, sP): Calculate pairwise distances between all 3D points, correctly taking into account periodic B.C.

inverseMapPartIndicesToSubhaloIDs(sP, indsType, ptName, debug=False, flagFuzz=True, SubhaloLenType=None, SnapOffsetsSubhalo=None): For a particle type ptName and snapshot indices for that type indsType, compute the subhalo ID to which each particle index belongs. Optional: SubhaloLenType (from groupcat) and SnapOffsetsSubhalo (from groupCatOffsetListIntoSnap()), otherwise loaded on demand. If flagFuzz is True (default), particles in FoF fuzz are marked as outside any subhalo, otherwise they are attributed to the closest (prior) subhalo.

inverseMapPartIndicesToHaloIDs(sP, indsType, ptName, GroupLenType=None, SnapOffsetsGroup=None, debug=False): For a particle type ptName and snapshot indices for that type indsType, compute the halo/fof ID to which each particle index belongs. Optional: GroupLenType (from groupcat) and SnapOffsetsGroup (from groupCatOffsetListIntoSnap()), otherwise loaded on demand.

subhaloIDListToBoundingPartIndices(sP, subhaloIDs, groups=False, strictSubhalos=False): For a list of subhalo IDs, return a dictionary with an entry for each partType, whose value is a 2-tuple of the particle index range bounding the members of the parent groups of this list of subhalo IDs. Indices are inclusive in the sense of sP.snapshotSubset(). If groups==True, then the input IDs are assumed to be actually group (FoF) IDs, and we return bounding ranges for them. If strictSubhalos==True, then do not use parent groups to bound subhalo members, instead return exact bounding index ranges.

cenSatSubhaloIndices(sP=None, gc=None, cenSatSelect=None): Return a tuple of three sets of indices into the group catalog for subhalos: centrals only, centrals & satellites together, and satellites only.

subboxSubhaloCat(sP, sbNum): Generate/return the SubboxSubhaloList catalog giving the intersection of fullbox subhalos with the subbox volumes across time, using the merger trees, and some interpolated properties of relevant subhalos at each subbox snapshot.

subboxSubhaloCatExtend(sP, sbNum, redo=False): Extend the SubboxSubhaloList catalog with particle-level interpolated properties (i.e. M* < 30pkpc) of the relevant subhalos at each subbox snapshot. Separated into second step since this is a heavy calculation, restartable.

temet.cosmo.xray module

Generate x-ray emissivity tables using AtomDB/XSPEC and apply these to gas cells.

integrate_to_common_grid(bins_in, cont_in, bins_out): Convert a ‘compressed’ APEC spectrum into a normally binned one, by interpolating from an input (bins,cont) pair to (bins_out).

apec_convert_tables(): Load APEC tables (currently v3.0.9) and convert to a more suitable format for later use.

class xrayEmission(sP, instrument=None, order=3, use_apec=False)

Bases: object

Use pre-computed APEC or XSPEC table to derive x-ray emissivities for simulation gas cells.

slice(instrument, metal=None, norm=None, temp=None, redshift=None): Return a 1D slice of the table specified by a value in all other dimensions (only one input can, and must, remain None).

emis(instrument, metal, temp, norm=None, redshift=0.0)

Interpolate the x-ray table, return fluxes [erg/s/cm^2], luminosities [10^44 erg/s], or counts [1/s]. Input gas properties can be scalar or np.array(), in which case they must have the same size.

Parameters:

instrument (str) – name of the requested dataset.
metal (float or ndarray[float]) – metallicity [log solar].
norm (float or ndarray[float]) – usual XSPEC normalization [log cm^-5], only if use_apec == False.
temp (float or ndarray[float]) – boltzmann constant * temperature [log keV].
redshift (float)

Returns: emis (ndarray[float]): x-ray emission per cell.

calcGasEmission(sP, instrument, subhaloID=None, indRange=None, tempSfCold=True): Compute x-ray emission, either (i) flux [erg/s/cm^2], (ii) luminosity [erg/s], or (iii) counts for a particular observational setup [count/s], always linear, for gas particles in the whole snapshot, optionally restricted to an indRange. tempSfCold : set temperature of SFR>0 gas to cold phase temperature (1000 K, i.e. no x-ray emission) instead of the effective EOS temp.

plotXrayEmissivities(): Debug plots of the x-ray emissivity table trends with (z,band,T,Z).

compare_tables_single_cell(): Debugging: compare XSPEC-based and APEC_based tables, for single cell results.

compare_tables(): Debugging: compare XSPEC-based and APEC-based tables.

temet.cosmo.zooms module

Analysis and helpers specifically for zoom resimulations in cosmological volumes.

pick_halos(): Testing.

calculate_contamination(sPzoom, rVirFacs=[1, 2, 3, 4, 5, 10], verbose=False): Calculate number of low-res DM within each rVirFac*rVir distance, as well as the minimum distance to any low-res DM particle, and a radial profile of contaminating particles.

contamination_profile(): Check level of low-resolution contamination (DM particles) in zoom run. Plot radial profile.

contamination_compare_profiles(): Compare contamination radial profiles between runs.

contamination_mindist(): Plot distribution of contamination minimum distances, and trend with halo mass.

contamination_mindist2(): For a virtual box, plot contamination min dist histograms and trends.

sizefacComparison(): Compare SizeFac 2,3,4 runs (contamination and CPU times) in the testing set.

parentBoxVisualComparison(haloID, conf=0)

Make a visual comparison (density projection images) between halos in the parent box and their zoom realizations.

Parameters:

haloID (int) – the zoom halo ID, at the final redshift (z=0).
variant (str) – the zoom variant.
conf (int) – the plotting configuration.
snap (int) – if not the final snapshot, plot at some redshift other than z=0.

zoomBoxVis(sPz=None, conf=0): Make a visualization of a zoom simulation, without using/requiring group catalog information.

plot_timeevo(): Diagnostic plots: group catalog properties for one halo vs time (no merger trees).

temet.cosmo.zoomsVirtualBox module

Creation of ‘virtual (uniform) boxes’ from a number of zoom resimulations.

maskBoxRegionsByHalo(): Compute spatial volume fractions of zoom runs via discrete convex hull type approach.

combineZoomRunsIntoVirtualParentBox(snap=99): Combine a set of individual zoom simulations into a ‘virtual’ parent simulation, i.e. concatenate the output/group* and output/snap* of these runs. Process a single snapshot, since all are independent. Note that we write exactly one output groupcat file per zoom halo, and exactly two output snapshot files.

testVirtualParentBoxGroupCat(snap=99): Compare all group cat fields (1d histograms) vs TNG300-1 to check unit conversions, etc.

testVirtualParentBoxSnapshot(snap=99): Compare all snapshot fields (1d histograms) vs TNG300-1 to check unit conversions, etc.

check_groupcat_property(): Compare TNG300 vs TNG-Cluster property.

check_particle_property()

fix_tracer_offsets(): Fix Group and Subhalo TracerOffsetType fields, which had overflow in original virtual box.

Module contents

Analysis and post-processing specialized for cosmological hydrodynamical simulations.