diff --git a/nbody/numba/nbody/barnes_hut_array/numba_functions.py b/nbody/numba/nbody/barnes_hut_array/numba_functions.py index 0913627..f380bf0 100644 --- a/nbody/numba/nbody/barnes_hut_array/numba_functions.py +++ b/nbody/numba/nbody/barnes_hut_array/numba_functions.py @@ -9,7 +9,7 @@ def buildTree(center0, box_size0, child, cell_center, cell_radius, particles): nbodies = particles.shape[0] for ip in range(nbodies): center = center0.copy() - box_size = box_size0.copy() + box_size = box_size0 x, y = particles[ip, :2] cell = 0 @@ -44,15 +44,15 @@ def buildTree(center0, box_size0, child, cell_center, cell_radius, particles): ncell += 1 child[childIndex] = nbodies + ncell center[:] = cell_center[cell] - box_size[:] = .5*cell_radius[cell] + box_size = .5*cell_radius[cell] if (oldchildPath&1): - center[0] += box_size[0] + center[0] += box_size else: - center[0] -= box_size[0] + center[0] -= box_size if ((oldchildPath>>1)&1): - center[1] += box_size[1] + center[1] += box_size else: - center[1] -= box_size[1] + center[1] -= box_size oldchildPath = 0 if particles[npart, 0] > center[0]: @@ -110,9 +110,9 @@ def buildTree(center0, box_size0, child, cell_center, cell_radius, particles): # acc += np.array([F * dx, F * dy]) # else: # dist = np.sqrt(dx**2 + dy**2) -# #print(dist, localNode[depth], self.cell_radius[child - self.nbodies][0], -# # self.cell_radius[child - self.nbodies][0]/dist) -# if dist != 0 and cell_radius[child - nbodies][0]/dist <.5: +# #print(dist, localNode[depth], self.cell_radius[child - self.nbodies], +# # self.cell_radius[child - self.nbodies]/dist) +# if dist != 0 and cell_radius[child - nbodies]/dist <.5: # #print('dist', dx, dy) # F = gamma_si*mass[child]/(dist*dist*dist) # acc += np.array([F * dx, F * dy]) @@ -146,7 +146,7 @@ def computeForce(nbodies, child_array, center_of_mass, mass, cell_radius, p): dx = center_of_mass[child, 0] - pos[0] dy = center_of_mass[child, 1] - pos[1] dist = np.sqrt(dx**2 + dy**2) - if dist != 0 and cell_radius[child - nbodies][0]/dist <.5: + if dist != 0 and cell_radius[child - nbodies]/dist <.5: Fx, Fy = force(pos, center_of_mass[child], mass[child]) acc[0] += Fx acc[1] += Fy diff --git a/nbody/numba/nbody/barnes_hut_array/quadTree.py b/nbody/numba/nbody/barnes_hut_array/quadTree.py index aa06103..ea37be2 100644 --- a/nbody/numba/nbody/barnes_hut_array/quadTree.py +++ b/nbody/numba/nbody/barnes_hut_array/quadTree.py @@ -9,12 +9,12 @@ def __init__(self, bmin, bmax, size): self.bmin = np.asarray(bmin) self.bmax = np.asarray(bmax) self.center = .5*(self.bmin + self.bmax) - self.box_size = (self.bmax - self.bmin) + self.box_size = (self.bmax - self.bmin).max() self.ncell = 0 self.cell_center = np.zeros((2*size+1, 2)) - self.cell_radius = np.zeros((2*size+1, 2)) + self.cell_radius = np.zeros(2*size+1) self.cell_center[0] = self.center - self.cell_radius[0] = self.box_size + self.cell_radius[0] = 0.5*self.box_size def buildTree(self, particles): self.ncell = numba_functions.buildTree(self.center, self.box_size, self.child, self.cell_center, self.cell_radius, particles) diff --git a/nbody/python/nbody/barnes_hut_array/quadTree.py b/nbody/python/nbody/barnes_hut_array/quadTree.py index 7b1de44..e317c02 100644 --- a/nbody/python/nbody/barnes_hut_array/quadTree.py +++ b/nbody/python/nbody/barnes_hut_array/quadTree.py @@ -8,17 +8,17 @@ def __init__(self, bmin, bmax, size): self.bmin = np.asarray(bmin) self.bmax = np.asarray(bmax) self.center = .5*(self.bmin + self.bmax) - self.box_size = (self.bmax - self.bmin) + self.box_size = (self.bmax - self.bmin).max() self.ncell = 0 self.cell_center = np.zeros((2*size+1, 2)) - self.cell_radius = np.zeros((2*size+1, 2)) + self.cell_radius = np.zeros(2*size+1) self.cell_center[0] = self.center - self.cell_radius[0] = self.box_size + self.cell_radius[0] = 0.5*self.box_size def buildTree(self, particles): for ip, p in enumerate(particles): center = self.center.copy() - box_size = self.box_size.copy() + box_size = self.box_size x, y = p[:2] cell = 0 @@ -55,13 +55,13 @@ def buildTree(self, particles): center[:] = self.cell_center[cell] box_size = .5*self.cell_radius[cell] if (oldchildPath&1): - center[0] += box_size[0] + center[0] += box_size else: - center[0] -= box_size[0] + center[0] -= box_size if ((oldchildPath>>1)&1): - center[1] += box_size[1] + center[1] += box_size else: - center[1] -= box_size[1] + center[1] -= box_size oldchildPath = 0 if particles[npart][0] > center[0]: @@ -125,7 +125,7 @@ def computeForce(self, p): dx = self.center_of_mass[child, 0] - pos[0] dy = self.center_of_mass[child, 1] - pos[1] dist = np.sqrt(dx**2 + dy**2) - if dist != 0 and self.cell_radius[child - self.nbodies][0]/dist <.5: + if dist != 0 and self.cell_radius[child - self.nbodies]/dist <.5: F = force(pos, self.center_of_mass[child], self.mass[child]) acc += F else: @@ -146,4 +146,4 @@ def __str__(self): s += 2*indent + 'particules: {p}\n'.format(p=cellElements[np.logical_and(0<=cellElements, cellElements=self.nbodies]-self.nbodies) - return s \ No newline at end of file + return s diff --git a/nbody/python/nbody/init.py b/nbody/python/nbody/init.py index dbfe81b..d6e4657 100644 --- a/nbody/python/nbody/init.py +++ b/nbody/python/nbody/init.py @@ -1,5 +1,5 @@ import numpy as np -from .physics import gamma_1 +from .physics import gamma_1, gamma_si def init_solar_system(): bodies = np.array([[ 0, 0, 0, 0], #sun @@ -33,7 +33,7 @@ def getOrbitalVelocity(xb, yb, mb, xs, ys): dist = np.sqrt(r[0] * r[0] + r[1] * r[1]) # Based on the distance from the sun calculate the velocity needed to maintain a circular orbit - v = np.sqrt(gamma_1 * mb / dist) + v = np.sqrt(gamma_si * mb / dist) # Calculate a suitable vector perpendicular to r for the velocity of the tracer vxs = ( r[1] / dist) * v @@ -65,7 +65,7 @@ def init_collisions(blackHole): nstars = b['stars'] rad = b['radstars'] - r = 0.1 + .8 * (rad * np.random.rand(nstars)) + r = 0.3 + .8 * (rad * np.random.rand(nstars)) a = 2*np.pi*np.random.rand(nstars) tmp_mass = 0.03 + 20*np.random.rand(nstars) x = b['coord'][0] + r*np.sin(a) @@ -75,8 +75,8 @@ def init_collisions(blackHole): particles[ind:ind+nstars, 0] = x particles[ind:ind+nstars, 1] = y - particles[ind:ind+nstars, 2] = 1e2 * (vx + vxb) - particles[ind:ind+nstars, 3] = 1e2 * (vy + vyb) + particles[ind:ind+nstars, 2] = vx + vxb + particles[ind:ind+nstars, 3] = vy + vyb mass[ind:ind+nstars] = tmp_mass ind += nstars