temet.util

Utilities, units, numba-accelerated algorithms, and other algorithmic implementations.

temet.util.boxRemap

Implementation of ‘a volume and local structure preserving mapping of periodic boxes’ (Carlson, White 2010).

http://mwhite.berkeley.edu/BoxRemap/.

class Plane(p, n)

Bases: object

Plane class for half-space representation.

normal()

Normal vector.

test(x, y, z)

Point-plane comparison: >0, <0 or ==0 depending on if point is above, below, or on the plane.

class Cell(ix=0, iy=0, iz=0)

Bases: object

Cell class for holding faces/planes.

contains(x, y, z)

Check if cell contains point (x,y,z).

class Cuboid(u1=(1, 0, 0), u2=(0, 1, 0), u3=(0, 0, 1))

Bases: object

Main cuboid remapping class.

UnitCubeTest(P)

Return +1, 0, or -1 if the unit cube is above, below, or intersecting the plane.

GetCells()

Create and return flattened, ndarray representation of cells and faces/planes.

GetN123()

Return remapping numbers.

Transform(x, y, z)

Transform a point (x,y,z) in [0,1] cubical space to remapped cuboid space.

InverseTransform(r1, r2, r3)

Inverse transform a point (r1,r2,r3) in remapped cuboid space to [0,1] cubical space.

CuboidTransformPoint(x, y, z, nCellsTot, nFacesTot, faceOff, faceLen, cell_ix, cell_iy, cell_iz, face_a, face_b, face_c, face_d, n1, n2, n3)

Do Cuboid transform of a single (x,y,z) point given the flattened input configuration. Return 3-tuple.

CuboidTransformArray(pos_in, pos_out, nCellsTot, nFacesTot, faceOff, faceLen, cell_ix, cell_iy, cell_iz, face_a, face_b, face_c, face_d, n1, n2, n3)

Do Cuboid transform of a [N,3] coordinate array, given the flattened input configuration.

Write results to pos_out.

findCuboidRemapInds(remapRatio, nPixels=None)

Find the closest remapping matrix (3x3) to achieve the input remapping ratio of [x,y,z] relative extents.

Also input nPixels [x,y] of the final image to be made, in order to calculate newBoxSize.

remapPositions(sP, pos, remapRatio, nPixels)

Remap an array of points from a cubic periodic box into a remapped cuboid volume.

Input: [N,3] coordinate array from the original cubic periodic domain of side-length sP.boxSize, into a non-periodic, rectangular domain with the relative shape of remapRatio in the (x,y,z) dimensions.

temet.util.dataConvert

Various data exporters/converters, between different formats, etc.

concatSubboxFilesAndMinify()

Minify subbox snapshots, removing unwanted fields and resaving concatenated into a smaller number of chunks.

groupCutoutFromSnap(run='tng')

Create a [full] subhalo/fof cutout from a snapshot (as would be done by the Web API).

tracerCutoutFromTracerTracksCat()

Create a subhalo cutout of tracers from the full postprocessing/tracer_tracks/ catalog.

makeSnapHeadersForLHaloTree()

Copy chunk 0 of each snapshot only and the header only (for LHaloTree B-HaloTree).

makeSnapSubsetsForPostProcessing()

Copy snapshot chunks reducing to needed fields for tree/post-processing calculations.

finalizeSubfindHBTGroupCat(snap, prep=False)

Finish Subfind-HBT post-processing of an original Subfind simulation.

Extract the new group ordered IDs from the rewritten snapshot, place them into the catalog itself, and cross-match to the original snapshot IDs, to write the index mapping into the catalog. Finally, rewrite the catalog into a single file, and add fof and subhalo cross-matching to original catalog. To run: (i) first create minimal snaps with makeSnapSubsetsForPostProcessing() above, (ii) run Gadget-4 SubfindHBT, (iii) run finalizeSubfindHBTGroupCat(prep=True), (iv) create SubLinkHBT trees, (v) run makeSubgroupOffsetsIntoSublinkTreeGlobal(basePath, treeName=’SubLinkHBT’), (vi) run finalizeSubfindHBTGroupCat() to finish.

makeSdssSpecObjIDhdf5()

Transform some CSV files into a HDF5 for SDSS objid -> specobjid mapping.

createEmptyMissingGroupCatChunk()

Fix some issue in 4503 variation.

combineAuxCatSubdivisions()

Combine a subdivision of a pSplit auxCat calculation.

snapRedshiftsTxt()

Output a text-file of snapshot redshifts, etc.

tngVariantsLatexOrWikiTable(variants='all', fmt='wiki')

Output latex-syntax table describing the TNG model variations runs.

splitSingleHDF5IntoChunks(snap=151)

Split a single-file snapshot/catalog/etc HDF5 into a number of roughly equally sized chunks.

combineMultipleHDF5FilesIntoSingle()

Combine multiple groupcat file chunks into a single HDF5 file.

convertVoronoiConnectivityVPPP(stage=1, thisTask=0)

Read the Voronoi mesh data from Chris Byrohl using his vppp (voro++ parallel) approach, save to HDF5.

exportSubhalosBinary()

Export a very minimal group catalog to a flat binary format (for WebGL/Explorer3D).

exportHierarchicalBoxGrids(sP, partType='gas', partField='mass', nCells=(32, 64, 128, 256, 512), haloID=None, haloSizeRvir=2.0, retData=False, memoryReturn=False)

Export one or more 3D uniform Cartesian grids, of different resolutions, to a flat binary format.

Parameters:

  • sP (simParams) – simulation instance.

  • partType (str) – particle type, e.g. ‘gas’, ‘dm’, ‘star’.

  • partField (str) – particle field to grid, e.g. ‘mass’, ‘dens’, ‘temp’, ‘metallicity’, etc.

  • nCells (list[int]) – one or more grid sizes (number of cells per linear dimension)

  • haloID (int or None) – if None then use full box, otherwise fof-restricted particle load.

  • haloSizeRvir (float) – if haloID is specified, then this gives the box side-length in rvir units.

  • retData (bool) – if True, immediately return python object (list of ndarrays) of data.

  • memoryReturn (bool) – if True, an actual file is not written, and instead a bytes array is returned.

Returns:

Either a list, raw binary inside a io.BytesIO, or no return if an actual file is written to disk.

Note

This function is used to pre-compute the grids used in the Explorer3D WebGL volume rendering interface, as well as on-the-fly grid generation for halo-scope volume rendering.

exportIltisCutout(sim, haloIDs, emLine='O  7 21.6020A', haloSizeRvir=2.0, sfrEmisFac=1)

Export a snapshot cutout (global scope) around a given halo, for use in ILTIS-RT.

Parameters:

  • sim (simParams) – simulation instance.

  • haloIDs (list[int]) – list of FoF halo IDs to center the cutouts on.

  • emLine (str) – the emission line to save the emissivity of, which also specifies the ion for the density field.

  • haloSizeRvir (float) – the box side-length in rvir units.

  • sfrEmisFac (float) – if not unity, then multiplicative (boost) factor by which to change the emissivities for star-forming gas only, e.g. to generate simple models for bright central sources.

Returns:

None. An actual file is written to disk.

temet.util.dataConvertSim

Conversion of data (snapshots/catalogs) between different cosmological simulations.

convertGadgetICsToHDF5(aip=False)

Convert a Gadget-1/2 binary format ICs (dm-only, only pos/vel/IDs) into HDF5 format (keep original ordering).

convertMillenniumSubhaloCatalog(snap=63)

Convert a subhalo catalog (‘sub_tab_NNN.X’ files), custom binary format of Millennium to Illustris-like HDF5.

convertMillennium2SubhaloCatalog(snap=67)

Convert a subhalo catalog (‘subhalo_tab_NNN.X’ files), custom binary format of Millennium-2 to TNG-like HDF5.

convertMillenniumSnapshot(snap=63)

Convert a complete Millennium snapshot (+IDS) into Illustris-like group-ordered HDF5 format.

convertMillennium2Snapshot(snap=67)

Convert a complete Millennium-2 snapshot into TNG-like group-ordered HDF5 format.

Note all snapshots except 4-7 (inclusive) are already group-ordered.

convertEagleSnapshot(snap=20)

Convert an EAGLE simulation snapshot (HDF5) to a TNG-like snapshot (field names, units, etc).

convertSimbaSnapshot(snap=151)

Convert an SIMBA simulation snapshot (HDF5) to a TNG-like snapshot (field names, units, etc).

fixSimbaSMBHs(snap=112)

SIMBA has some SMBHs at exactly the same position, which breaks Subfind.

testSimba(redshift=0.0)

Compare all snapshot fields (1d histograms) vs TNG to check unit conversions, etc.

testSimbaCat(redshift=0.0)

Compare all group cat fields (1d histograms) to check unit conversions, etc.

temet.util.gfmExternalFiles

Create external data files required for GFM runs (photometrics, cooling, etc).

makeStellarPhotometricsHDF5_BC03()

Create stellar_photometrics.hdf5 file using BC03 models, as used for Illustris and IllustrisTNG runs.

Bands: UBVK (Buser U,B3,V,IR K filter + Palomar200 IR detectors + atmosphere.57) in Vega, griz (sdss) in AB.

Notes

temet.util.helper

General helper functions, small algorithms, basic I/O, etc.

nUnique(x)

Return number of unique elements in input numpy array x.

isUnique(x)

Does input array contain only unique values?

closest(array, value)

Return closest element of array to input value.

array_equal_nan(a, b)

As np.array_equal(a,b) but allowing NaN==NaN.

evenlySample(sequence, num, logSpace=False)

Return num samples from sequence roughly equally spaced.

contiguousIntSubsets(x)

Return a list of index pairs corresponding to contiguous integer subsets of the input array.

Final index of each pair is exclusive as in numpy syntax, so to obtain all the elements of a range from a numpy array, x[ranges[0][0]:ranges[0][1]].

logZeroSafe(x, zeroVal=1.0)

