17.01 Solution to the Illustris problems
Here we consider a halo at the three redshifts z=0, z=0.5 and z=1 and study the mass distribution.
#required packages
import numpy as np
import matplotlib.pyplot as plt
import requests
import json
import h5py
%matplotlib inline
These two functions facilitate the communication with the illustris server:
(they are provided on the Illustris project web page)
def get(path, params=None):
# make HTTP GET request to path
headers = {"api-key": "2566d8dd2bcf9aefbb3d8b01080a877b"}
r = requests.get(path, params=params, headers=headers)
# raise exception if response code is not HTTP SUCCESS (200)
r.raise_for_status()
if r.headers['content-type'] == 'application/json':
return r.json() # parse json responses automatically
if 'content-disposition' in r.headers:
filename = r.headers['content-disposition'].split("filename=")[1]
with open(filename, 'wb') as f:
f.write(r.content)
return filename # return the filename string
return r
def getsnapshot(id, redshift):
params = {'dm': 'Coordinates'}
url = 'http://www.illustris-project.org/api/Illustris-1/snapshots/z=' + str(redshift) + '/subhalos/' + str(id)
subhalo = get(url) # read the parameters of the subhalo in the json format
return subhalo, get(url+"/cutout.hdf5", params=params) # read the cutout of the subhalo into memory
1. Determine masses
redshifts = [0.0, 0.5, 1.0]
id = 75 # 75, 1030, 543681
dm = 6.3e6
for z in redshifts:
sub, filename = getsnapshot(id, z) # get snapshot with ID=id and given redshift
with h5py.File(filename) as f:
dx = f['PartType1']['Coordinates'][:,0] - sub['pos_x']
dy = f['PartType1']['Coordinates'][:,1] - sub['pos_y']
dz = f['PartType1']['Coordinates'][:,2] - sub['pos_z']
rr = np.sqrt(dx**2+dy**2+dz**2)
print('mass inside 50 kpc at z={0:.1f} is M={1:5.2f} x 10^10 Msun'.format(z, len(rr[rr<50])*dm/1e10))
Aside: Plot particle distributions
colors = ['red', 'green', 'blue']
redshifts = [0.0, 0.5, 1.0]
id = 75 # 75, 1030, 543681
plt.figure(figsize=(15,5))
col = 1
for z in redshifts:
sub, filename = getsnapshot(id, z)
ax = plt.subplot(1, 3, col)
with h5py.File(filename) as f:
dx = f['PartType1']['Coordinates'][:,0] - sub['pos_x']
dy = f['PartType1']['Coordinates'][:,1] - sub['pos_y']
dz = f['PartType1']['Coordinates'][:,2] - sub['pos_z']
print('x=', sub['pos_x'], 'y=', sub['pos_y'], 'z=', sub['pos_z'])
ax.scatter(dx, dy, c=colors[col-1], alpha=0.3)
ax.legend(['z='+str(z)])
ax.set_xlim(-200, 200)
ax.set_ylim(-200, 200)
col += 1
2. Determine the density profile $\frac{\mathrm{d}N}{\mathrm{d}V}$.
The density profile (objects per spherical shell) can be determined by binning the separation array [r1,r2,r3,...] for all dark matter particles with plt.hist.
nbins = 20
maxradius = 200
xp = 0*np.ones((3,nbins))
xpsq = 0*xp
e = 0*xp
count= 0*xp
redshifts = [0.0, 0.5, 1.0]
id = 1030 # 75, 1030, 543681
i = 0
for z in redshifts:
sub, filename = getsnapshot(id, z)
with h5py.File(filename) as f:
dx = f['PartType1']['Coordinates'][:,0] - sub['pos_x']
dy = f['PartType1']['Coordinates'][:,1] - sub['pos_y']
dz = f['PartType1']['Coordinates'][:,2] - sub['pos_z']
rr = np.sqrt(dx**2 + dy**2 + dz**2) # radii of all dm particles
num, x, ax = plt.hist(rr, bins=nbins, range=(0,maxradius), alpha=0.3)
# num: numbers in each bin
# x: position
xp[i,:] = (x[1:] + x[:-1])/2. # the middle of each bin
xpsq[i,:] = (4.*np.pi/3)*(x[1:]**3-x[:-1]**3) # spherical shell volume
count[i,:] = num # counts
e[i,:] = np.sqrt(num) # Poisson error bar
i += 1
Work out the occupation in each bin, work out volumes, work out Poisson error bars.
for i in range(3):
plt.errorbar(xp[i,:], count[i,:]/xpsq[i,:], yerr=e[i,:]/xpsq[i,:])
plt.xscale('log')
plt.yscale('log')
plt.xlim(1, 200)
plt.ylim(1e-6, 1)
plt.legend(['z=0.0','z=0.5','z=1.0'])
plt.xlabel('Radius / kpc')
plt.ylabel('dN/dV')
3. We model the density profiles using analytical representations :
from scipy.optimize import curve_fit
def nfw(x, a, rs):
return a/((x/rs)*(1. + (x/rs))**4)
# example: z=0 -> i=0
x = xp[0,:]
y = count[0,:]/xpsq[0,:]
err = e[0,:]/xpsq[0,:]
# initial guestimates for the fit
a = 0.1
rs = 12.
keep = err>0
popt, pcov = curve_fit(nfw, x[keep], y[keep], sigma=err[keep])
plt.loglog(x, y, label='Illustris')
plt.loglog(x, nfw(x,*popt), label='NFW profile')
plt.legend()
18.01 Exercise for object-oriented programming
import numpy as np
class Particle(object):
def __init__(self, mass, x, y, z, vx, vy, vz):
self.mass = mass
self.x = x
self.y = y
self.z = z
self.vx = vx
self.vy = vy
self.vz = vz
def kinetic_energy(self):
return 0.5 * self.mass * np.sqrt(self.vx ** 2 +
self.vy ** 2 +
self.vz ** 2)
def distance(self, other):
return np.sqrt((self.x - other.x) ** 2 +
(self.y - other.y) ** 2 +
(self.z - other.z) ** 2)
def move(self, dt):
self.x = self.x + self.vx * dt
self.y = self.y + self.vy * dt
self.z = self.z + self.vz * dt
class ChargedParticle(Particle):
def __init__(self, mass, x, y, z, vx, vy, vz, charge):
Particle.__init__(self, mass, x, y, z, vx, vy, vz)
self.charge = charge
p1 = Particle(2., 1., 2., 3., 1., -3., -1.)
print (p1.kinetic_energy())
p1.move(1.)
print(p1.x, p1.y, p1.z)
p1.move(1.)
print(p1.x, p1.y, p1.z)
p2 = ChargedParticle(3., 4., 3., -2., 1., 1., 1., -1.)
p2.distance(p1)