Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/ic/generate_spherical_model_nfw.py
6c0821250383master

Last Change
Raw File

generate_spherical_model_nfw.py
View Options

#!/usr/bin/env python
import Ptools as pt
from numpy import *
from pNbody import ic
from pNbody import profiles
from pNbody import libutil
from scipy import optimize
import sys
# parameters for the profile
n = 2**17
rs = 20.
rmax = 200.
M = 1.
# parameters for the profile generation
Neps_des = 10. # number of des. points in eps
ng = 256 # number of division to generate the model
rc = 0.1 # default rc (if automatic rc fails) length scale of the grid
# random params
seed = 1
# param for the plot
dR = 0.1
nr = int(rmax/dR)
pr_fct = profiles.nfw_profile
mr_fct = profiles.nfw_mr
ic_fct = ic.nfw
args = (rs,)
# compute grid parameters
Rs,rc,eps,Neps,g,gm = ic.ComputeGridParameters(n,args,rmax,M,pr_fct,mr_fct,Neps_des,rc,ng)
# create the model
ic_args = (n,)+args+(rmax,None,Rs,seed,'snap.dat','gadget')
nb = ic.nfw(n,rs,rmax,dR,Rs,"nfw.dat",ftype='gadget')
nb.write()
# check
nbs = nb.selectc(nb.rxyz()<Rs[1])
print "eps =",Rs[1]
print "des part. in eps",Neps
print "eff part. in eps",nbs.nbody
print "rc",rc
#sys.exit()
###########################
#
# plot density and Mr
#
###########################
# central density
rho0 = M/apply(mr_fct,(rmax,)+args)
# compute the profile
r = arange(dR,rmax,dR)
args = args + (rho0,)
##########################
# density
##########################
pt.subplot(2,1,1)
Rho = apply(pr_fct,(r,)+args)
r_nb,Rho_nb = nb.mdens(nb=nr,rm=rmax)
pt.plot(r,Rho,'-r')
pt.plot(r_nb,Rho_nb,'-b')
pt.semilogx()
pt.semilogy()
pt.xlabel('Radius')
pt.ylabel('Density')
##########################
# mr
##########################
pt.subplot(2,1,2)
Mr = apply(mr_fct,(r,)+args)
r_nb,Mr_nb = nb.Mr_Spherical(nr=nr,rmin=0,rmax=rmax)
pt.plot(r,Mr,'-r')
pt.plot(r_nb,Mr_nb,'-b')
pt.xlabel('Radius')
pt.ylabel('M(r)')
pt.show()
c4science · Help