pNbody tutorial

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()