Create a log-normal mock dark-matter distribution

%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np

import powerbox as pbox
pbox.__version__

'0.5.5'

In this demo, we create a mock dark-matter distribution, based on the cosmological power spectrum. To generate the power-spectrum we use the hmf code (https://github.com/steven-murray/hmf).

import hmf
hmf.__version__

'3.0.3'

The box can be set up like this:

from scipy.interpolate import InterpolatedUnivariateSpline as spline
import numpy as np

# Set up a MassFunction instance to access its cosmological power-spectrum
mf = hmf.MassFunction(z=0)

# Generate a callable function that returns the cosmological power spectrum.
spl = spline(np.log(mf.k),np.log(mf.power),k=2)
power = lambda k : np.exp(spl(np.log(k)))

# Create the power-box instance. The boxlength is in inverse units of the k of which pk is a function, i.e.
# Mpc/h in this case.
pb = pbox.LogNormalPowerBox(N=256, dim=3, pk = power, boxlength= 500., seed=1234)

We will be using the delta_x quantity a couple of times throughout this demo, so we save it to a variable here. Otherwise, powerbox will recalculate it each time it is called.

deltax = pb.delta_x()

Now we can make a plot of a slice of the density field:

plt.imshow(np.mean(deltax[:100,:,:],axis=0),extent=(0,100,0,100))
plt.colorbar()
plt.show()

And we can also compare the power-spectrum of the output field to the input power:

p_k, kbins = pbox.get_power(deltax,pb.boxlength)
plt.plot(mf.k,mf.power,label="Input Power")
plt.plot(kbins,p_k,label="Sampled Power")
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.show()

Furthermore, we can sample a set of discrete particles on the field and plot them:

particles = pb.create_discrete_sample(nbar=0.003,min_at_zero=True)

plt.figure(figsize=(8,8))
plt.scatter(particles[:,0],particles[:,1],s=np.sqrt(100./particles[:,2]),alpha=0.2)
plt.xlim(0,500)
plt.ylim(0,500)
plt.show()

Or plot them in 3D!

from mpl_toolkits.mplot3d import Axes3D

fig = plt.figure(figsize=(10,10))
ax = fig.add_subplot(111, projection='3d')

ax.scatter(particles[:,0], particles[:,1], particles[:,2],s=1,alpha=0.2)
plt.show()

Then check that the power-spectrum of the sample matches the input:

p_k_sample, kbins_sample = pbox.get_power(particles, pb.boxlength,N=pb.N)

plt.plot(mf.k,mf.power,label="Input Power")
plt.plot(kbins_sample,p_k_sample,label="Sampled Power Discrete")
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.show()