Homec4science
Diffusion pNbody (master) Browse

pNbody/src/PyGadget/test_SphEvaluate.py
021ca11fcae0master

Last Change
Raw File

test_SphEvaluate.py
View Options

#!/usr/bin/env python
from mpi4py import MPI
from pNbody import ic
from numpy import *
import gadget
import Ptools as pt
import sys
import time
random.seed(MPI.COMM_WORLD.Get_rank()) # 2 for two points
n = 10000
rc = 2
rmax = 100
nb = ic.plummer(n,1,1,1,eps=rc,rmax=rmax,ftype='gadget')
nb.rename('plummer.dat')
nb.write()
gadget.InitMPI() # init MPI
gadget.InitDefaultParameters() # init default parameters
gadget.Info()
params = {}
params['ErrTolTheta'] = 0.7
params['DesNumNgb'] = 50
params['MaxNumNgbDeviation'] = 1
params['UnitLength_in_cm'] = 3.085e+21
params['UnitMass_in_g'] = 4.435693e+44
params['UnitVelocity_in_cm_per_s'] = 97824708.2699
params['SofteningGas'] = rc/100.
params['SofteningHalo'] = rc/100.
params['SofteningDisk'] = rc/100.
params['SofteningBulge'] = rc/100.
params['SofteningStars'] = rc/100.
params['SofteningBndry'] = rc/100.
params['SofteningGasMaxPhys'] = rc/100.
params['SofteningHaloMaxPhys'] = rc/100.
params['SofteningDiskMaxPhys'] = rc/100.
params['SofteningBulgeMaxPhys'] = rc/100.
params['SofteningStarsMaxPhys'] = rc/100.
params['SofteningBndryMaxPhys'] = rc/100.
gadget.SetParameters(params)
params = gadget.GetParameters()
gadget.LoadParticles(array(nb.npart),nb.pos,nb.vel,nb.mass,nb.num,nb.tpe)
nb.rho = gadget.GetAllDensities()
nb.hsml = gadget.GetAllHsml()
nb.pos = gadget.GetAllPositions()
nb.obs = nb.x()**2 # define an observable
x = arange(-50,50,5)
y = zeros(len(x),float32)
z = zeros(len(x),float32)
pos = transpose(array([x,y,z])).astype(float32)
r = sqrt(x**2 + y**2 + z**2)
hsml = ones(len(pos)).astype(float32) * nb.hsml.mean()
print "start SphEvaluate"
rho,hsml=gadget.InitHsml(pos,hsml)
obs=gadget.SphEvaluate(pos,hsml,nb.obs)
###############################
# plot
###############################
if MPI.COMM_WORLD.Get_rank()==0:
color = 'r'
if MPI.COMM_WORLD.Get_rank()==1:
color = 'b'
pt.scatter(nb.x(),nb.obs,s=1,color='k')
pt.scatter(x,obs,s=10,color=color)
#pt.figure()
#pt.scatter(nb.rxyz(),nb.rho,s=1)
#pt.scatter(r,rho,s=10,color='r')
##pt.semilogx()
##pt.semilogy()
pt.show()
c4science · Help