我一直在编写一个递归四叉树构造函数以用于某些n体仿真,但是我当前的实现似乎无法正常工作,并且在进行大量调试之后,我很沮丧。尽管所有调试检查似乎都给出了应有的结果,但它给出的结果显然是错误的。有人可以帮忙吗?这是我的代码:
import numpy as np
import matplotlib.pyplot as plt
from numba import jit, autojit
import matplotlib.animation as ani
from mpl_toolkits.mplot3d import Axes3D
from time import sleep
mu, sigma = 0, 0.5
size = 50
X = np.array([np.random.normal(mu, sigma, size),np.random.normal(mu, sigma, size)])
A = np.zeros([size,2],dtype='float64')
V = np.zeros([size,2],dtype='float64')
M = np.random.normal(1, 0, size)
quad_list = np.zeros([4,1])
def quadtree(p,n,x,y,w,h):
global quad_list
L = len(p[0])
if L>1:
print(len(p[0]))
midx = (w/2+x)
midy = (h/2+y)
px = p[0]
py = p[1]
plx, ply = px[px<midx], py[px<midx]
prx, pry = px[px>midx], py[px>midx]
p1 = np.array([plx[ply>midy],ply[ply>midy]])
p2 = np.array([prx[pry>midy],pry[pry>midy]])
p3 = np.array([prx[pry<midy],pry[pry<midy]])
p4 = np.array([plx[ply<midy],ply[ply<midy]])
quad_list = np.append(quad_list,np.array([x,y,w,h]))
quadtree(p1,n+1,x,y+h/2,w/2,h/2)
quadtree(p2,n+1,x+w/2,y+h/2,w/2,h/2)
quadtree(p3,n+1,x,y,w/2,h/2)
quadtree(p4,n+1,x,y,w/2,h/2)
else:
quad_list = np.append(quad_list,np.array([x,y,w,h]))
quadtree(X,0,-2,-2,4,4)
plt.scatter(X[0],X[1],c='black')
out = np.zeros([4,int(len(quad_list)/4)])
for i in range(0,int(len(quad_list)),4):
for j in range(4):
out[j,int(i/4)] = quad_list[i+j]
for n in range(int(len(out[0,:]))):
print(n)
plt.plot([out[0,n],out[0,n]], [out[1,n],out[1,n]+out[3,n]])
plt.plot([out[0,n]+out[2,n],out[0,n]+out[2,n]], [out[1,n],out[1,n]+out[3,n]])
plt.plot([out[0,n],out[0,n]+out[2,n]], [out[1,n],out[1,n]])
plt.plot([out[0,n],out[0,n]+out[2,n]], [out[1,n]+out[3,n],out[1,n]+out[3,n]])
plt.show()
谢谢!(注意:我正在制作一个更优化的版本,其中quad_list
是列表)
尝试更改行:
quadtree(p3,n+1,x,y,w/2,h/2)
to
quadtree(p3,n+1,x+w/2,y,w/2,h/2)