temet.cosmo
Analysis and post-processing specialized for cosmological hydrodynamical simulations.
temet.cosmo.cloudy
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:
objectUse 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 (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, nThreads=16)
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)orx = 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]
redshift (float or None) – redshift, if we are interpolating in redshift space.
nThreads (int) – number of threads to use for interpolation (1=serial).
- Returns:
ionization fraction per cell [log].
- Return type:
ndarray
- calcGasMetalAbundances(sP, element, ionNum, indRange=None, solarAbunds=False, solarMetallicity=False, sfGasTemp='cold', parallel=False)
Compute abundance mass fraction of the given metal ion for all gas, 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.
solarAbunds (bool) – assume solar abundances (metal ratios), thereby ignoring GFM_Metals field.
solarMetallicity (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.
parallel (bool) – if True, use parallel snapshot loads (i.e. not already inside a parallelized load).
- 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:
objectUse 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 (only one input can remain None).
- Parameters:
- 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, solarAbunds=False, solarMetallicity=False, sfGasTemp='cold')
Compute luminosity of line emission for a given ‘line’, for all gas, 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.
solarAbunds (bool) – assume solar abundances (metal ratios), thereby ignoring GFM_Metals field.
solarMetallicity (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
Run grids of CLOUDY photo-ionization models for ion abundances, emissivities, or cooling rates.
- getEmissionLines()
Return the list of emission lines (
lineListabove) 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, includes 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
Galaxy clustering statistics, e.g. two-point correlation functions.
- 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 correspondingcTypeto be specified.cType (2-tuple or None) – if
colorBinis 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 correspondingmTypeto be specified.mType (str or None) – if
mstarBinis not None, then should be a string defining the stellar mass field (ofplot.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 betweenrad[i]:rad[i+1].xi_err (ndarray): one sigma errors derived from the covariance (None if
jackKnifeNumSubis None).covar (ndarray): 2d covariance matrix (None if
jackKnifeNumSubis 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 within some distance.
- Parameters:
- 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 a sample selection of: {colorBin, cType, mstarBin, mType}. Specifications for these are as in 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.
temet.cosmo.color
Calculations for optical stellar light of galaxies and galaxy colors.
- loadColors(sP, quantName)
Load 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 from snapshot or auxcats, and return magnitudes, colors, or complete data.
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 GFM_StellarPhotometrics or SubhaloStellarPhotometrics into a specified color.
Choose the right elements and handle any necessary conversions.
- calcSDSSColors(bands, redshiftRange=None, eCorrect=False, kCorrect=False, petro=False)
Load the SDSS data files and compute a requested color.
Optionally restrict to a given galaxy redshift range, correcting for extinction, and/or doing a K-correction.
temet.cosmo.cooling
Analysis and plotting related to Grackle-based cooling functions/tables.
- grackle_cooling(densities, metallicity, uvb='fg20_shielded', redshift=0.0, ssm=0, PE=False, plot=False)
Run a single-cell cooling model (vs time) with Grackle.
This will initialize a single cell at a given temperature, iterate the cooling solver for a fixed time, and output the temperature vs. time.
- grackle_equil(ssm=3)
Plot equilibrium temperature curve as a function of density (varying UVBs/SSMs), at fixed Z and z.
- grackle_equil_vs_Zz()
Plot equilibrium temperature curve as a function of density (varying Z, z) at fixed UVB/SSM.
temet.cosmo.hydrogen
The UVB and hydrogen states, including self-shielding and HI/H2 fractions (following Rahmati and Simeon Bird).
- 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.
- photoCrossSecGray(freq, J_nu, ion)
Compute a gray i.e. frequency/spectrum-averaged cross section.
- photoRate(freq, J_nu, ion)
Compute the photoionization rate of (H, HeI, HeII), and relevant photochemical rates, given the UVB (spectrum).
- 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.
- Parameters:
- Returns:
energy density U [erg/cm^3].
- Return type:
- uvbPhotoionAtten(log_hDens, log_temp, redshift)
Compute the reduction in the photoionisation rate at an energy of 13.6 eV at a given density and temp.
Uses 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.
- Parameters:
- Returns:
ratio of attenuated to unattenuated photoionisation rate. gamma_UVB_z (float): unattenuated photoionisation rate at this redshift [1/s].
- Return type:
photUVBratio (ndarray[float])
- 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.
Note: only meaningful when nH > 0.1. From Bird+ (2014) Eqn 4, e.g. the pressure-based model of Blitz & Rosolowsky (2006).
- Parameters:
nH (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.
We use the calculations of Rahmati+ (2012) for the neutral fractions as a function ofdensity.
- Parameters:
gas (dict) – gas fields (if None, will be loaded from snapshot).
sP (
simParams) – simulation instance.total (bool) – return total hydrogen mass.
totalNeutral (bool) – return total neutral hydrogen mass (HI + H2).
totalNeutralSnap (bool) – return total neutral hydrogen mass based on snapshot field.
atomic (bool) – return atomic hydrogen mass only (HI).
molecular (bool or str) – if True, return molecular hydrogen mass only (H2) using default model. If str, use specified model (e.g. ‘BL06’).
indRange (tuple) – index range to load gas cells from snapshot if gas is None.
- Returns:
hydrogen mass in the specified state for each gas cell [code units e.g. 10^10 Msun/h].
- Return type:
ndarray[float]
- calculateCDDF(N_GridVals, binMin, binMax, binSize, sP, depthFrac=1.0)
Calculate the CDDF (column density distribution function) f(N) given a set of column densities.
For example, HI or metal column densities [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].
sP (
simParams) – simulation instance.depthFrac (float) – is the fraction of sP.boxSize over which the projection was done (for dX).
temet.cosmo.kCorr
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 (str) – name of the filter to calculate K-correction for, e.g. ‘u’, ‘g’, ‘J2’.
redshifts (ndarray[float]) – redshifts of galaxies, should be between 0.0 and 0.5.
color_name (str) – name of the color, e.g. ‘u - g’, ‘g - r’, ‘ui’.
color_values (ndarray[float]) – values of galaxy color, as specified in color_name.
- Returns:
K-correction in specified filter for given redshift and color.
temet.cosmo.lightcone
Cosmological simulations - generating and working with lightcones (w/ Andres Aramburo-Garcia).
- get_cone(sP, group, config, snap_index)
Transform 3D periodic positions into lightcone geometry.
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.
- lightcone3DtoSkyCoords(pos, vel, sP, velType)
Transform 3D positions and velocities into (ra,dec,redshift) for a lightcone geometry.
pos and vel represent particle or galaxy positions in the periodic cube. The observer is assumed to be placed at the origin of the cube, (0,0,0), and the view direction is hardcoded.
temet.cosmo.mergertree
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:
- 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:
- treeFieldnames(sP, treeName='SubLink')
Load names of fields available in a mergertree.
- mpbPositionComplete(sP, id, extraFields=None)
Load a MPB that includes a complete SubhaloPos at all snapshots.
The filled version of SubhaloPos 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.
- quantMPB(sim, subhaloInd, quants, add_ghosts=False, z_vals=None, smooth=False)
Return particular quantit(ies) from a MPB.
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.
- plot_tree(sP, subhaloID, saveFilename, treeName='SubLink', dpi=100, ctName='inferno', output_fmt='png')
Visualize a full merger tree of a given subhalo.
- Parameters:
sP (
simParams) – simulation instance.subhaloID (int) – subhalo ID at sP.snap to plot the tree for.
saveFilename (str or BytesIO or None) – filename or buffer to save the plot to, or None to return image array. 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 purposes.
treeName (str) – which merger tree to use? ‘SubLink’ or ‘SubLink_gal’.
dpi (int) – resolution of the output image.
ctName (str) – colortable name for the color quantity.
output_fmt (str) – output format when saveFilename is a buffer (e.g. ‘pdf’, ‘png’, ‘jpg’).
temet.cosmo.skirt
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.
- Parameters:
- 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.
Looks for existing output files in the current working directory.
temet.cosmo.stellarPop
Stellar population synthesis, evolution, photometrics (FSPS).
- class sps(sP, iso='padova07', imf='chabrier', dustModel='cf00_res_conv', order=3, redshifted=False, emlines=False)
Bases:
objectUse 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']
- 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], possibly redshifted.
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)
Attenuate a spectrum for a given redshift.
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 luminosities in the given band, using snapshot-stored values or on-the-fly calculations.
Note that wind is returned as NaN luminosity, assuming it is filtered out elsewhere, e.g. in gridBox().
- Returns:
linear luminosities in units of [Lsun/Hz].
- Return type:
ndarray
- 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)
Spatially resolved dust attenuation on stellar spectra.
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 redshifting and apparent vs. absolute magnitudes.
Check 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_check_rawspec()
Check spectral tables.
temet.cosmo.time_evo
Efficient analysis for time evolution of (sub)halos across snapshots.
- subhalo_subbox_overlap(sP, sbNum, subInds, verbose=False)
Determine intersection with a subhalo selection and evolving tracks through a given subbox.
- halosTimeEvo(sP, haloInds, haloIndsSnap, centerPos, minSnap, maxSnap)
Derive properties for one or more halos at all subbox/normal snapshots.
Halos are defined by their evolving centerPos locations. minSnap/maxSnap define the range over which to consider each halo. sP can be a fullbox or subbox, which sets the data origin. One save file is made per halo.
- halosTimeEvoSubbox(sP, sbNum, sel, selInds)
Record several properties for one or more halos at each subbox snapshot.
Halos are specified by selInds, which index the result of subbox_subbox_overlap() which intersects the SubboxSubhaloList catalog with the simple mass selection returned by halo_selection().
- halosTimeEvoFullbox(sP, haloInds)
Record several properties for one or more halos at each full box snapshot.
Use SubLink MPB for positioning, extrapolating back to snapshot zero.
temet.cosmo.util
Helper functions related to cosmo box simulations.
- redshiftToSnapNum(redshifts=None, times=None, sP=None, recalculate=False, load=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)
Match snapshots across multiple simulations at a common set of redshifts.
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 boundary conditions for positions.
To do so, add boxSize to any negative points, subtract boxSize from any points outside box.
- correctPeriodicPosBoxWrap(vecs, sP)
Determine if an array of positions spans a periodic boundary and, if so, wrap them (for plotting).
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)
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)
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)
Compute the subhalo ID that each particle/cell belongs to.
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)
Compute the halo ID that each particle/cell belongs to.
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.
- subhaloIDsToBoundingPartIndices(sP, subhaloIDs, groups=False, strictSubhalos=False)
For a list of subhalo IDs, identfy the particle index that bounds all their members.
- Parameters:
- Returns:
- 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 as in snapshotSubset().
- Return type:
- cenSatSubhaloIndices(sP=None, gc=None, cenSatSelect=None)
Return a tuple of three sets of subhalo IDs: centrals only, centrals & satellites, and satellites only.
- subboxSubhaloCat(sP, sbNum)
Generate the SubboxSubhaloList catalog giving the intersection of fullbox subhalos with subboxes vs time.
Use the merger trees from the fullbox and interpolate positions to subbox times. Determine interpolate properties of relevant subhalos at each subbox snapshot
- subboxSubhaloCatExtend(sP, sbNum, redo=False)
Extend the SubboxSubhaloList catalog with custom (interpolated) properties.
For example, 30 pkpc stellar masses. Separated into second step since this is a heavy calculation, restartable.
- subsampleRandomSubhalos(sP, maxPerDex, mstarMinMax, mstar=None, cenOnly=False)
Sub-sample subhalos IDs randomly, uniformly across some quantity (can be generalized) such as stellar mass.
- Parameters:
sP (
simParams) – simulation instance.maxPerDex (int) – maximum number of subhalos to select per dex in the quantity.
mstarMinMax (2-tuple) – minimum and maximum values of the quantity to consider (e.g. stellar mass).
mstar (ndarray) – optional pre-loaded array of the quantity values for all subhalos..
cenOnly (bool) – if True, only consider central subhalos when mstar is not provided.
temet.cosmo.xray
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.
Interpolate 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:
objectUse 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 can 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:
- Returns:
x-ray emission per cell.
- Return type:
emis (ndarray[float])
- calcGasEmission(sP, instrument, subhaloID=None, indRange=None, tempSfCold=True)
Compute x-ray emission for all gas, optionally restricted to an index range.
Compute 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()
Plot 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
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 statistics of low-resolution dark matter particles in/near the halo.
Compute: 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 between halos in the parent box and their zoom realizations.
- 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
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.
We concatenate the output/group* and output/snap* of these runs, and 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()
Debug: compare TNG300 vs TNG-Cluster particle property.