Take log10 of input variable or array, keeping zeros at some value.

logZeroMin(x)

Take log10 of input variable, setting zeros to 100 times less than the minimum.

logZeroNaN(x)

Take log10, setting zeros to NaN and leaving NaN as NaN (same as default behavior, but suppress warnings).

last_nonzero(array, axis, invalid_val=-1)

Return the indices of the last nonzero entries of the array, along the given axis.

iterable(x)

Wrap input as needed so that a for loop can iterate over this object correctly.

Protect against non-list/non-tuple (e.g. scalar or single string) values.

rebin(x, shape)

Resize input array x, must be 2D, to new shape, probably smaller, by averaging.

reportMemory()

Return current Python process memory usage in GB.

numPartToChunkLoadSize(numPart)

For a given snapshot size, in terms of total particle count, decide on a good chunk loading size.

Goal: a reasonable compute/memory balance.

tail(fileName, nLines)

Wrap linux tail command line utility.

cache(_func=None, *, overwrite=False)

Decorator to cache the return (dict) of a function. Cache filename depends on arguments.

running_median(X, Y, nBins=100, binSize=None, binSizeLg=None, skipZeros=False, percs=None, minNumPerBin=10, mean=False, weights=None)

Create a adaptive median line of a (x,y) point set using some number of bins.

running_median_clipped(X, Y_in, nBins=100, minVal=None, maxVal=None, binSize=None, skipZerosX=False, skipZerosY=False, clipPercs=(0, 90))

Create a constant-bin median line of a (x,y) point set, clipping outliers above/below clipPercs.

running_median_sub(X, Y, S, nBins=100, binSize=None, skipZeros=False, sPercs=(25, 50, 75), percs=(16, 84), minNumPerBin=10)

Create a adaptive median line of a (x,y) point set, by slicing in a third value on percentiles.

Use some number of bins, where in each bin only the sub-sample of points obtained by slicing a third value (S) above and/or below one or more percentile thresholds is used.

running_sigmawindow(X, Y, windowSize=None)

Create an local/adaptive estimate of the stddev of a (x,y) point set using a sliding window.

running_histogram(X, nBins=100, binSize=None, normFac=None, skipZeros=False)

Create a adaptive histogram of a (x) point set using some number of bins.

dist_theta_grid(size, nPixels)

Compute impact parameter and angle for every pixel of a map.

Parameters:

  • size (float) – physical size of the map [pkpc].

  • nPixels (int or 2-tuple) – number of pixels along each axis.

shrinking_center(xyz, boxSize, frac_stop=0.1, drop_frac_per_iter=0.05)

Shrinking center algorithm: iteratively search for a center position given a [N,3] coordinate set.

replicateVar(childCounts, subsetInds=None)

Given a number of children for each parent, replicate the parent indices for each child.

subset_inds : still need to walk the full child_counts, but only want parent indices of a subset.

pSplit(array, numProcs, curProc)

Divide work for embarassingly parallel problems.

pSplitRange(indrange, numProcs, curProc, inclusive=False)

As pSplit(), but accept a 2-tuple of [start,end] indices and return a new range subset.

If inclusive==True, then assume the range subset will be used e.g. as input to snapshotSubseet(), which unlike numpy convention is inclusive in the indices.

getIDIndexMapSparse(ids)

Return an array which maps ID->indices within dense, disjoint subsets which can be sparse in the entire ID range.

Within each subset i of size binsize array[ID-minids[i]+offset[i]] is the index of the original array ids where ID is found (assumes no duplicate IDs).

getIDIndexMap(ids)

Return an array such that array[ID-min(ids)] is the index of the original array ids where ID is found.

The return has size max(ids)-min(ids). Assumes a one to one mapping, not repeated indices as in the case of parentIDs for tracers).

trapsum(xin, yin)

Trapezoidal rule numerical quadrature.

leastsq_fit(func, params_init, args=None)

Wrap scipy.optimize.leastsq() by making the error function and handling returns.

If args is not None, then the standard errors are also computed, but this ASSUMES that args[0] is the x data points and args[1] is the y data points.

least_squares_fit(func, params_init, params_bounds, args=None)

Wrap least_squares() using the Trust Region Reflective method for fitting with parameter constraints.

reducedChiSq(sim_x, sim_y, data_x, data_y, data_yerr=None, data_yerr_up=None, data_yerr_down=None)

Compute reduced (i.e.) mean chi squared between a simulation ‘line’ and observed points with errors.

sgolay2d(z, window_size, order, derivative=None)

Szalay-golay filter in 2D using FFT convolutions.

kde_2d(x, y, xrange, yrange)

Simple 2D KDE calculation using scipy gaussian_kde.

gaussian_filter_nan(zz, sigma)

Filter 2D array zz by a Gaussian with input sigma, ignoring any NaN entries.

mvbe(points, tol=0.005)

Find the ‘minimum volume bounding ellipsoid’ of a given set of points (iterative).

Returns the ellipse equation in “center form” (x-c).T * A * (x-c) = 1 Based on MVEE in Matlab by Nima Moshtagh.

binned_stat_2d(x, y, c, bins, range_x, range_y, stat='median')

Replacement of binned_statistic_2d for mean_nan or median_nan.

binned_statistic_weighted(x, values, statistic, bins, weights=None, weights_w=None)

Compute a binned statistic with optional weights.

If weights == None, straight passthrough to scipy.stats.binned_statistic(). Otherwise, compute once for values*weights, again for weights alone, then normalize and return. If weights_w is not None, apply this np.where() result to the weights array.

lowess(x, y, x0, degree=1, l=1, robust=False)

