Source code for pynbody.analysis.halo

"""

halo
====

Functions for dealing with and manipulating halos in simulations.


"""

from .. import filt, util, config, array, units, transformation
from . import cosmology, _com, profile
import numpy as np
import math
import logging
logger = logging.getLogger('pynbody.analysis.halo')


[docs]def center_of_mass(sim): """ Return the centre of mass of the SimSnap """ mtot = sim["mass"].sum() p = np.sum(sim["mass"] * sim["pos"].transpose(), axis=1) / mtot # otherwise behaviour is numpy version dependent p.units = sim["pos"].units # only return position to be consistent with other functions in halo.py return p
[docs]def center_of_mass_velocity(sim): """ Return the center of mass velocity of the SimSnap """ mtot = sim["mass"].sum() v = np.sum(sim["mass"] * sim["vel"].transpose(), axis=1) / mtot # otherwise behaviour is numpy version dependent v.units = sim["vel"].units return v
[docs]def shrink_sphere_center(sim, r=None, shrink_factor=0.7, min_particles=100, verbose=False, num_threads = config['number_of_threads'],**kwargs): """ Return the center according to the shrinking-sphere method of Power et al (2003) **Input**: *sim* : a simulation snapshot - this can be any subclass of SimSnap **Optional Keywords**: *r* (default=None): initial search radius. This can be a string indicating the unit, i.e. "200 kpc", or an instance of :func:`~pynbody.units.Unit`. *shrink_factor* (default=0.7): the amount to shrink the search radius by on each iteration *min_particles* (default=100): minimum number of particles within the search radius. When this number is reached, the search is complete. *verbose* (default=False): if True, prints out the diagnostics at each iteration. Useful to determine whether the centering is zeroing in on the wrong part of the simulation. """ if r is None: # use rough estimate for a maximum radius # results will be insensitive to the exact value chosen r = (sim["x"].max() - sim["x"].min()) / 2 elif isinstance(r, str) or issubclass(r.__class__, units.UnitBase): if isinstance(r, str): r = units.Unit(r) r = r.in_units(sim['pos'].units, **sim.conversion_context()) mass = np.asarray(sim['mass'], dtype='double') pos = np.asarray(sim['pos'], dtype='double') if shrink_factor == 1.0: tol = sim['eps'].in_units(sim['pos'].units, **sim.conversion_context()).min()*0.1 com = _com.shrink_sphere_center(pos, mass, min_particles, shrink_factor, r, num_threads) com = _com.move_sphere_center(pos, mass, min_particles, shrink_factor, r, tol) else: com = _com.shrink_sphere_center(pos, mass, min_particles, shrink_factor, r, num_threads) logger.info("Final SSC=%s", com) return array.SimArray(com, sim['pos'].units)
[docs]def virial_radius(sim, cen=None, overden=178, r_max=None, rho_def='matter'): """Calculate the virial radius of the halo centered on the given coordinates. The default is here defined by the sphere centered on cen which contains a mean density of overden * rho_M_0 * (1+z)^3. **Input**: *sim* : a simulation snapshot - this can be any subclass of SimSnap, especially a halo. **Optional Keywords**: *cen* (default=None): Provides the halo center. If None, assumes that the snapshot is already centered. *rmax (default=None): Maximum radius to start the search. If None, inferred from the halo particle positions. *overden (default=178): Overdensity corresponding to the required halo boundary definition. 178 is the virial criterion for spherical collapse in an EdS Universe. Common possible values are 200, 500 etc *rho_def (default='matter'): Physical density used to define the overdensity. Default is the matter density at the redshift of the simulation. An other choice is "critical" for the critical density at this redshift. """ if r_max is None: r_max = (sim["x"].max() - sim["x"].min()) else: if cen is not None: sim = sim[filt.Sphere(r_max, cen)] else: sim = sim[filt.Sphere(r_max)] r_min = 0.0 if cen is not None: tx = transformation.inverse_translate(sim, cen) else: tx = transformation.null(sim) if rho_def == 'matter': ref_density = sim.properties["omegaM0"] * cosmology.rho_crit(sim, z=0) * (1.0 + sim.properties["z"]) ** 3 elif rho_def == 'critical': ref_density = cosmology.rho_crit(sim, z=sim.properties["z"]) else: raise ValueError(rho_def + "is not a valid definition for the reference density") target_rho = overden * ref_density logger.info("target_rho=%s", target_rho) with tx: sim = sim[filt.Sphere(r_max)] with sim.immediate_mode: mass_ar = np.asarray(sim['mass']) r_ar = np.asarray(sim['r']) """ #pure numpy implementation rho = lambda r: np.dot( mass_ar, r_ar < r) / (4. * math.pi * (r ** 3) / 3) #numexpr alternative - not much faster because sum is not threaded def rho(r) : r_ar; mass_ar; # just to get these into the local namespace return ne.evaluate("sum((r_ar<r)*mass_ar)")/(4.*math.pi*(r**3)/3) """ rho = lambda r: util.sum_if_lt(mass_ar,r_ar,r)/(4. * math.pi * (r ** 3) / 3) result = util.bisect(r_min, r_max, lambda r: target_rho - rho(r), epsilon=0, eta=1.e-3 * target_rho, verbose=False) return result
def potential_minimum(sim): i = sim["phi"].argmin() return sim["pos"][i].copy()
[docs]def hybrid_center(sim, r='3 kpc', **kwargs): """ Determine the center of the halo by finding the shrink-sphere -center inside the specified distance of the potential minimum """ try: cen_a = potential_minimum(sim) except KeyError: cen_a = center_of_mass(sim) return shrink_sphere_center(sim[filt.Sphere(r, cen_a)], **kwargs)
[docs]def index_center(sim, **kwargs): """ Determine the center of mass based on specific particles. Supply a list of indices using the ``ind`` keyword. """ try: ind = kwargs['ind'] return center_of_mass(sim[ind]) except KeyError: raise RuntimeError("Need to supply indices for centering")
[docs]def vel_center(sim, mode=None, cen_size="1 kpc", retcen=False, move_all=True, **kwargs): """Use stars from a sphere to calculate center of velocity. The size of the sphere is given by the ``cen_size`` keyword and defaults to 1 kpc. **Keyword arguments:** *mode*: reserved for future use; currently ignored *move_all*: if True (default), move the entire snapshot. Otherwise only move the particles in the halo passed in. *retcen*: if True only return the velocity center without moving the snapshot (default = False) """ logger.info("Finding halo velocity center...") if move_all: target = sim.ancestor else: target = sim cen = sim.star[filt.Sphere(cen_size)] if len(cen) < 5: # fall-back to DM cen = sim.dm[filt.Sphere(cen_size)] if len(cen) < 5: # fall-back to gas cen = sim.gas[filt.Sphere(cen_size)] if len(cen) < 5: # very weird snapshot, or mis-centering! raise ValueError("Insufficient particles around center to get velocity") vcen = (cen['vel'].transpose() * cen['mass']).sum(axis=1) / \ cen['mass'].sum() vcen.units = cen['vel'].units if config['verbose']: logger.info("vcen=%s", vcen) if retcen: return vcen else: return transformation.v_translate(target, -vcen)
[docs]def center(sim, mode=None, retcen=False, vel=True, cen_size="1 kpc", move_all=True, wrap=False, **kwargs): """ Determine the center of mass of the given particles using the specified mode, then recenter the particles (of the entire ancestor snapshot) accordingly Accepted values for *mode* are *pot*: potential minimum *com*: center of mass *ssc*: shrink sphere center *ind*: center on specific particles; supply the list of particles using the ``ind`` keyword. *hyb*: for sane halos, returns the same as ssc, but works faster by starting iteration near potential minimum or a function returning the COM. **Other keywords:** *retcen*: if True only return the center without centering the snapshot (default = False) *ind*: only used when *mode=ind* -- specifies the indices of particles to be used for centering *vel*: if True, translate velocities so that the velocity of the central 1kpc (default) is zeroed. Other values can be passed with cen_size. *move_all*: if True (default), move the entire snapshot. Otherwise only move the particles in the halo passed in. *wrap*: if True, pre-centre and wrap the simulation so that halos on the edge of the box are handled correctly. Default False. """ global config if mode is None: mode = config['centering-scheme'] try: fn = {'pot': potential_minimum, 'com': center_of_mass, 'ssc': shrink_sphere_center, 'hyb': hybrid_center, 'ind': index_center}[mode] except KeyError: fn = mode if move_all: target = sim.ancestor else: target = sim if wrap: # centre on something within the halo and wrap target = transformation.inverse_translate(target, sim['pos'][0]) target.sim.wrap() if retcen: return fn(sim, **kwargs) else: cen = fn(sim, **kwargs) tx = transformation.inverse_translate(target, cen) if vel: velc = vel_center(sim, cen_size=cen_size, retcen=True) tx = transformation.inverse_v_translate(tx, velc) return tx
[docs]def halo_shape(sim, N=100, rin=None, rout=None, bins='equal'): """ Returns radii in units of ``sim['pos']``, axis ratios b/a and c/a, the alignment angle of axis a in radians, and the rotation matrix for homeoidal shells over a range of N halo radii. **Keyword arguments:** *N* (100): The number of homeoidal shells to consider. Shells with few particles will take longer to fit. *rin* (None): The minimum radial bin in units of ``sim['pos']``. Note that this applies to axis a, so particles within this radius may still be included within homeoidal shells. By default this is taken as rout/1000. *rout* (None): The maximum radial bin in units of ``sim['pos']``. By default this is taken as the largest radial value in the halo particle distribution. *bins* (equal): The spacing scheme for the homeoidal shell bins. ``equal`` initialises radial bins with equal numbers of particles, with the exception of the final bin which will accomodate remainders. This number is not necessarily maintained during fitting. ``log`` and ``lin`` initialise bins with logarithmic and linear radial spacing. Halo must be in a centered frame. Caution is advised when assigning large number of bins and radial ranges with many particles, as the algorithm becomes very slow. """ #-----------------------------FUNCTIONS----------------------------- # Define an ellipsoid shell with lengths a,b,c and orientation E: def Ellipsoid(r, a,b,c, E): x,y,z = np.dot(np.transpose(E),[r[:,0],r[:,1],r[:,2]]) return (x/a)**2 + (y/b)**2 + (z/c)**2 # Define moment of inertia tensor: MoI = lambda r,m: np.array([[np.sum(m*r[:,i]*r[:,j]) for j in range(3)]\ for i in range(3)]) # Splits 'r' array into N groups containing equal numbers of particles. # An array is returned with the radial bins that contain these groups. sn = lambda r,N: np.append([r[i*int(len(r)/N):(1+i)*int(len(r)/N)][0]\ for i in range(N)],r[-1]) # Retrieves alignment angle: almnt = lambda E: np.arccos(np.dot(np.dot(E,[1.,0.,0.]),[1.,0.,0.])) #-----------------------------FUNCTIONS----------------------------- if (rout == None): rout = sim.dm['r'].max() if (rin == None): rin = rout/1E3 posr = np.array(sim.dm['r'])[np.where(sim.dm['r'] < rout)] pos = np.array(sim.dm['pos'])[np.where(sim.dm['r'] < rout)] mass = np.array(sim.dm['mass'])[np.where(sim.dm['r'] < rout)] rx = [[1.,0.,0.],[0.,0.,-1.],[0.,1.,0.]] ry = [[0.,0.,1.],[0.,1.,0.],[-1.,0.,0.]] rz = [[0.,-1.,0.],[1.,0.,0.],[0.,0.,1.]] # Define bins: if (bins == 'equal'): # Each bin contains equal number of particles mid = sn(np.sort(posr[np.where((posr >= rin) & (posr <= rout))]),N*2) rbin = mid[1:N*2+1:2] mid = mid[0:N*2+1:2] elif (bins == 'log'): # Bins are logarithmically spaced mid = profile.Profile(sim.dm, type='log', ndim=3, rmin=rin, rmax=rout, nbins=N+1)['rbins'] rbin = np.sqrt(mid[0:N]*mid[1:N+1]) elif (bins == 'lin'): # Bins are linearly spaced mid = profile.Profile(sim.dm, type='lin', ndim=3, rmin=rin, rmax=rout, nbins=N+1)['rbins'] rbin = 0.5*(mid[0:N]+mid[1:N+1]) # Define b/a and c/a ratios and angle arrays: ba,ca,angle = np.zeros(N),np.zeros(N),np.zeros(N) Es = [0]*N # Begin loop through radii: for i in range(0,N): # Initialise convergence criterion: tol = 1E-3 count = 0 # Define initial spherical shell: a=b=c = rbin[i] E = np.identity(3) L1,L2 = rbin[i]-mid[i],mid[i+1]-rbin[i] # Begin iterative procedure to fit data to shell: while True: count += 1 # Collect all particle positions and masses within shell: r = pos[np.where((posr < a+L2) & (posr > c-L1*c/a))] inner = Ellipsoid(r, a-L1,b-L1*b/a,c-L1*c/a, E) outer = Ellipsoid(r, a+L2,b+L2*b/a,c+L2*c/a, E) r = r[np.where((inner > 1.) & (outer < 1.))] m = mass[np.where((inner > 1.) & (outer < 1.))] # End iterations if there is no data in range: if (len(r) == 0): ba[i],ca[i],angle[i],Es[i] = b/a,c/a,almnt(E),E logger.info('No data in range after %i iterations' %count) break # Calculate shape tensor & diagonalise: D = list(np.linalg.eig(MoI(r,m)/np.sum(m))) # Purge complex numbers: if isinstance(D[1][0,0],complex): D[0] = D[0].real ; D[1] = D[1].real logger.info('Complex numbers in D removed...') # Compute ratios a,b,c from moment of intertia principles: anew,bnew,cnew = np.sqrt(abs(D[0])*3.0) # The rotation matrix must be reoriented: E = D[1] if ((bnew > anew) & (anew >= cnew)): E = np.dot(E,rz) if ((cnew > anew) & (anew >= bnew)): E = np.dot(np.dot(E,ry),rx) if ((bnew > cnew) & (cnew >= anew)): E = np.dot(np.dot(E,rz),rx) if ((anew > cnew) & (cnew >= bnew)): E = np.dot(E,rx) if ((cnew > bnew) & (bnew >= anew)): E = np.dot(E,ry) cnew,bnew,anew = np.sort(np.sqrt(abs(D[0])*3.0)) # Keep a as semi-major axis and distort b,c by b/a and c/a: div = rbin[i]/anew anew *= div bnew *= div cnew *= div # Convergence criterion: if (np.abs(b/a-bnew/anew) < tol) & (np.abs(c/a-cnew/anew) < tol): if (almnt(-E) < almnt(E)): E = -E ba[i],ca[i],angle[i],Es[i] = bnew/anew,cnew/anew,almnt(E),E break # Increase tolerance if convergence has stagnated: elif (count%10 == 0): tol *= 5. # Reset a,b,c for the next iteration: a,b,c = anew,bnew,cnew return [array.SimArray(rbin, sim.d['pos'].units), ba, ca, angle, Es]