In [1]:

    

%pylab inline


    

    




Populating the interactive namespace from numpy and matplotlib

In [2]:

    

import sys
sys.path.insert(0, '../')


    

In [3]:

    

from fastpm import core
import numpy
from fastpm.core import leapfrog

from pmesh.pm import ParticleMesh
from nbodykit.cosmology import Planck15
from kdcount import KDTree


    

In [4]:

    

from fastpm.glass import Solver


    

In [5]:

    

from nbodykit.cosmology.background import MatterDominated


    

In [ ]:

    

# NEED very high resolution to meet the shot noise limit if we want to
# start at z=99.


    

In [42]:

    

pm = ParticleMesh(BoxSize=1., Nmesh=[32, 32, 32])

def test_solver(B=2, spread=3, N=3, T=4):
    solver= Solver(pm, B)
    stages = numpy.linspace(0, numpy.pi * 2 * (N + 0.25), int(T * (N + 0.25) + 1))

    Q = pm.generate_uniform_particle_grid()
    state = solver.run(3333, Q, stages, spread=spread, N=N)
    
    return state


    

In [43]:

    

from nbodykit.lab import *


    

In [50]:

    

print('------')

state = test_solver(spread=8, N=9, B=2)
r = FFTPower(state.to_catalog(attrs={'BoxSize' :  state.pm.BoxSize}), mode='1d', Nmesh=128)
kref = 2 * pi / pm.BoxSize[0] * pm.Nmesh[0]// 4
powerlaw = (r.power['k'] ** 4 / kref ** 4)* r.power.sel(method='nearest', k=[kref])['power'].real 
plot(r.power['k'], r.power['power'] / powerlaw, label='8-9-2', lw=3)


#state = test_solver(spread=5, N=3, B=2)
#r = FFTPower(state.to_catalog(attrs={'BoxSize' :  state.pm.BoxSize}), mode='1d', Nmesh=128)
#plot(r.power['k'], r.power['power']  / powerlaw, label='5-3-2', ls=':')

#state = test_solver(spread=3, N=6, B=2)
#r = FFTPower(state.to_catalog(attrs={'BoxSize' :  state.pm.BoxSize}), mode='1d', Nmesh=128)
#plot(r.power['k'], r.power['power']  / powerlaw, label='3-6-2', ls=':')

#state = test_solver(spread=1, N=3, B=2)
#r = FFTPower(state.to_catalog(attrs={'BoxSize' :  state.pm.BoxSize}), mode='1d', Nmesh=128)
#plot(r.power['k'], r.power['power']  / powerlaw, label='1-3-2', ls=':')


state = test_solver(spread=3, N=3, B=2)
r = FFTPower(state.to_catalog(attrs={'BoxSize' :  state.pm.BoxSize}), mode='1d', Nmesh=128)
plot(r.power['k'], r.power['power'] / powerlaw, label='3-3-2 default', lw=3)

print('------')
plot(r.power['k'], powerlaw / powerlaw, label='k**4', ls='--', color='k')
plot(r.power['k'], Planck15.clone(P_k_max=r.power['k'].max() * 2).get_pklin(r.power['k'], z=99) / powerlaw, label='z=99', ls=':', color='black')
plot(r.power['k'], 1  / (pm.Nmesh / pm.BoxSize).prod()  / powerlaw, label='shot-noise', ls='-', color='black')
axvline(2 * pi / pm.BoxSize[0] * pm.Nmesh[0] / 2, label='grid nyquist', ls='--', color='grey')

legend(ncol=2)
yscale('log')
#ylim(0.1, 100)
xscale('log')
xlabel('k [kpc/h]')
ylabel('P / k**4 law')


    

    




------






    




/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/ipykernel_launcher.py:7: RuntimeWarning: invalid value encountered in true_divide
  import sys
/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/ipykernel_launcher.py:25: RuntimeWarning: invalid value encountered in true_divide
/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/ipykernel_launcher.py:28: RuntimeWarning: invalid value encountered in true_divide






    




------






    




/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/ipykernel_launcher.py:29: RuntimeWarning: invalid value encountered in true_divide
/home/yfeng1/anaconda3/install/lib/python3.6/site-packages/ipykernel_launcher.py:30: RuntimeWarning: divide by zero encountered in true_divide






    
Out[50]:





Text(0,0.5,'P / k**4 law')






    

In [ ]:

    

state = test_solver(spread=3, N=3.5, B=2, T=4)
mask = state.Q[:, 1] < 130
plot(state.X[:, 0][mask] % pm.BoxSize[0], state.X[:, 2][mask] % pm.BoxSize[0], '.')
state = test_solver(spread=3, N=3.5, B=2, T=16)
plot(state.X[:, 0][mask] % pm.BoxSize[0], state.X[:, 2][mask] % pm.BoxSize[0], 'x')


    

In [56]:

    

plot(r.power['k'], r.power['power'])
plot(r.power['k'], r.power['k'] ** 4 * r.power.sel(method='nearest', k=[1.])['power'].real)
axhline(1  / (pm.Nmesh / pm.BoxSize).prod())
axvline(2 * pi / pm.BoxSize[0] * pm.Nmesh[0] / 2)
yscale('log')
xscale('log')


    

In [46]:

    

h, b, b = histogram2d(state.X[:, 0] % pm.BoxSize[0], state.X[:, 1] % pm.BoxSize[1], bins=32)


    

In [47]:

    

hist(h.flat)


    

    
Out[47]:





(array([  1.,  11.,  75., 216., 294., 253., 126.,  40.,   7.,   1.]),
 array([22. , 24.1, 26.2, 28.3, 30.4, 32.5, 34.6, 36.7, 38.8, 40.9, 43. ]),
 <a list of 10 Patch objects>)






    

In [ ]: