Exercise solutions Day 05 (April 12th)¶
Python course summer 2024 University of Heidelberg
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.
In [3]:
#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:
In [1]:
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,num):
params = {'dm':'Coordinates'}
url = "http://www.illustris-project.org/api/Illustris-1/snapshots/" + str(num) + "/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
In [4]:
redshifts=[1.0,0.5,0.0]
snapshots=[85,103,135]
id=75 # 1030 75 543681
dm=6.3e6
for i in reversed(range(3)):
sub,filename=getsnapshot(id,snapshots[i]) # get snapshot with ID=id and given snapshot number
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 10kpc at z=',redshifts[i],
' is M={0:5.2f} x 10^10 Msun'.format(len(rr[rr<10])*dm/1e10))
Aside: Plot particle distributions
In [9]:
colors=['red','green','blue']
redshifts=[1.0,0.5,0.0]
snapshots=[85,103,135]
dict={0:[0,0],1:[0,1],2:[1,0]} # this is used to associate sub-plots with the indices 0,1,2
id=75 # 1030 75 543681
g, axarr = plt.subplots(2,2) # a plot object with axes as described in axarr
for i in range(3):
sub,filename=getsnapshot(id,snapshots[i]) # get snapshot with ID=id and given snapshot number
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'])
axarr[dict[i][0],dict[i][1]].scatter (dx,dy,c=colors[i],alpha=0.3)
axarr[dict[i][0],dict[i][1]].legend(['z='+str(redshifts[i])])
axarr[dict[i][0],dict[i][1]].set_xlim(-100,100)
axarr[dict[i][0],dict[i][1]].set_ylim(-100,100)
2. Determine the density profile $\frac{dN}{dV}$.
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
In [5]:
nbins=20
maxradius=50
xp = 0*np.ones((3,nbins))
xpsq = 0*xp
e = 0*xp
count= 0*xp
id=1030 # 1030 75 543681
snapshots=[85,103,135]
for i in range(3): # loop through three snapshots/redshifts
sub,filename=getsnapshot(id,snapshots[i]) # get json header and data file
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
""" create a histogram with 20 bins
num = counts
and x = bin boundaries
are used """
num,x,ax=plt.hist(rr,bins=nbins,range=(0,maxradius))
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
Work out the occupation in each bin, work out volumes, work out Poisson error bars.
In [6]:
%matplotlib inline
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,100)
plt.ylim(1e-6,1)
plt.legend(['z=1.0','z=0.5','z=0.0'])
plt.xlabel('kpc')
Out[6]:
3. We model the density profiles using analytical representations :
In [7]:
from scipy.optimize import curve_fit
def dark(x,a,rs):
return a/((x/rs)*(1+(x/rs))**4)
In [8]:
# example: z=0
x=xp[2,:]
y=count[2,:]/xpsq[2,:]
err=e[2,:]/xpsq[i,:]
a=0.1
rs=12
In [9]:
keep = err>0
popt, pcov = curve_fit(dark, x[keep], y[keep], sigma=err[keep])
In [10]:
plt.loglog(x, y,x,dark(x,*popt))
plt.legend(['illustris','dark'])
Out[10]:
18.01 Exercise for programming classes
In [10]:
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
In [11]:
p1 = Particle(2., 1., 2., 3., 1., -3, -1)
In [12]:
p1 = Particle(2., 1., 2., 3., 1., -3, -1)
print (p1.kinetic_energy())
In [13]:
p1.move(1.)
print(p1.x, p1.y, p1.z)
In [14]:
p1.move(1.)
print(p1.x, p1.y, p1.z)
In [15]:
import numpy as np
p2 = ChargedParticle(3., 4., 3., -2., 1, 1, 1, -1.)
p2.distance(p1)
Out[15]:
19.01
In [16]:
diogenes=['Cede','paulum','e','sole']
dic={'Cede':'Go','paulum':'a little','e':'out of','sole':'the sun'}
answer = [dic[word] for word in diogenes]
print (answer)
19.02
In [17]:
# using the solution from exercise 06.02 on day 02
def digit_sum(x):
ans=0
for digit in str(x):
ans=ans+int(digit)
return ans
answer=[num for num in range(0,1001) if digit_sum(num)==7]
print (answer)