from pNbody import *

nb = Nbody('disk.dat',ftype='gadget')
nb.info()
nb.show()


# compute acceleration

nb.getTree()
nb.Tree.info()
eps = 0.25
nb.TreeAccel(array([[0,0,0]]).astype(Float32),eps)
acc = nb.TreeAccel(nb.pos,eps)


# compute potential
R = arange(0,200,1.0).astype(Float32)
pos = concatenate((R,zeros(len(R)),zeros(len(R))))
pos.shape = (3,len(R))
pos = transpose(pos).astype(Float32)

pot = nb.TreePot(pos,eps)

import SM
g = SM.plot()
g.erase()
g.limits(R,pot)
g.box()
g.connect(R,pot)
g.show()


# compute sph
nb.InitSphParameters(33,3)
nb.ComputeSph()
nb.Density
nb.Hsml


val = nb.SphEvaluate(nb.Vr())
rot = nb.SphEvaluate('rot')
div = nb.SphEvaluate('div')


# for a specific place
Hsml = array([100]).astype(Float32)
pos  = array([[0,0,0]]).astype(Float32)
Density,Hsml = nb.ComputeDensityAndHsml(pos,Hsml)
Density,Hsml
Density,Hsml = nb.ComputeDensityAndHsml(pos,Hsml)
Density,Hsml
nb.SphEvaluate(nb.Vr(),pos=pos,vel=None,hsml=Hsml)



# gwin -t gadget disk.dat
self.nb.getTree()
self.nb.acc = self.nb.TreeAccel(self.nb.pos,0.25)

# parallelisme 
mpd &
mpirun -np 2 gwin -t gadget disk.dat



