Homec4science
Diffusion pNbody (master) Browse

pNbody/examples/profiles/generate_spherical_profiles.py
6c0821250383master

Last Change
Raw File

generate_spherical_profiles.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.
gamma = 1.5
a = 2.
b = 3.
rmax = 200.
M = 1.
# parameters for the profile generation
eps_des = 0.1 # desired eps (not used)
Neps_des = 10. # number of des. poin in eps
ng = 256 # number of division to generate the model
rc = 0.1 # default rc (if automatic rc fails)
# random params
seed = 1
# param for the plot
dR = 0.1
nr = int(rmax/dR)
# model
#pr_fct = profiles.hernquist_profile
#mr_fct = profiles.hernquist_mr
#ic_fct = ic.hernquist
#args = (rs,)
#pr_fct = profiles.burkert_profile
#mr_fct = profiles.burkert_mr
#ic_fct = ic.burkert
#args = (rs,)
#pr_fct = profiles.nfw_profile
#mr_fct = profiles.nfw_mr
#ic_fct = ic.nfw
#args = (rs,)
#pr_fct = profiles.nfwg_profile
#mr_fct = profiles.nfwg_mr
#ic_fct = ic.nfwg
#args = (rs,gamma)
pr_fct = profiles.generic2c_profile
mr_fct = profiles.generic2c_mr
ic_fct = ic.generic2c
args = (rs,a,b)
# compute grid parameters
Rs,eps,Neps = 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 = apply(ic_fct,ic_args)
nb.write()
# check
nbs = nb.selectc(nb.rxyz()<eps)
print "rc",rc,Rs[1]
print "des part. in eps",Neps
print "eff part. in eps",nbs.nbody
print "eps/eps_des",eps/eps_des
###########################
#
# 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