Python Plummer 模型中的随机初始速度

Question

我想在体积内随机分布 N 个粒子，使其满足 Plummer 势能分布。我试图从 Hut 的 "The Art of Computational Science" 开始工作，它有一个描述，但我似乎无法实现它。我与 Hut 的不同之处在于我需要为每个粒子提供 3 个速度分量。这是我到目前为止所做的：

f=0
g=0.1
if g >f*f*(1-f*f)**3.5:
    f=np.random.uniform(0,1,N)
    g=np.random.uniform(0,0.1,N)

vel_x= f*np.sqrt(2)*(1+x*x)**(-0.25)
vel_y= f*np.sqrt(2)*(1+y*y)**(-0.25)
vel_z= f*np.sqrt(2)*(1+z*z)**(-0.25)

vel = np.zeros((N,3))
vel[:,0]=vel_x
vel[:,1]=vel_y
vel[:,2]=vel_z

但是，当我运行 Hut 描述的能量检查时，动能为 ~0.147 N body 单位，此代码失败。任何关于我哪里出错的建议将不胜感激

Answer 1

您可能误读了 Hut 书中的 Ruby 代码，因为它也生成 3 维速度向量：

x = 0.0
y = 0.1
while y > x*x*(1.0-x*x)**3.5
  x = frand(0,1)
  y = frand(0,0.1)
end
velocity = x * sqrt(2.0) * ( 1.0 + radius*radius)**(-0.25)
theta = acos(frand(-1, 1))
phi = frand(0, 2*PI)
b.vel[0] = velocity * sin( theta ) * cos( phi )
b.vel[1] = velocity * sin( theta ) * sin( phi )
b.vel[2] = velocity * cos( theta )

第一部分生成|v|通过从速度分布中进行拒绝采样。第二部分在 space（极坐标）中生成随机方向，代码的最后一部分从极坐标转换为笛卡尔坐标。

您的代码做了完全不同的事情。您应该改写上面 Python 中显示的代码片段，例如：

f = 0.0
g = 0.1
while g > f*f*(1.0-f*f)**3.5:
  f = np.random.uniform(0,1)
  g = np.random.uniform(0,0.1)
velocity = f * np.sqrt(2.0) * (1.0 + radius*radius)**(-0.25)
theta = np.arccos(np.random.uniform(-1, 1))
phi = np.random.uniform(0, 2*np.pi)
vel[n,0] = velocity * np.sin(theta) * np.cos(phi)
vel[n,1] = velocity * np.sin(theta) * np.sin(phi)
vel[n,2] = velocity * np.cos(theta)

代码可能被矢量化，但实际上它没有什么意义，因为拒绝采样不可矢量化（它可能并且很可能会对每个样本进行不同次数的迭代）。

Python Plummer 模型中的随机初始速度

Random initial velocities within Plummer model in Python

python

physics