Secure your code as it's written. Use Snyk Code to scan source code in minutes - no build needed - and fix issues immediately.
def SolarNeighborhood(r1=units.Unit("5 kpc"), r2=units.Unit("10 kpc"), height=units.Unit("2 kpc"), cen=(0, 0, 0)):
"""
Convenience function that returns a filter which selects particles
in a disc between radii `r1` and `r2` and thickness `height`.
"""
x = Disc(max(r1, r2), height, cen) & ~Disc(min(r1, r2), height, cen)
x._descriptor = "Solar Neighborhood"
return x
"""
if center:
angmom.faceon(sim)
if min:
min_r = min
else:
min_r = sim['rxy'].min()
if max:
max_r = max
else:
max_r = sim['rxy'].max()
if isinstance(pretime, str):
pretime = units.Unit(pretime)
diskstars = sim.star[filt.Disc(max_r, disk_height)]
youngstars = np.where(diskstars['tform'].in_units("Myr") >
sim.properties['time'].in_units(
"Myr", **sim.conversion_context())
- pretime.in_units('Myr'))[0]
pro = profile.Profile(diskstars[youngstars], type=bin_spacing,
nbins=nbins, min=min_r, max=max_r)
r = pro['rbins'].in_units(r_units)
fourierprof = pro['fourier']
a2 = fourierprof['amp'][2]
if clear:
p.clf()
width = _units.Unit(width)
width = width.in_units(sim['pos'].units, **sim.conversion_context())
width = float(width)
pixel_size = width / vector_resolution
X, Y = np.meshgrid(np.arange(-width / 2 + pixel_size/2, width / 2 + pixel_size/2, pixel_size ),
np.arange(-width / 2 + pixel_size/2, width / 2 + pixel_size/2, pixel_size))
im = image(sim, width=width, **kwargs)
if mode == 'quiver':
if scale is None:
Q = p.quiver(X, Y, vx, vy, color=vector_color, edgecolor=edgecolor)
else:
Q = p.quiver(X, Y, vx, vy, scale=_units.Unit(scale).in_units(
sim['vel'].units), color=vector_color, edgecolor=edgecolor)
if quiverkey:
qk = p.quiverkey(Q, key_x, key_y, key_unit.in_units(sim['vel'].units, **sim.conversion_context()),
r"$\mathbf{" + key_unit.latex() + "}$", labelcolor=key_color, color=key_color, fontproperties={'size': 24})
if quiverkey_bg_color is not None:
qk.text.set_backgroundcolor(quiverkey_bg_color)
elif mode == 'stream' :
Q = p.streamplot(X, Y, vx, vy, color=vector_color, density=density)
# RS - if a axis object is passed in, use the right limit call
if subplot:
p.set_xlim(-width/2, width/2)
p.set_ylim(-width/2, width/2)
else:
p.xlim(-width/2, width/2)
p.ylim(-width/2, width/2)
*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.
"""
import os
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')
com = _com.shrink_sphere_center(pos, mass, min_particles, shrink_factor, r)
logger.info("Final SSC=%s", com)
return array.SimArray(com, sim['pos'].units)
simulation is in comoving units and you specify a different
redshift z. Specifically, the physical density for the redshift
you specify is calulated, but expressed as a comoving density *at
the redshift of the snapshot*. This is intentional behaviour."""
if z is None:
z = f.properties['z']
if unit is None:
try:
unit = f.dm["mass"].units / f.dm["pos"].units ** 3
except units.UnitsException:
unit = units.NoUnit()
if hasattr(unit, "_no_unit"):
unit = units.Unit("Msol kpc^-3 a^-3")
omM = f.properties['omegaM0']
omL = f.properties['omegaL0']
h0 = f.properties['h']
a = 1.0 / (1.0 + z)
H_z = _a_dot(a, h0, omM, omL) / a
H_z = units.Unit("100 km s^-1 Mpc^-1") * H_z
rho_crit = (3 * H_z ** 2) / (8 * math.pi * units.G)
return rho_crit.ratio(unit, **f.conversion_context())
def _get_units_from_description(self, description, expectedCgsConversionFactor=None):
arr_units = units.Unit('1.0')
conversion = 1.0
for unitname in self._hdf_unitvar.keys():
power = 1.
if unitname in description:
sstart = description.find(unitname)
if sstart > 0:
if description[sstart - 1] == "/":
power *= -1.
if len(description) > sstart + len(unitname):
# Just check we're not at the end of the line
if description[sstart + len(unitname)] == '^':
## Has an index, check if this is negative
if description[sstart + len(unitname) + 1] == "-":
power *= -1.
power *= float(
description[sstart + len(unitname) + 2:-1].split()[0]) ## Search for the power
self.max = units.Unit(kwargs['max']).ratio(x.units,
**sim.conversion_context())
else:
self.max = kwargs['max']
else:
self.max = np.max(x)
if kwargs.has_key('bins'):
self.nbins = len(kwargs['bins']) - 1
elif kwargs.has_key('nbins'):
self.nbins = kwargs['nbins']
else:
self.nbins = 100
if kwargs.has_key('min'):
if isinstance(kwargs['min'], str):
self.min = units.Unit(kwargs['min']).ratio(x.units,
**sim.conversion_context())
else:
self.min = kwargs['min']
else:
self.min = np.min(x[x > 0])
if kwargs.has_key('bins'):
self._properties['bin_edges'] = kwargs['bins']
self.min = kwargs['bins'].min()
self.max = kwargs['bins'].max()
elif type == 'log':
self._properties['bin_edges'] = np.logspace(
np.log10(self.min), np.log10(self.max), num=self.nbins + 1)
elif type == 'lin':
self._properties['bin_edges'] = np.linspace(
self.min, self.max, num=self.nbins + 1)
gv = gpro['rotation_curve_spherical'].in_units(v_units)
dv = dpro['rotation_curve_spherical'].in_units(v_units)
sv = spro['rotation_curve_spherical'].in_units(v_units)
else:
gv = gpro['v_circ'].in_units(v_units)
dv = dpro['v_circ'].in_units(v_units)
sv = spro['v_circ'].in_units(v_units)
p.plot(r, gv, "--", label="gas")
p.plot(r, dv, label="dark")
p.plot(r, sv, linestyle="dotted", label="star")
else:
p.plot(r, v, **kwargs)
if yrange:
p.axis(
[min_r, units.Unit(max_r).in_units(r.units), yrange[0], yrange[1]])
if not axes:
p.xlabel("r / $" + r.units.latex() + "$", fontsize='large')
p.ylabel("v$_c / " + v.units.latex() + '$', fontsize='large')
if legend:
p.legend(loc=0)
if (filename):
logger.info("Saving %s", filename)
p.savefig(filename)
return r, v
@Profile.profile_property
def sb(self,band='v') :
# At 10 pc (distance for absolute magnitudes), 1 arcsec is 10 AU=1/2.06e4 pc
# In [5]: (np.tan(np.pi/180/3600)*10.0)**2
# Out[5]: 2.3504430539466191e-09
# 1 square arcsecond is thus 2.35e-9 pc^2
sqarcsec_in_bin = self._binsize.in_units('pc^2') / 2.3504430539466191e-09
bin_luminosity = 10.0**(-0.4*self['magnitudes,'+band])
#import pdb; pdb.set_trace()
surfb = -2.5*np.log10(bin_luminosity / sqarcsec_in_bin)
surfb = array.SimArray(surfb, units.Unit('1'))
surfb.sim = self.sim
return surfb
it is used in place of the redshift in output f."""
if z is None:
z = f.properties['z']
a = 1. / (1. + z)
omegam0 = f.properties['omegaM0']
omegal0 = f.properties['omegaL0']
b, X = _lingrowthfac(z, omegam0, omegal0, return_norm=True)
I = _lingrowthintegrand(a, omegam0)
term1 = -(1.5 * omegam0 * a ** -3) * b / \
math.sqrt(1. - omegam0 + omegam0 * a ** -3)
term2 = (2.5 * omegam0) * hzoverh0(a, omegam0) ** 2 * a * I / X
res = units.h * (term1 + term2) * 100. * units.Unit("km s^-1 Mpc^-1")
return res.in_units(unit, **f.conversion_context())