Locally smoothed regression with the LOWESS algorithm (https://github.com/arokem/lowess).

x - values of x for which f(x) is input, with shape (ndim,y.size) y - values of f(x) at these points (1d) x0 - x values to output LOWESS estimate for (can be the same as x) degree - degree of smoothing functions (0=locally constant, 1=linear, …) l - metric window size for the kernel (scalar, or array of y.size) robust - if True, apply the robustification procedure from Cleveland (1979), pg 831

weighted_std_binned(x, vals, weights, bins)

Weighted standard deviation within bins.

For a given set of bins (edges), histogram x into those bins, and then compute and return the standard deviation (unbiased) of vals weighted by weights, per-bin. Assumes ‘reliability’ (non-random) weights, following https://en.wikipedia.org/wiki/Weighted_arithmetic_mean#Weighted_sample_variance.

bincount(x, dtype)

Same behavior as np.bincount() except can specify dtype different than x.dtype to save memory.

periodicDistsN(pos1, pos2, BoxSize, squared=False)

Compute periodic distance between each (x,y,z) coordinate in pos1 vs. the corresponding point in pos2.

Either pos1 and pos2 have the same shapes, and are matched pairwise, or pos1 is a tuple (i.e. reference position).

periodicDistsIndexed(pos1, pos2, indices, BoxSize)

Compute periodic distance between each (x,y,z) coordinate in pos1 vs. the corresponding point in pos2.

Here pos1.shape[0] != pos2.shape[0], e.g. in the case that pos1 are group centers, and pos2 are particle positions. Then indices gives, for each pos2, the index of the corresponding pos1 element to compute the distance to, e.g. the group ID of each particle. Return size is the length of pos2.

crossmatchHalosByCommonIDs(ids1, lengths1, ids2, lengths2)

For each object in the first set, find the best-matching object in the second set, based on common IDs.

For two sets of objects (e.g. halos or subhalos) which contain member IDs, finding the best matching index of the second set, for each object in each first set, coresponding to the object in the second set with the largest number of shared IDs.

faddeeva985(x, y)

Real arguments represent z = x + im*y and return only the real part.

xypairs_to_np(s)

Convert a string of comma-separated x,y pairs, each separated by a newline, into two separate arrays.

curRepoVersion()

Return a hash of the current state of the mercurial python repo.

temet.util.job_monitor

Helpers for working with SLURM jobs.

checkVisJobs()

Categorize a large job set into running/completed/failed and automatically re-submit jobs which have failed.

submitExpandedJobs(jobNum)

Jobs to build image pyramid for explorer.

temet.util.match

Helper function to efficiently match elements between two arrays.

match(ar1, ar2, firstSorted=False, parallel=True, debug=False)

Calculate the common elements between two arrays and their respective indices.

Returns index arrays i1,i2 of the matching elements between ar1 and ar2.While the elements of ar1 must be unique, the elements of ar2 need not be. For every matched element of ar2, the return i1 gives the index in ar1 where it can be found. For every matched element of ar1, the return i2 gives the index in ar2 where it can be found. Therefore, ar1[i1] = ar2[i2]. The order of ar2[i2] preserves the order of ar2. Therefore, if all elements of ar2 are in ar1 (e.g. ar1=all TracerIDs in snap, ar2=set of TracerIDs to locate) then ar2[i2] = ar2. The approach is one sort of ar1 followed by bisection search for each element of ar2, therefore O(N_ar1*log(N_ar1) + N_ar2*log(N_ar1)) ~= O(N_ar1*log(N_ar1)) complexity so long as N_ar2 << N_ar1.

match1(ar1, ar2, uniq=False, debug=False)

Calculate the common elements between two arrays and their respective indices (unused version 1).

My version of numpy.in1d with invert=False. Return is a ndarray of indices into ar1, corresponding to elements which exist in ar2. Meant to be used e.g. as ar1=all IDs in snapshot, and ar2=some IDs to search for, where ar2 could be e.g. ParentID from the tracers, in which case they are generally not unique (multiple tracers can exist in the same parent).

match2(ar1, ar2, debug=False)

Calculate the common elements between two arrays and their respective indices (unused version 2).

My alternative version of numpy.in1d with invert=False, which is more similar to calcMatch(). Return is two ndarrays. The first is indices into ar1, the second is indices into ar2, such that ar1[inds1] = ar2[inds2]. Both ar1 and ar2 are assumed to contain unique values. Can be used to e.g. crossmatch between two TracerID sets from different snapshots, or between some ParentIDs and ParticleIDs of other particle types. The approach is a concatenated mergesort of ar1,ar2 combined, therefore O( (N_ar1+N_ar2)*log(N_ar1+N_ar2) ) complexity.

temet.util.parallelSort

argsort(*args, **kwargs)

Overloaded function.

  1. argsort(arg0: numpy.ndarray[numpy.int64]) -> numpy.ndarray[numpy.int64]

An argsort-like function (int32 input).

  1. argsort(arg0: numpy.ndarray[numpy.int64]) -> numpy.ndarray[numpy.int64]

An argsort-like function (int64 input).

  1. argsort(arg0: numpy.ndarray[numpy.uint64]) -> numpy.ndarray[numpy.int64]

An argsort-like function (uint32 input).

  1. argsort(arg0: numpy.ndarray[numpy.uint64]) -> numpy.ndarray[numpy.int64]

An argsort-like function (uint64 input).

sort(*args, **kwargs)

Overloaded function.

  1. sort(arg0: numpy.ndarray[numpy.float32]) -> None

A function that sorts in-place (float).

  1. sort(arg0: numpy.ndarray[numpy.float64]) -> None

A function that sorts in-place (double).

  1. sort(arg0: numpy.ndarray[numpy.int16]) -> None

A function that sorts in-place (int16).

  1. sort(arg0: numpy.ndarray[numpy.int64]) -> None

A function that sorts in-place (int32).

  1. sort(arg0: numpy.ndarray[numpy.int64]) -> None

A function that sorts in-place (int64).

  1. sort(arg0: numpy.ndarray[numpy.uint16]) -> None

A function that sorts in-place (uint16).

  1. sort(arg0: numpy.ndarray[numpy.uint64]) -> None

A function that sorts in-place (uint32).

  1. sort(arg0: numpy.ndarray[numpy.uint64]) -> None

A function that sorts in-place (uint64).

temet.util.rotation

Find rotation matrices (moment of inertia tensors) to place galaxies edge-on/face-on, do coordinate rotations.

meanAngMomVector(sP, subhaloID, shPos=None, shVel=None)

Calculate the mean angular momentum 3-vector, for rotation and projection into disk face/edge-on views.

Of either the star-forming gas or the inner stellar component.

momentOfInertiaTensor(sP, gas=None, stars=None, rHalf=None, shPos=None, subhaloID=None, useStars=True, onlyStars=False)

Calculate the moment of inertia tensor (3x3 matrix) for a subhalo or halo.

Given a load of its member gas and stars (at least within 2*rHalf==shHalfMassRadStars) and center position shPos. If useStars == True, then switch to stars if not enough SFing gas present, otherwise never use stars. If onlyStars == True, use stars alone to determine.

rotationMatricesFromInertiaTensor(I)

Calculate 3x3 rotation matrix by a diagonalization of the moment of inertia tensor.

Note the resultant rotation matrices are hard-coded for projection with axes=[0,1] e.g. along z.

rotationMatrixFromVec(i_v_in, target_vec=None)

Calculate 3x3 rotation matrix to align input vec with a target vector.

By default this is the z-axis, such that with vec the angular momentum vector of the galaxy, an (x,y) projection will yield a face on view, and an (x,z) projection will yield an edge on view.

rotationMatrixFromAngleDirection(angle, direction)

Calculate 3x3 rotation matrix for input angle about an axis defined by the input direction.

Parameters:

  • angle (float) – rotation angle in degrees.

  • direction (ndarray[float][3] or list[float][3]) – 3-vector defining direction of rotation, about the origin.

rotationMatrixFromAngle(angle)

Calculate 2x2 rotation matrix for input angle, in degrees, CCW from the positive x-axis.

rotateCoordinateArray(sP, pos, rotMatrix, rotCenter, shiftBack=True)

Rotate a [N,3] array of Coordinates about rotCenter according to rotMatrix.

perspectiveProjection(n, f, l, r, b, t, pos, hsml, axes)

Transforms coordinates using a perspective projection, taking into account finite sizes.

Based on the Perspective Projection Matrix using the ratio of similar triangles (http://www.songho.ca/opengl/gl_projectionmatrix.html).

The truncated pyramid frustrum is defined by:

n (float): The distance to the near plane along the line of sight direction. f (float): The distance to the far plane along the line of sight direction. [l, r] ([float, float]): The range of “x-axis” coordinates along the near plane. [b, t] ([float, float]): The range of “y-axis” coordinates along the near plane.

(l,b)/(r,t) thus correspond to the bottom-left/top-right corners of the frustum projected onto the near plane.

The Perspective Projection Matrix computed from this set of parameters is thereafter used to transform:

pos (ndarray[float][N,3]): array of 3-coordinates for the particles; camera is situated at z=0. hsml (ndarray[float][N]): smoothing lengths. axes (list[int][2]): the axis of the projection plane, e.g. [0,1] for a projection along the z-axis.

Returns:

Transformed coordinates. Since the z-component is always mapped to the near

plane in a perspective projection, only the x and y components are transformed here. The z component is left as is for use in, e.g., sphMap(). As a result, the transformed coordinates cannot be unprojected.

tHsml (ndarray[float][N]): Transformed hsml; equal to original hsml along the near plane, and scales

inversely with projection distance farther away from camera.

Return type:

tPos (ndarray[float][N,3])

Notes

  • The tranformed values of hsml will be negative if the point is behind the camera.

ellipsoidfit(pos_in, mass, scalerad, rin, rout, weighted=False)

Iterative algorithm to derive the shape (axes ratios) of a set of particles/cells.

Their distribution is given by pos_in and mass, within a radial shell defined by rin, rout. Positions should be unitless, normalized by scalerad (a scaling factor, could be e.g. r200 or half mass radius). Originally from Eddie Chua.

temet.util.simParams

The simParams class encapsulates all meta-data and details for a particular simulation.

In addition, it can hold specifications to point to a specific part of a given simulation, for instance a particular snapshot/redshift or halo/subhalo. We often start analysis by creating an instance of this class as:

sim = temet.sim(run='tng100-1', redshift=0.0)

To analyze a new or custom simulation which is not defined in this file, you can instead pass its path, and the simulation metadata will be automatically loaded.

sim = temet.sim('/virgo/simulations/IllustrisTNG/TNG50-1/', redshift=0.0)

This sim object can then be passed to many analysis routines, which then automatically know (i) where to find the data files, and (ii) which snapshot to analyze. Furthermore, many of the most common data loading functions are “attached” to this sim object, such that the two follow calls are functionally identical:

masses = load.snapshot.snapshotSubset(sim, 'gas', 'mass')
masses = sim.snapshotSubset('gas', 'mass')

The second providing a convenient short-hand. If you were using the public data release scripts alone, these two calls would also be identical to:

import illustris_python as il
basePath = 'sims.TNG/TNG100-1/output/'
snap = 99
masses = il.snapshot.loadSubset(basePath, snap, 'gas', fields=['Masses'])

Note that upon creation, each simParams object instantiates and fills a util.units class which is automatically tailored to this simulation and its unit system. For example

dists = sim.snapshotSubset('gas', 'rad', haloID=1234)
dists_kpc = sim.units.codeLengthToKpc(dists)

converts the original dists array from code units to physical kpc. While “code units” typically implies ckpc/h, this particular simultion could be instead using a Mpc-unit system, and this difference is handled transparently, as are different values of h for different simulations.

To add a new simulation (suite), one of the existing examples should be copied and adapted as needed.

class simParams(run, res=None, variant=None, redshift=None, time=None, snap=None, hInd=None, haloInd=None, subhaloInd=None, arepoPath=None, simName=None)

Bases: object

Simulation instance: metadata, helper functions, and (optionally) a snapshot/redshift specification.

simPath = ''

path (root) containing ‘output’ directory with simulation snapshots and group catalogs

arepoPath = ''

path to Arepo binary, Config.sh and param.txt files of this run

derivPath = ''

path to put derivative files (“data.files/”)

cachePath = ''

path to put cache files (“data.files/cache/”)

postPath = ''

path to put postprocessed files (“postprocessing/”)

simName = ''

label to add to plot legends (e.g. “Illustris-2”, “TNG300-1”)

simNameAlt = ''

alternative label for simulation names (e.g. “L75n1820FP”, “L205n1250TNG_DM”)

run = ''
res = 0
variant = ''
groupOrdered = None
comoving = True
boxSize = 0.0
targetGasMass = 0.0
gravSoft = 0.0
omega_m = None
omega_L = None
omega_k = 0.0
omega_b = None
HubbleParam = None
mpcUnits = False
nTypes = 6
numSnaps = 0
subbox = None
subboxCen = None
subboxSize = None
levelmin = 0
levelmax = 0
zoomLevel = 0
hInd = None
zoomShift = None
zoomShiftPhys = None
sP_parent = None
trMassConst = 0.0
trMCPerCell = 0
trMCFields = None
trVelPerCell = 0
haloInd = None
subhaloInd = None
refPos = None
refVel = None
colors = None
marker = None
data = None
metals = None
eEOS = False
star = False
BHs = False
winds = False
redshift = None
time = None
snap = None
scan_simulation(inputPath, simName=None)

Fill simulation parameters automatically, based on path.

lookup_simulation(res=None, run=None, variant=None, redshift=None, time=None, snap=None, hInd=None, haloInd=None, subhaloInd=None)

Fill parameters based on inputs (hardcoded).

auxCatSplit(field, nThreads=1, **kwargs)

Automatically do a pSplit auxCat calculation.

setRedshift(redshift=None)

Update sP based on new redshift.

setSnap(snap=None)

Update sP based on new snapshot.

ptNum(partType)

Return particle type number (in snapshots) for input partType string.

Allows different simulations to use arbitrary numbers for each type (so far they do not).

isPartType(ptToCheck, ptToCheckAgainst)

Check whether the two particle types are the same.

Return either True or False depending on if ptToCheck is the same particle type as ptToCheckAgainst. For example, ptToCheck could be ‘star’, ‘stars’, or ‘stellar’ and ptToCheckAgainst could then be e.g. ‘stars’. The whole point is to remove any hard-coded dependencies on numeric particle types. Otherwise, you would naively check that e.g. partTypeNum(ptToCheck)==4. This can now vary for different simulations (so far it does not).

copy()

Return a deep copy of this simparams, which can be manipulated/changed without affecting the original.

property name

Simulation name.

property isZoom

Is a zoom simulation?

property isZoomOrVirtualBox

Is a zoom simulation or TNG-Cluster (which has dmlowres particles in full box)?

property isDMO

Is a dark-matter-only simulation?

property isSubbox

Is a subbox ‘simulation’ i.e. variant?

property hasMergerTree

Has the default merger tree (SubLink) available?

property partTypes

Return a list of particle type names contained in this simulation, excluding tracers.

property parentBox

Return a sP corresponding to the parent volume, at the same redshift (fullbox for subbox only for now).

property dmoBox

Return a sP corresponding to the DMO/Dark analog volume, at the same redshift.

subboxSim(sbNum)

Return a sP corresponding to the given subbox number, at the same redshift.

property numMetals

Return number of tracked metal species in this simulation.

property scalefac

Current cosmological scale factor.

property tage

Current age of the universe [Gyr].

property tlookback

Current lookback time [Gyr].

property boxSizeCubicPhysicalMpc

Return box size in cubic physical Mpc at the current redshift.

property boxSizeCubicComovingMpc

Return box size in cubic comoving Mpc.

property boxLengthDeltaRedshift

The redshift interval dz corresponding to traversing one box side-length.

property dz

The redshift interval dz corresponding to traversing one box side-length.

property boxLengthComovingPathLength

The comoving pathlength i.e. ‘absorption distance’ (dX) corresponding to one box side-length.

property dX

The comoving pathlength i.e. ‘absorption distance’ (dX) corresponding to one box side-length.

property zoomSubhaloID

Return the subhalo ID of the target subhalo in the zoom suite. In general, not known, and may not be zero.

property dmParticleMass

Return dark matter particle mass (scalar constant) in code units.

property numHalos

Return number of FoF halos / groups in the group catalog at this sP.snap.

property numSubhalos

Return number of Subfind subhalos in the group catalog at this sP.snap.

property subhaloIndsCen

Return indices of central subhalos in the group catalog at this sP.snap.

property subhaloIndsSat

Return indices of satellite subhalos in the group catalog at this sP.snap.

property numPart

Return number of particles/cells of all types at this sP.snap.

property simCode

Return simulation code used to produce snapshots. Currently only differentiates between AREPO and SWIFT.

property cpuHours

Return CPU core hours to z=0 for this simulation.

property config

Return Config.sh variable dict (from snapshot metadata).

property params

Return param.txt variable dict (from snapshot metadata).

temet.util.sphMap

Interpolation of scattered point sets onto a uniform grid using the SPH spline kernel deposition method.

sphMap(pos, hsml, mass, quant, axes, boxSizeImg, boxSizeSim, boxCen, nPixels, ndims, hsml_1=None, colDens=False, nThreads=16, multi=False, posTarget=None, maxIntProj=False, minIntProj=False, refGrid=None)

Calculate gridded maps of scattered pointsets using the SPH spline kernel.

Simultaneously calculate a gridded map of projected mass and some other mass weighted quantity (e.g. temperature) with the sph spline kernel. If quant=None, the map of mass is returned, optionally converted to a column density map if colDens=True. If quant is specified, the mass-weighted map of the quantity is additionally returned. If len(nPixels) is three, compute grid and not image.

Parameters:

  • pos (ndarray[float][N,3 or N,2]) – array of 3-coordinates for the particles (or 2-coords only, to ignore z-axis).

  • hsml (ndarray[float][N]) – smoothing lengths to use for the particles, same units as pos.

  • mass (ndarray[float][N or 1]) – mass to deposit, can be scalar value if all particles have the same mass.

  • quant (ndarray[float][N]) – some quantity to calculate a mass-weighted projection of (None=skip).

  • axes (list[int][2]) – two axis indices defining the projection direction, e.g. [0,1] for x,y (project along z-axis).

  • boxSizeImg (list[float][3]) – the physical size the image should cover, same units as pos.

  • boxSizeSim (list[float][3]) – the physical size of the simulation box for periodic wrapping (0=non periodic).

  • boxCen (list[float][3]) – (x,y,z) coordinates of the center of the imaging box, same units as pos.

  • nPixels (list[int][2 or 3]) – number of pixels in x,y directions for output image or number of cells in x,y,z directions for output grid.

  • ndims (int) – number of dimensions of simulation (1,2,3), to set SPH kernel coefficients.

  • hsml_1 (ndarray[float][N]) – if not None, do asymmetric sph kernel projection with hsml==hsml_0 (x), and hsml_1 (y).

  • colDens (bool) – if True, normalize each grid value by its area (default=False).

  • nThreads (int) – do multithreaded calculation (mem required=nThreads times more).

  • multi (bool) – if True, return the tuple (dens,quant) instead of selecting one, under the assumption that we are using sphMap() in a chunked mode and have to normalize later.

  • posTarget (ndarray[float][N,3]) – array of 3-coordinates for a set of targets (e.g. star particle locations). If not None, then both pos and posTarget will be sorted on the 3rd (projection direction) coordinate, and mapping will be front to back, such that the grid value at the position (2D+depth) of each posTarget is estimated. In this case, the return is of shape [N(posTarget)] giving the projected density or mass weighted quantity along the line of sight to each posTarget.

  • maxIntProj (bool) – perform a maximum intensity projection (MIP) instead of the usual weighting.

  • minIntProj (bool) – perform a minimum intensity projection, instead of the usual weighting.

  • refGrid (ndarray[float][nPixels[0],nPixels[1]]) – if not None, calculation is the sum of squared differences between the mass-weighted particle quantity and the value of this reference grid at its pixel location (i.e. to calculate standard deviation per pixel).

Returns:

computed 2D or 3D array, depositing either mass or mass-weighted quant. If multi is True, then a 2-tuple of the (dens_map, quant_map) arrays separately.

Return type:

ndarray[float][nPixels]

Notes

  • the transpose of _calcSphMap() is taken such that with default plotting approaches e.g. axes=[0,1] gives imshow(return[i,j]) with x and y axes corresponding correctly to code coordinates. This does not apply if len(nPixels) == 3 and a 3D grid is computed. No transpose is taken in that case.

  • the entries of boxSizeImg, boxCenter, and boxSizeSim correspond to [axes[0],axes[1],axes[2]], meaning pos[:,axes[0]] is compared against the first entries boxSizeImg[0] and boxCenter[0] and not compared against e.g. boxSizeImg[axes[0]].

sphMapWholeBox(pos, hsml, mass, quant, axes, nPixels, sP, colDens=False, nThreads=8, multi=False, posTarget=None, sliceFac=1.0)

Wrap sphMap() specialized to projecting an entire cosmological/periodic box (or subbox).

Specifically, (ndims,boxSizeSim,boxSizeImg,boxCen) are derived from sP, and nPixels should be input as a single scalar which is then assumed to be square.

sphGridWholeBox(sP, pos, hsml, mass, quant, nCells=32)

Wrap sphMap(), produce a 3D grid of the given pt/field combination for an entire box.

temet.util.subfind

Implementation of (serial) subfind algorithm.

load_custom_dump(sP, GrNr, phase1=False)

Load groups_{snapNum}/fof_{fofID}.{taskNum} custom binary data dump.

load_snapshot_data(sP, GrNr)

Load the FoF particles from an actual snapshot for testing.

treeSearchIndicesIterate(P, xyz, h_guess, nNGB, boxSizeSim, next_node, tree_nodes)

Helper routine for subfind(), see below.

Note: no nNGBDev, instead we terminate if we ever find >=nNGB, and the return is sorted by distance.

buildFullTree(P, boxSizeSim, xyzMin, xyzMax, extent, ForceSoftening)

As above, but minimal and JITed.

subfind_treeevaluate_potential(target, P, ForceSoftening, next_node, tree_nodes, boxHalf, boxSizeSim)

Evaluate gravitational potential using the tree.

subfind_unbind_calculate_potential(num, ud, loc_P, ForceSoftening, NextNode, TreeNodes, boxhalf, boxsize, PS, G, atime)

Loop to parallelize.

subfind_unbind_calculate_potential_weak(num, ud, loc_P, ForceSoftening, NextNode, TreeNodes, boxhalf, boxsize, PS, G, atime, weakly_bound_limit)

Loop to parallelize.

subfind_unbind(P, SphP, PS, ud, num, vel_to_phys, H_of_a, G, atime, boxsize, SofteningTable, ForceSoftening, xyzMin, xyzMax, extent, central_flag)

Unbinding.

subfind(P, PS, SphP, StarP, BHP, atime, H_of_a, G, boxsize, SofteningTable, ForceSoftening)

Run serial subfind. (Offs = 0).

subfind_properties_1(candidates)

Determine the parent subhalo for each candidate.

get_time_difference_in_Gyr(a0, a1)

Time difference between two cosmological scale factors (from timestep.c).

assign_stellar_photometrics(i, P, StarP, mags, LogMetallicity_bins, LogAgeInGyr_bins, TableMags, atime)

Interpolate for stellar magnitudes (from GFM/stellar_photometrics.c).

subfind_properties_2(candidates, Tail, Next, P, PS, SphP, StarP, BHP, atime, H_of_a, G, boxsize, LogMetallicity_bins, LogAgeInGyr_bins, TableMags, SofteningTable, GrNr)

Determine the properties of each subhalo.

subfind_particle_order(P, PS)

Prepare a particle-level output order according to our new PS (subgroup) assignment.

set_softenings(P, SphP, sP, snapLoaded=False)

Generate SofteningTable following grav_softening.c, and set P[].SofteningType values (based on sizes/masses).

load_gfm_stellar_photometrics()

Load photometrics table for interpolation later.

run_subfind_snapshot(sP, GrNr)

Run complete Subfind algorithm on a single group, using the saved snaphshot itself.

run_subfind_customfof0save_phase1(sP, GrNr=0)

Run complete Subfind algorithm on custom FOF0 save files, phase one (TNG50-1).

run_subfind_customfof0save_phase2(sP, GrNr=0)

Run complete Subfind algorithm on custom FOF0 save files, phase two (TNG50-1).

verify_results(sP, GrNr=0)

Check results vs. actual group catalog for a test run where the subgroups of FOF0 are computed in AREPO.

rewrite_groupcat(sP, GrNr=0)

Rewrite a group catalog which is missing FOF0 subhalos using the phase2 (final) results.

add_so_quantities(sP, GrNr=0)

FoF0 is missing SO quantities (Group_R_Crit200, etc). Need to derive now.

rewrite_snapshot(sP, GrNr=0)

Rewrite a snapshot which is missing FOF0 subhalos using the phase2 (final) results.

rewrite_particle_level_cat(sP, filename, partType)

Shuffle the FoF0 segment of a dataset from a full snapshot into the new order.

compare_subhalos_all_quantities(snap_start=67)

Plot diagnostic histograms.

run_subfind(snap)

Main driver.

benchmark()

Benchmark.

temet.util.tpcf

Two-point correlation functions (pairwise distances).

tpcf(pos, radialBins, boxSizeSim, weights=None, pos2=None, weights2=None, nThreads=32)

Calculate and simultaneously histogram the results of a two-point auto correlation function.

Approach: compute all the pairwise (periodic) distances in pos. 3D only.

pos[N,3] : array of 3-coordinates for the galaxies/points radialBins[M] : array of (inner) bin edges in radial distance (code units) boxSizeSim[1] : the physical size of the simulation box for periodic wrapping (0=non periodic) weights[N] : if not None, then use these scalars for a weighted correlation function pos2[L,3] : if not None, then cross-correlation instead of auto, of pos vs. pos1 weights2[L] : must be None if weights is None, and vice versa

return is xi(r),DD,RR where xi[i]=(DD/RR-1) is computed between radialBins[i:i+1] (size == M-1)

quantReductionInRad(pos_search, pos_target, radial_bins, quants, reduce_op, boxSizeSim, nThreads=8)

Calculate a reduction operation on one or more quantities for each search point.

In each case, consider all target points falling within a 3D periodic search radius.

pos_search[N,3] : array of 3-coordinates for the galaxies/points to search from pos_target[M,3] : array of the 3-coordiantes of the galaxies/points to search over radial_bins[M] : array of bin edges in radial distance (code units) quants[M]/[M,P] : 1d or P-d array of quantities, one per pos_target, to process reduce_op[str] : one of ‘min’, ‘max’, ‘sum’ boxSizeSim[1] : the physical size of the simulation box for periodic wrapping (0=non periodic)

return is reduced_quants[N,M-1]/[N,M-1,P]

temet.util.treeSearch

Adaptive estimation of a smoothing length (radius of sphere enclosing N nearest neighbors) using oct-tree.

buildFullTree(pos, boxSizeSim, treePrec=None, verbose=False)

Helper. See below.

calcHsml(pos, boxSizeSim, posSearch=None, posMask=None, nNGB=32, nNGBDev=1, nDims=3, weighted_num=False, treePrec='single', tree=None, nThreads=32, nearest=False, verbose=False)

Calculate a characteristic size (i.e. smoothing length) given a set of input particle coordinates.

The size is defined as the radius of the sphere (or circle in 2D) enclosing the nNGB nearest neighbors. If posSearch==None, then pos defines both the neighbor and search point sets, otherwise a radius for each element of posSearch is calculated by searching for nearby points in pos.

Parameters:

  • pos (ndarray[float][N,3]/[N,2]) – array of 3-coordinates for the particles (or 2-coords for 2D).

  • boxSizeSim (float) – the physical size of the simulation box for periodic wrapping (0=non periodic).

  • posSearch (ndarray[float][M,3]) – search sites.

  • posMask (ndarray[bool][N]) – if not None, bool mask, only True entries are considered in search.

  • nNGB (int) – number of nearest neighbors to search for in order to define HSML.

  • nNGBDev (int) – allowed deviation (+/-) from the requested number of neighbors.

  • nDims (int) – number of dimensions of simulation (1,2,3), to set SPH kernel coefficients.

  • weighted_num (bool) – if True, use the SPH kernel weighted number of neighbors, instead of real number.

  • treePrec (str) – construct the tree using ‘single’ or ‘double’ precision for coordinates.

  • tree (list or None) – if not None, should be a list of all the needed tree arrays (pre-computed), i.e the exact return of :py:func`~util.treeSearch.buildFullTree`.

  • nThreads (int) – do multithreaded calculation (on treefind, while tree construction remains serial).

  • nearest (bool) – if True, then instead of returning hsml values based on nNGB, return indices and distances to single closest match only.

  • verbose (bool) – print timing and warning messages.

Returns:

derived smoothing length for each input point (if not nearest). Instead if nearest == True, then a 2-tuple of (dists,indices).

Return type:

ndarray[float or int]

calcQuantReduction(pos, quant, hsml, op, boxSizeSim, posSearch=None, treePrec='single', tree=None, nThreads=32)

Reduction operation on all particles/cells within a fixed search distance around each position.

Calculate a reduction of a given quantity on all the particles within a fixed search distance hsml around each pos. If posSearch==None, then pos defines both the neighbor and search point sets, otherwise a reduction at the location of each posSearch is calculated by searching for nearby points in pos.

Parameters:

  • pos (ndarray[float][N,3]/[N,2]) – array of 3-coordinates for the particles (or 2-coords for 2D).

  • quant (ndarray[float][N]) – array of quantity values (i.e. mass, temperature).

  • hsml (ndarray[float][1 or N]) – array of search distances, or scalar value if constant.

  • op (str) – ‘min’, ‘max’, ‘mean’, ‘kernel_mean’, ‘sum’, ‘count’.

  • boxSizeSim (float) – the physical size of the simulation box for periodic wrapping (0=non periodic).

  • posSearch (ndarray[float][N,3]) – search coordinates (optional).

  • tree (list or None) – if not None, should be a list of all the needed tree arrays (pre-computed), i.e the exact return of buildFullTree().

  • treePrec (str) – construct the tree using ‘single’ or ‘double’ precision for coordinates.

  • nThreads (int) – do multithreaded calculation (on treefind, while tree construction remains serial).

Returns:

reduction operation applied for each input.

Return type:

ndarray[float]

calcParticleIndices(pos, posSearch, hsmlSearch, boxSizeSim, posMask=None, treePrec='single', tree=None)

Find and return the actual particle indices (indexing pos, hsml) within the search radii of posSearch locations.

Serial by construction, since we do only one search.

Parameters:

  • pos (ndarray[float][N,3]/[N,2]) – array of 3-coordinates for the particles (or 2-coords for 2D).

  • posSearch (ndarray[float][3]) – search position.

  • hsmlSearch (float) – search distance.

  • boxSizeSim (float) – the physical size of the simulation box for periodic wrapping (0=non periodic).

  • posMask (ndarray[bool][N]) – if not None, then only True entries are considered in the search.

  • treePrec (str) – construct the tree using ‘single’ or ‘double’ precision for coordinates.

  • tree (list or None) – if not None, should be a list of all the needed tree arrays (pre-computed), i.e the exact return of buildFullTree().

Returns:

list of indices into pos of the neighbors within the search distance of the search position.

Return type:

ndarray[int]

temet.util.units

Various unit conversion and cosmological calculation utilities.

Includes conversion from ‘code units’ to physical units, conversions between different unit systems, and derivations of common halo properties (e.g. virial temperature) and cosmological quantities (e.g. age at a given redshift).

class units(sP=None)

Bases: object

Contains static methods which perform various unit conversions.

Can also be instantiated with a redshift/sP, in which case contains the relevant unit system and redshift-dependent constants.

UnitMass_in_g = 1.989e+43
UnitVelocity_in_cm_per_s = 100000.0
UnitVelocity_str = 'km/s'
boltzmann = 1.38065e-16
boltzmann_JK = 1.38065e-23
boltzmann_keV = 11604505.0
planck_erg_s = 6.626e-27
planck_eV_s = 4.135668e-15
hc_kev_ang = 12.3981689
mass_proton = 1.672622e-24
mass_electron = 9.1095e-28
gamma = 1.666666667
hydrogen_massfrac = 0.76
helium_massfrac = 0.24
mu = 0.6
Gravity = 6.6738e-08
H0_h1_s = 3.24078e-18
Z_solar = 0.0127
L_sun = 3.839e+33
Msun_in_g = 1.98892e+33
c_cgs = 29979000000.0
c_km_s = 299790.0
c_kpc_Gyr = 306595.0
sigma_thomson = 6.6524e-25
electron_charge = 4.8032e-10
rydberg_ang = 0.00109737
rydberg_freq = 3289840000000000.0
Habing = 5.29e-14
CourantFac = 0.3
BH_eps_r = 0.2
BH_eps_f_high = 0.1
BH_f_thresh = 0.05
N_SNII = 0.0118008
E_SNII51 = 1.0
winds_tau = 0.1
winds_fZ = 0.25
winds_Zref = 0.002
winds_gamma_Z = 2.0
winds_e = 3.6
winds_kappa = 7.4
winds_vmin = 350.0
PhysDensThresh = 0.000754654
sh03_nH_thresh = 0.13
sh03_A0 = 573.0
sh03_T_SN = 57300000.0
sh03_T_c = 1000.0
sh03_beta = 0.22578
sh03_t_star = 3.27665
s_in_yr = 31556930.0
pc_in_cm = 3.08568e+18
km_in_cm = 100000.0
kpc_in_ly = 3261.56
arcsec_in_rad = 4.84814e-06
ster_in_arcsec2 = 42500000000.0
ang_in_cm = 1e-08
erg_in_J = 1e-07
erg_in_kev = 624150000.0
eV_in_erg = 1.602e-12
UnitLength_in_cm = 3.085678e+21
UnitLength_str = 'ckpc/h'
UnitMass_str = '10$^{10}$ M$_{\\rm sun}$/h'
UnitTime_in_s = None
UnitDensity_in_cgs = None
UnitPressure_in_cgs = None
UnitEnergy_in_cgs = None
UnitTemp_in_cgs = None
UnitMass_in_Msun = None
UnitTime_in_yr = None
s_in_kyr = None
s_in_Myr = None
s_in_Gyr = None
kpc_in_km = None
Mpc_in_cm = None
kmS_in_kpcYr = None
kmS_in_kpcGyr = None
kpc_in_cm = None
msunKpc3_in_gCm3 = None
H0 = None
G = None
rhoCrit = None
rhoCrit_msunMpc3 = None
H0_kmsMpc = None
mag2cgs = None
c_ang_per_sec = None
f_b = None
Hubble = None
rhoBack = None
H2_z_fact = None
Omega_z = None
H_z = None
H_of_a = None
rhoCrit_z = None
scalefac = None
codeMassToMsun(mass)

Convert mass from code units (10**10 msun/h) to (msun).

codeMassToLogMsun(mass)

Convert mass from code units (10**10 msun/h) to (log msun).

msunToCodeMass(mass)

Convert mass in [msun] to code units.

logMsunToCodeMass(mass)

Convert mass in (log msun) to code units.

codeMassToVirTemp(mass, meanmolwt=None, log=False)

Convert from halo mass in code units to virial temperature in Kelvin, at the specified redshift.

(Barkana & Loeb (2001) eqn.26).

codeMassOverTimeToMsunPerYear(mass_over_time_val)

Convert a code [mass/time] value (e.g. BH_Mdot) into [Msun/yr]. The usual 10.22 factor.

codeBHMassToMdotEdd(mass, eps_r=None)

Convert a code mass (of a blackhole) into dM/dt_eddington in [Msun/yr].

Also available directly as ‘BH_MdotEddington’ field.

codeBHMassToLumEdd(mass)

Convert a code mass (of a blackhole) into its Eddington luminosity [erg/s].

codeBHMassMdotToBolLum(mass, mdot, basic_model=False, obscuration=False)

Convert a (BH mass, BH mdot) pair to a bolometric luminosity [erg/s].

BH_chi(M_BH)

Return chi(M_BH) threshold for the fiducial Illustris/TNG model parameters. M_BH in Msun.

codeBHValsToFeedbackMode(bh_mass, bh_mdot, bh_mdot_bondi, bh_mdot_edd)

Return the feedback mode (0 = low/kinetic, 1 = high/quasar) of the BH, based on its mass and mdot.

codeBHMassMdotToInstantaneousEnergy(bh_mass, bh_mdot, bh_density, bh_mdot_bondi, bh_mdot_edd)

Convert the instantaneous mass/mdot of a BH to the instantaneous energy injection rate [erg/s].

All inputs are in code units. The energy being injected follows the (mode-dependent) feedback model of the simulation.

codeMetallicityToWindSpecificEnergy(metal_code)

Convert the metallicity of a gas cell into the wind specific energy parameter (see TNG methods).

codeSfrZToWindEnergyRate(sfr_msunyr, metal_code)

Convert the SFR [Msun/yr] of a gas cell into dot{E}_SN of total wind energy available [10^51 erg/s].

codeSfrZToWindMomentumRate(sfr_msunyr, metal_code, dm_veldisp)

Convert the SFR [Msun/yr] of a gas cell into dot{p}_SN of total wind momentum available [10^51 g*cm/s^2].

sigmaDMToWindVel(dm_veldisp)

Convert a code (3D) velocity dispersion (i.e. SubfindVelDisp), into the wind launch velocity [km/s].

codeSfrZSigmaDMToWindMassLoading(sfr_msunyr, metal_code, dm_veldisp)

Convert a gas cell SFR [Msun/yr], metallicity [code], and 3D vel disp [km/s] into the wind mass loading.

The return is dimensionless linear.

densToSH03TwoPhase(dens_in, sfr)

Convert a gas cell density to values corresponding to the two-phase state of the sub-cell gas.

According to the Springel & Hernquist (2003) sub-grid ISM pressurization model.

Parameters:

  • dens_in (float or ndarray) – hydrogen physical gas number density [H atoms/cm^3]

  • sfr (float or ndarray) – star formation rate [Msun/yr]

Returns: a 2-tuple composed of

  • x (ndarray): the cold-phase (i.e. cold clouds) mass fraction, which is density-dependent, and always between zero and unity.

  • T_h (ndarray): the hot-phase temperature, which is density-dependent, and increases from ~1e5 to ~1e8 K.

logMsunToVirTemp(mass, meanmolwt=None, log=False)

Convert halo mass (in log msun, no little h) to virial temperature at specified redshift.

codeLengthToComovingKpc(x)

Convert length/distance in code units to comoving kpc.

codeLengthToKpc(x)

Convert length/distance in code units to physical kpc.

codeLengthToMpc(x)

Convert length/distance in code units to physical Mpc.

codeLengthToPc(x)

Convert length/distance in code units to physical parsec.

codeLengthToComovingMpc(x)

Convert length/distance in code units to comoving Mpc.

codeLengthToCm(x)

Convert length/distance in code units to cgs (cm).

codeAreaToKpc2(x)

Convert an area [length^2] in code units to physical kpc^2.

codeAreaToMpc2(x)

Convert an area [length^2] in code units to physical Mpc^2.

codeVolumeToCm3(x)

Convert a volume [length^3] in code units to physical cm^3 (cgs).

codeVolumeToKpc3(x)

Convert a volume [length^3] in code units to physical kpc^3.

codeVolumeToMpc3(x)

Convert a volume [length^3] in code units to physical Mpc^3.

physicalKpcToCodeLength(x)

Convert a length in [pkpc] to code units [typically ckpc/h].

lightyearsToCodeLength(x)

Convert a length in [lightyears] to code units.

physicalMpcToCodeLength(x)

Convert a length in [pMpc] to code units [typically ckpc/h].

particleCodeVelocityToKms(x)

Convert velocity field (for cells/particles, not group properties) into km/s.

groupCodeVelocityToKms(x)

Convert velocity vector (for groups, not subhalos nor particles) into km/s.

subhaloCodeVelocityToKms(x)

Convert velocity vector (for subhalos, not groups nor particles) into km/s.

subhaloSpinToKpcKms(x)

Convert spin vector (for subhalos, not groups nor particles) into kpc km/s.

particleCodeBFieldToGauss(b)

Convert magnetic field 3-vector (for cells) into Gauss, input b is PartType0/MagneticField.

particleCodeDivBToGaussPerKpc(divb)

Convert magnetic field divergence into [Gauss/kpc] physical, input is PartType0/MagneticFieldDivergence.

codePotentialToEscapeVelKms(pot)

Convert Potential [(km/s)^2/a] into an escape velocity [km/s].

particleAngMomVecInKpcKmS(pos, vel, mass, haloPos, haloVel)

Calculate particle angular momentum 3-vector in [Msun*kpc km/s].

Takes input arrays of pos,vel,mass and the halo CM position and velocity to compute relative to. Includes Hubble correction.

particleSpecAngMomMagInKpcKmS(pos, vel, mass, haloPos, haloVel, log=False)

Calculate particle specific angular momentum magnitude in [kpc km/s].

particleRadialVelInKmS(pos, vel, haloPos, subhaloVel)

Calculate particle radial velocity in [km/s] (negative=inwards).

Takes input arrays of pos,vel and the halo CM position and velocity to compute relative to. Includes Hubble correction.

particleRelativeVelInKmS(vel, subhaloVel)

Calculate particle velocity magnitude in [km/s], relative to a given reference frame motion.

This reference is typically SubhaloVel. If not, must be in physical units.

codeDensToPhys(dens, scalefac=None, cgs=False, numDens=False, msunpc3=False, totKpc3=False)

Convert mass density comoving->physical and add little_h factors.

Unless overridden by a parameter option, the default return units are \([10^{10} M_\odot/\rm{kpc}^3]\).

Parameters:

  • dens (array[float]) – density in code units, should be \([10^{10} M_\odot/h / (ckpc/h)^3]\) = \([10^{10} M_\odot h^2 / ckpc^3]\).

  • scalefac (array[float]) – if provided, use this scale factor instead of the current snapshot value. Can be a single value, or a different value per density.

  • cgs (bool) – if True, return units are [g/cm^3].

  • numDens (bool) – if True and cgs == True, return units are [1/cm^3].

  • msunpc3 (bool) – if True, return units are [Msun/pc^3].

  • totKpc3 (bool) – if True, return units are [[orig units]/kpc^3].

Returns:

densities in physical units, as specified above.

Return type:

array[float]

physicalDensToCode(dens, cgs=False, numDens=False)

Convert mass density in physical units to code units (comoving + w/ little h factors, in unit system).

Input: dens in [msun/kpc^3] or [g/cm^3 if cgs==True] or [1/cm^3 if cgs==True and numDens==True]. Output: dens in [10^10 Msun/h / (ckpc/h)^3] = [10^10 Msun h^2 / ckpc^3] comoving.

codeColDensToPhys(colDens, cgs=False, numDens=False, msunKpc2=False, totKpc2=False)

Convert a mass column density [mass/area] from comoving -> physical and remove little_h factors.

Unless overridden by a parameter option, the default return units are \([10^{10} M_{\odot}/kpc^2]\).

Parameters:

  • colDens (ndarray) – column densities in code units, which should be \([10^{10} M_{\odot}/h / (ckpc/h)^2]\) = \([10^{10} M_{\odot} * h / ckpc^2].\)

  • cgs (bool) – if True, return units in [g/cm^2].

  • numDens (bool) – if True and cgs == True, return units in [1/cm^2], which is in fact [H atoms/cm^2].

  • msunKpc2 (bool) – return units in [Msun / kpc^2].

  • totKpc2 (bool) – return units in [[orig units] / kpc^2].

Returns:

physical column densities, units depending on the above.

Return type:

ndarray[float32]

UToTemp(u, xe, log=False)

Convert (U,xe) pair in code units to temperature in Kelvin.

TempToU(temp, xe, log=False)

Convert temperature in Kelvin to InternalEnergy (u) in code units.

coolingRateToCGS_unused(coolrate)

Convert code units (du/dt) to erg/s/g (cgs, specific). Unused.

coolingRateToCGS(code_dens, code_gfmcoolrate)

Convert cooling/heating rate to specific CGS [erg/s/g].

Input is PartType0/[Masses,Density,GFM_CoolingRate].

powellEnergyTermCGS(code_dens, code_divb, code_b, code_vel, code_vol)

The ‘Powell heating/cooling’ energy source term (rightmost in Eqn. 21 Pakmor & Springel arxiv:1212.1452).

coolingTimeGyr(code_dens, code_gfmcoolrate, code_u)

Calculate a cooling time [Gyr] from gas cell properties.

Inputs: three snapshot values (i.e. code units) of Density, GFM_CoolingRate, InternalEnergy.

tracerEntToCGS(ent, log=False)

Fix cosmological/unit system in TRACER_MC[MaxEnt], output in cgs [K cm^2].

calcXrayLumBolometric(sfr, u, xe, mass, dens, temp=None, log=False)

Estimate bolometric X-ray luminosity of gas [10^30 erg/s], assuming only free-free (bremsstrahlung) emission.

Note: based only on gas density and temperature, and also assumes simplified (primordial) high-temp cooling function, and only free-free emission contribution from T>10^6 Kelvin gas. All inputs in code units.

opticalDepthLineCenter(transition, dens_cgs, temp_K, cellsize_code)

Derive an optical depth tau_0 = n * sigma_0 * L for a given transition.

Assume assuming the frequency is at line center (neglecting the Voigt profile shape). dens_cgs is the volume number density of the species of relevance [1/cm^3], temp is the cell temperature [linear K], and cellsize_code is the usual radius of the cell [code units], which we take as L/2.

sfrToHalphaLuminosity(sfr)

Convert SFR from code units (Msun/yr) into H-alpha line luminosity [linear 10^30 erg/s].

Just the usual linear conversion from Kennicutt.

gasSfrMetalMassToS850Flux(sfr, metalmass, temp, dens, ismCut=True)

Convert SFR [code units, Msun/yr] and metal mass [Msun] into 850 micron submm flux [linear mJy].

Simple model of Hayward+ (2013b) assuming a constant dust-to-metal ratio. Also enforce Torrey+12 ‘ISM’ constraint using temp [log K] and dens [code units, comoving].

calcEntropyCGS(u, dens, log=False)

Calculate entropy as P/rho^gamma, converting rho from comoving to physical. Return [K cm^2].

calcPressureCGS(u, dens, log=False)

Calculate pressure as (gamma-1)*u*rho in physical ‘cgs’ [K/cm^3] units.

calcMagneticPressureCGS(b, log=False)

Calculate magnetic pressure as B^2/8/pi in physical ‘cgs’ [K/cm^3] units.

calcKineticEnergyDensityCGS(dens_code, vel_kms, log=False)

Calculate kinetic energy density (KE/volume = 1/2 * mv^2 / volume) in ‘cgs’ [K/cm^3] units.

Inputs: dens_code in code units, vel_kms in physical km/s.

calcSoundSpeedKmS(u, dens, log=False)

Calculate sound speed as sqrt(gamma*Pressure/Density) in physical km/s.

soundSpeedFromTemp(temp, log=False)

Calculate sound speed given temperature [in Kelvin] in physical km/s.

calcSunyaevZeldovichYparam(mass, xe, temp)

Calculate per-cell (thermal) SZ y-parameter (McCarthy+2014 Eqn 2, Roncarelli+2007 Eqn 5, Kay+2012 Eqn 12).

Parameters:

  • mass (ndarray[float]) – gas cell masses [code units].

  • xe (ndarray[float]) – gas electron number density fraction [code units, i.e. dimensionless linear].

  • temp (ndarray[float]) – gas temperature [linear K].

Returns:

y-parameter in units of area [pkpc^2].

Return type:

ndarray[float]

calcKineticSZYParam(mass, xe, vel_los)

Calculate per-cell kinetic SZ y-parameter (e.g. Dolag+16 Eqn. 4, Altamura+23 Eqn. 2).

Parameters:

  • mass (ndarray[float]) – gas cell masses [code units].

  • xe (ndarray[float]) – gas electron number density fraction [code units, i.e. dimensionless linear].

  • vel_los (ndarray[float]) – line-of-sight velocity [km/s].

Returns:

y-parameter in units of area [pkpc^2].

Return type:

ndarray[float]

codeDensToCritRatio(rho, baryon=False, log=False, redshiftZero=False)

Normalize code density by the critical (total/baryonic) density at some redshift.

If redshiftZero, normalize by rho_crit,0 instead of rho_crit(z).

critRatioToCodeDens(ratioToCrit, baryon=False)

Convert a ratio of the critical density at some redshift to a code density.

codeMassToVirEnt(mass, log=False)

Given a total halo mass, return a S200 (e.g. Pvir/rho_200crit^gamma).

codeMassToVirVel(mass)

Given a total halo mass [in code units], return a virial velocity (V200) in physical [km/s].

codeMassToVirRad(mass)

Given a total halo mass [in code units], return a virial radius (r200) in physical [kpc].

codeM200R200ToV200InKmS(m200, r200)

Given a (M200,R200) pair in code units for a FoF group, compute V200 in physical [km/s].

avgEnclosedDensityToFreeFallTime(rho_code)

Convert a mass density (code units) to gravitational free-fall time [linear Gyr].

rho_code Should be the mean density interior to the radius R of a test mass within a halo.

metallicityInSolar(metal, log=False)

Given a code metallicity (M_Z/M_total), convert to value with respect to solar.

codeTimeStepToYears(TimeStep, Gyr=False)

Convert a TimeStep/TimeStepHydro/TimeStepGrav for a comoving run to a physical time in years.

Note: these inputs are an integer times All.Timebase_interval.

scalefacToAgeLogGyr(scalefacs)

Convert scalefactors of formation (e.g. GFM_StellarFormationTime) to age in log(Gyr).

Uses the current age of the universe as specified by sP.redshift.

codeEnergyToErg(energy, log=False)

Convert energy from code units (unitMass*unitLength^2/unitTime^2) to [erg]. (for BH_CumEgy*).

codeEnergyRateToErgPerSec(energy_rate, log=False)

Convert energy/time from code units (unitEnergy/unitTime) to [erg/s]. (for Gas EnergyDissipation).

codeEnergyDensToErgPerCm3(energy_dens, log=False)

Convert energy/volume from -proper- code units to [erg/cm^3] (for RadiationEnergyDensity).

codeEnergyDensToHabing(energy_dens, log=False)

Convert energy/volume from -proper- code units to Habing units (for RadiationEnergyDensity).

lumToAbsMag(lum)

Convert from an input luminosity in units of [Lsun/Hz] to an AB absolute magnitude.

absMagToLuminosity(mag)

Convert from input AB absolute magnitudes to (linear) luminosity units of [Lsun/Hz].

absMagToApparent(absolute_mag, redshift=None)

Convert an absolute magnitude to apparent.

apparentMagToAbsolute(apparent_mag, redshift=None)

Convert an apparent magnitude to absolute.

photonWavelengthToErg(wavelength, redshift=None)

Convert a photon wavelength [rest-frame Ang] emitted at a redshift into photon energy [erg].

luminosityToFlux(lum, wavelength=None, redshift=None)

Convert a luminosity in [erg/s] to a flux [photon/s/cm^2] for e.g. line emission.

At a given wavelength in [Angstroms] if not None, from a source at the given redshift. If wavelength is None, then output units are an energy flux e.g. [erg/s/cm^2].

fluxToLuminosity(flux, redshift=None)

Convert a flux in [erg/s/cm^2] to a luminosity [erg/s], from a source at the given redshift.

fluxToSurfaceBrightness(flux, pxDimsCode, arcsec2=True, arcmin2=False, ster=False, kpc=False)

Convert a flux in e.g. [energy/s/cm^2] or [photon/s/cm^2] into a surface brightness.

At a given redshift and for a certain pixel scale. pxDimsCode is a 2-tuple of the x and y dimensions of the pixel in code units, i.e. [ckpc/h]^2. Output e.g.: [photon/s/cm^2/arcsec^2] if arcsec2 == True, otherwise possibly [flux/arcmin^2], [flux/ster], or [flux/kpc^2].

synchrotronPowerPerFreq(gas_B, gas_vol, watts_per_hz=True, log=False, telescope='SKA', eta=1.0, k=10, gamma_min=300, gamma_max=15000, alpha=1.7)

Calculate synchrotron power per unit frequency (simple model) for gas cells.

Inputs: gas_B and gas_vol the magnetic field 3-vector and volume, both in code units. Output units [Watts/Hz]. Default model parameter assumptions of Xu+ (2012)/Marinacci+ (2017).

redshiftToAgeFlat(z)

Calculate age of the universe [Gyr] at the given redshift (assuming flat cosmology).

Analytical formula from Peebles, p.317, eq 13.2.

ageFlatToRedshift(age)

Calculate redshift from age of the universe [Gyr] (assuming flat cosmology).

Inversion of analytical formula from redshiftToAgeFlat().

redshiftToLookbackTime(z)

Calculate lookback time from z=0 to redshift in [Gyr], assuming flat cosmology.

redshiftToComovingDist(z)

Convert redshift z to line of sight distance (in Mpc). Assumes flat.

redshiftToComovingVolume(z)

Calculate total comoving volume from the present to redshift z (all-sky) in [cGpc^3].

comovingVolumeDeltaRedshift(z1, z2, arcSeqSq=None, arcMinSq=None, arcDegSq=None, mpc3=False)

Calcuate comoving volume between z1 and z2 (all-sky).

If any of arcSeqSq, arcMinSq, or arcDegSq are input, then return for this solid angle on the sky instead of 4pi. If mpc3, return in units of [cMpc^3] instead of the default [cGpc^3].

redshiftToAngDiamDist(z)

Convert redshift z to angular diameter distance (in Mpc).

This equals the proper/physical transverse distance for theta=1 rad. Assumes flat. Peebles, p.325.

redshiftToLumDist(z)

Convert redshift z to luminosity distance (in Mpc).

This then allows the conversion between luminosity and a flux at that redshift.

arcsecToAngSizeKpcAtRedshift(ang_diam, z=None)

Convert an angle in arcseconds to an angular/transverse size (in proper/physical kpc) at redshift z.

Assumes flat cosmology.

arcsecToCodeLength(x_arcsec, z=None)

Convert an angle in arcseconds to an angular/tranverse size (in code length units).

degToAngSizeKpcAtRedshift(ang_diam_deg, z)

Convert an angle in degrees to a physical size [kpc] at redshift z.

codeLengthToAngularSize(length_codeunits, z=None, arcsec=True, arcmin=False, deg=False)

Convert a distance in code units (i.e. ckpc/h) to an angular scale in [arcsec] at a given redshift z.

Assumes flat cosmology. If arcmin or deg is True, then [arcmin] or [deg].

physicalKpcToAngularSize(length_kpc, z=None, arcsec=True, arcmin=False, deg=False)

Convert a distance in physical kpc to an angular scale in [arcsec] at the given redshift z.

Assumes flat cosmology. If arcmin or deg is True, then [arcmin] or [deg].

magsToSurfaceBrightness(vals, pxSizeCode)

Convert magnitudes (probably on a grid) to surface brightness values, given a constant pxSizeCode.

haloMassToOtherOverdensity(mass_orig, delta_orig=200, delta_new=500)

Convert a halo mass between two different spherical overdensity definitions, e.g. M200 to M500.

Assumes NFW density profile and the c-M relation from Bullock by default. Based on https://arxiv.org/abs/astro-ph/0203169 (appendix). Input halo mass (and output) in units of linear Msun.

m200_to_m500(halo_mass)

Convert a halo mass from Delta=200 to Delta=500.

m500_to_m200(halo_mass)

Convert a halo mass from Delta=500 to Delta=200.

particleCountToMass(N_part, baryon=True, boxLength=None)

Convert the cube-root of a total particle count (e.g. 512, 1820) to its average mass.

Uses the cosmology and volume of the box. If boxLength is specified, then should be a box side-length in cMpc/h units (e.g. 25, 205). Return in [Msun].

meanmolwt(Y, Z)

Mean molecular weight, from Monaco+ (2007) eqn 14, for hot halo gas.

Y = helium fraction (0.25) Z = metallicity (non-log metal mass/total mass)

temet.util.virtualSimFile

Creation and update of the virtual ‘simulation.hdf5’ file.

createVirtualSimHDF5()

Create a single ‘simulation.hdf5’ file which is made up of virtual datasets (HDF5 1.1x/h5py 2.9.x features).

Note: dataset details acquired from first chunk of last snapshot! Snapshot must be full, and first chunk must have at least one of every particle type! Note: run in simulation root dir, since we make relative path links.

supplementVirtualSimHDF5AddSnapField()

Add to existing ‘simulation.hdf5’ file (modify as needed, careful!).

supplementVirtualSimHDF5AddOrUpdateGroupcatField()

Add to existing ‘simulation.hdf5’ file (modify as needed, careful!).

supplementVirtualSimHDF5()

Add to existing ‘simulation.hdf5’ file (modify as needed, careful!).

temet.util.voronoi

Algorithms and methods related to the use and analysis of the Voronoi mesh.

loadSingleHaloVPPP(sP, haloID)

Load Voronoi connectivity information for a single FoF halo.

loadGlobalVPPP(sP)

Load global Voronoi connectivity information for a snapshot.

voronoiThresholdSegmentation(sP, haloID, propName, propThresh, propThreshComp)

Apply a Voronoi thresholded segmentation algorithm to all gas cells in a halo.

Identifies spatially connected collectins, in the sense that they are Voronoi natural neighbors, and which satisfy a threshold criterion on a particular gas property, e.g. log(T) < 4.5.

test_cen()

Test circumsphere center and coplanarity functions.

plotFacePolygon(cell_cen, neighbor_cen, vertices)

Plot helper.

voronoiSliceWithPlane()

Testing.

temet.util.voronoiRay

Algorithms and methods related to ray-tracing through a Voronoi mesh.

trace_ray_through_voronoi_mesh_with_connectivity(cell_pos, num_ngb, ngb_inds, offset_ngb, ray_pos_in, ray_dir, total_dl, boxSize, debug, verify, fof_scope_mesh)

Ray-trace a single ray through a Voronoi mesh, with pre-computed neighbor connectivity information.

The ray is specified by its starting location, direction, and length.

Parameters:

  • cell_pos (ndarray[float]) – Voronoi cell center positions [N,3].

  • num_ngb (array[int]) – Voronoi mesh connectivity, first return of util.voronoi.loadSingleHaloVPPP() or util.voronoi.loadGlobalVPPP() .

  • ngb_inds (array[int]) – Voronoi mesh connectivity, second return of util.voronoi.loadSingleHaloVPPP() or util.voronoi.loadGlobalVPPP() .

  • offset_ngb (array[int]) – Voronoi mesh connectivity, third return of util.voronoi.loadSingleHaloVPPP() or util.voronoi.loadGlobalVPPP() .

  • ray_pos_in (float[3]) – x,y,z coordinates of ray starting position.

  • ray_dir (float[3]) – normalized unit vector specifying ray direction.

  • total_dl (float) – pathlength to integrate ray (code units i.e. same as cell_pos).

  • boxSize (float) – simulation box size, for periodic boundaries (code units).

  • debug (bool) – if True, >=1, or >=2, print increasing verbose debugging information (disabled for @jit).

  • verify (bool) – if True, do brute-force distance verification every step of parent Voronoi cell.

  • fof_scope_mesh (bool) – if True, indicate that we have loaded and are working with a fof-scope set of cell data and mesh connectivity, i.e. not a correct not periodic mesh at the edges.

Returns:

a 2-tuple composed of

  • dx (ndarray[float]): per-cell path length, in order, for each intersected cell (code/input units).

  • ind (ndarray[float]): per-cell index, in order, for each intersected cell.

trace_ray_through_voronoi_mesh_treebased(cell_pos, NextNode, length, center, sibling, nextnode, ray_pos_in, ray_dir, total_dl, boxSize, debug, verify)

Ray-trace a single ray through a Voronoi mesh, using searches in a pre-computed oct-tree.

Parameters:

  • cell_pos (ndarray[float]) – Voronoi cell center positions [N,3].

  • NextNode (array[int]) – neighbor tree, first return of util.treeSearch.buildFullTree().

  • length (array[int]) – neighbor tree, second return of util.treeSearch.buildFullTree().

  • center (array[int]) – neighbor tree, third return of util.treeSearch.buildFullTree().

  • sibling (array[int]) – neighbor tree, fourth return of util.treeSearch.buildFullTree().

  • nextnode (array[int]) – neighbor tree, fifth return of util.treeSearch.buildFullTree().

  • ray_pos_in (float[3]) – x,y,z coordinates of ray starting position.

  • ray_dir (float[3]) – normalized unit vector specifying ray direction.

  • total_dl (float) – pathlength to integrate ray (code units i.e. same as cell_pos).

  • boxSize (float) – simulation box size, for periodic boundaries (code units).

  • debug (bool) – if True, >=1, or >=2, print increasing verbose debugging information (disabled for @jit).

  • verify (bool) – if True, do brute-force distance verification every step of parent Voronoi cell.

Returns:

a 2-tuple composed of

  • dx (ndarray[float]): per-cell path length, in order, for each intersected cell (code/input units).

  • ind (ndarray[float]): per-cell index, in order, for each intersected cell.

rayTrace(sP, ray_pos, ray_dir, total_dl, pos, quant=None, quant2=None, mode='full', nThreads=72, tree=None)

Tree-based ray-tracing through a Voronoi mesh.

For a given set of particle coordinates, assumed to be Voronoi cell center positions, perform a tree-based ray-tracing through this implicit Voronoi mesh, for a number of specified ray_pos initial ray positions, in a given direction, for a total pathlength. The operation while tracing is specified by mode.

Parameters:

  • sP (util.simParams.simParams) – simParams instance.

  • ray_pos (ndarray[float][M,3]) – starting ray positions in (x,y,z) space.

  • ray_dir (float[3]) – normalized unit vector specifying (the same) direction for all rays.

  • total_dl (float) – path length to integrate all rays (same units as pos).

  • pos (ndarray[float][N,3]) – array of 3-coordinates for the gas cells.

  • quant (ndarray[float][N]) – if not None, a quantity with the same size as pos to operate on.

  • quant2 (ndarray[float][N]) – if not None, a second quantity with the same size as pos, for e.g. weighting.

  • mode (str) – one of the following modes of operation: full: return full rays as a 4-tuple of (offsets, lengths, cell_dx, cell_inds) count, dx_sum, etc: return a summary statistic for each ray (single array of values)

  • nThreads (int) – do multithreaded calculation (mem required=nThreads times more).

  • tree (list or None) – if not None, should be a list of all the needed tree arrays (pre-computed), i.e the exact return of util.treeSearch.buildFullTree().