pNbody tutorial¶
- official website : http://obswww.unige.ch/~revaz/pNbody/
- tutorial : http://obswww.unige.ch/~revaz/pNbody/rst/Tutorial.html
Using pNbody in the interpreter¶
First, load the python module:
>>> from pNbody import *
Create an empty pNbody object and get info
>>> nb = Nbody()
>>> nb.info()
-----------------------------------
particle file : ['file.dat']
ftype : 'Nbody_default'
mxntpe : 6
nbody : 0
nbody_tot : 0
npart : [0, 0, 0, 0, 0, 0]
npart_tot : [0, 0, 0, 0, 0, 0]
mass_tot : 0.0
byteorder : 'little'
pio : 'no'
>>> nb.nbody
0
Create a pNbody object from positions:
>>> pos = ones((10,3),float32)
>>> nb = Nbody(pos=pos)
>>> nb.info()
-----------------------------------
particle file : ['file.dat']
ftype : 'Nbody_default'
mxntpe : 6
nbody : 10
nbody_tot : 10
npart : [10, 0, 0, 0, 0, 0]
npart_tot : [10, 0, 0, 0, 0, 0]
mass_tot : 1.0000001
byteorder : 'little'
pio : 'no'
len pos : 10
pos[0] : array([ 1., 1., 1.], dtype=float32)
pos[-1] : array([ 1., 1., 1.], dtype=float32)
len vel : 10
vel[0] : array([ 0., 0., 0.], dtype=float32)
vel[-1] : array([ 0., 0., 0.], dtype=float32)
len mass : 10
mass[0] : 0.1
mass[-1] : 0.1
len num : 10
num[0] : 0
num[-1] : 9
len tpe : 10
tpe[0] : 0
tpe[-1] : 0
>>> nb.nbody
10
>>> nb.pos
array([[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.],
[ 1., 1., 1.]], dtype=float32)
>>>
Create a pNbody object using the ic module:
>>> from pNbody import ic
>>> nb = ic.plummer(10000,1,1,1,1,rmax=20.,M=1,name='plummer.dat',ftype='gadget')
comovingintegration= False
>>> nb.write()
[10000, 0, 0, 0, 0, 0]
Read a pNbody from a file¶
>>> nb = Nbody('plummer.dat',ftype='gadget')
reading pos...
reading vel...
reading num...
comovingintegration= False
>>> nb.info()
-----------------------------------
particle file : ['plummer.dat']
ftype : 'Nbody_gadget'
mxntpe : 6
nbody : 10000
nbody_tot : 10000
npart : [10000, 0, 0, 0, 0, 0]
npart_tot : [10000, 0, 0, 0, 0, 0]
mass_tot : 1.0000535
byteorder : 'little'
pio : 'no'
len pos : 10000
pos[0] : array([ 0.08620726, -0.88255513, -0.68550044], dtype=float32)
pos[-1] : array([-0.68355829, 2.11995625, -0.13416442], dtype=float32)
len vel : 10000
vel[0] : array([ 0., 0., 0.], dtype=float32)
vel[-1] : array([ 0., 0., 0.], dtype=float32)
len mass : 10000
mass[0] : 9.9999997e-05
mass[-1] : 9.9999997e-05
len num : 10000
num[0] : 0
num[-1] : 9999
len tpe : 10000
tpe[0] : 0
tpe[-1] : 0
atime : 0.0
redshift : 0.0
flag_sfr : 0
flag_feedback : 0
nall : [10000, 0, 0, 0, 0, 0]
flag_cooling : 0
num_files : 1
boxsize : 0.0
omega0 : 0.0
omegalambda : 0.0
hubbleparam : 0.0
flag_age : 0
flag_metals : 0
nallhw : [0, 0, 0, 0, 0, 0]
flag_entr_ic : 0
critical_energy_spec: 0.0
>>> nb.nbody
10000
>>> nb.npart
[10000, 0, 0, 0, 0, 0]
Get particles properties¶
Particles properties are numpy arrays of size nb.nbody:
>>> nb.pos
array([[ 0.08620726, -0.88255513, -0.68550044],
[-0.71585971, 1.87051642, -0.16716339],
[-0.0029057 , -0.03153811, -0.03678056],
...,
[ 0.40189224, 0.11727811, 0.06398547],
[-0.06060761, -0.17355472, -0.50859088],
[-0.68355829, 2.11995625, -0.13416442]], dtype=float32)
>>> nb.vel
array([[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.],
...,
[ 0., 0., 0.],
[ 0., 0., 0.],
[ 0., 0., 0.]], dtype=float32)
>>> nb.mass
array([ 9.99999975e-05, 9.99999975e-05, 9.99999975e-05, ...,
9.99999975e-05, 9.99999975e-05, 9.99999975e-05], dtype=float32)
>>> nb.num
array([ 0, 1, 2, ..., 9997, 9998, 9999])
>>> nb.tpe
array([0, 0, 0, ..., 0, 0, 0])
Get the list of particles properties:
>>> nb.get_list_of_array()
['mass', 'num', 'pos', 'tpe', 'vel']
Display a model¶
>>> nb.show(size=(10,10))
>>> nb.show(size=(20,20),shape=(256,256),palette='rainbow4')
Get physical values¶
positions x coordinates:
>>> nb.x()
array([ 0.08620726, -0.71585971, -0.0029057 , ..., 0.40189224,
-0.06060761, -0.68355829], dtype=float32)
center of mass:
>>> nb.cm()
array([ 0.0013582 , -0.04152351, 0.01155748])
angular momentum:
>>> nb.Ltot()
array([ 0., 0., 0.], dtype=float32)
compute potential en acceleration:
>>> r = nb.rxyz()
>>> pot = nb.Pot(nb.pos,eps=0.1)
>>> acc = nb.Accel(nb.pos,eps=0.1)
alternatively, you can use a treecode algorithm to compute the potential and acceleration (which is much faster):
>>> pot = nb.TreePot(nb.pos,eps=0.1)
>>> acc = nb.TreeAccel(nb.pos,eps=0.1)
plot:
>>> import pylab as pt
>>> pt.scatter(r,pot)
>>> pt.show()