如何计算任何二维多边形的重力？
Max*_*Max 4 python algorithm physics
（为简单起见，在 2d 中工作）我知道由于重力作用在两个球体上的力
G(m1*m2/r**2)
是 我最初的想法是将圆圈装入物体中，这样重力产生的力将等于每个圆圈所受力的总和。例如（伪代码），
def gravity(pos1,shape):
circles = packCircles(shape.points)
force = 0
for each circle in circles:
distance = distanceTo(pos1,circle.pos)
force += newtonForce(distance,shape.mass,1) #1 mass of observer
return force
这会是一个可行的解决方案吗？如果是这样，我将如何高效快速地打包圆圈？如果没有，是否有更好的解决方案？
编辑：注意我想如何在特定点找到物体的力，所以需要计算圆和观察者之间的角度（和矢量相加）。它不同于计算施加的总力。
背景
有些解释会有点偏离主题，但我认为有必要帮助澄清评论中提出的一些事情，因为其中大部分内容有点违反直觉。
这种引力相互作用的解释取决于点质量的概念。假设您有两个处于孤立系统中的点质量，它们彼此相隔一定距离r 1，质量分别为m 1和m 2，
m 1产生的引力场由下式给出
其中G是万有引力常数，r是m 1和r的距离？是沿m 1和m 2之间的线的单位方向。
该场施加在m 2上的重力由下式给出
注意？- ? 重要的是，这适用于任何距离的任何两点质量。1
引力相互作用的场性质允许我们在计算由多重相互作用引起的净引力时采用叠加。考虑如果我们在前面的场景中添加另一个质量m 3，
然后在大容量的重力米2是简单地从由每个所创建的领域中的重力的总和等质量
与[R I，J = R J，I。这适用于任何分离处的任何质量数。这也意味着由质量集合创建的场可以通过向量 sum聚合，如果您更喜欢那种形式主义。
现在考虑我们是否有大量的点质量M聚集在一起，形成一个连续的、密度均匀的刚体。然后我们想计算单个空间上不同的点质量m上的重力，这是由于总质量M：
然后，我们可以考虑不同大小的质量区域（或体积）而不是考虑点质量，并对这些区域（或体积）对点质量的影响进行积分或求和。在二维情况下，引力的大小为
在哪里？是总质量的密度。2这相当于对每个微分质量Δdxdy产生的重力矢量场求和。这样的等价是非常重要的，因为这意味着任何一点质量远远不够外一个质量分布，重力由于质量分布几乎是完全相同一样，因为这将是大规模的点质量中号位于中心质量的质量分布。3 4
这意味着，为了非常好地近似，当涉及由于任何质量分布而计算引力场时，质量分布可以用分布质心处的等效质量点质量代替。这适用于任意数量的空间上不同的质量分布，无论这些分布是否构成刚体。此外，这意味着您甚至可以将分布组聚合为系统质心处的单个点质量。5 只要参考点距离足够远。
然而，为了找到由于任何点的质量分布而导致的点质量上的引力，对于形状和分离不可知方式的任何质量分布，我们必须通过对每个部分的贡献求和来计算该点的引力场的质量分布。6
回到问题
当然，对于任意多边形或多面体，解析解可能非常困难，因此使用求和要简单得多，算法方法也将类似地使用求和。
从算法上讲，这里最简单的方法实际上不是几何包装（圆形/球体或正方形/立方体）。使用打包并非不可能，但从数学上讲，这种方法存在重大挑战 - 最好采用依赖于更简单数学的方法。一种这样的方法是定义一个包含质量分布空间范围的网格，然后以网格点为顶点创建简单的（正方形/立方体或矩形/立方体）多边形或多面体。这将创建三种多边形或多面体：
- 那些不包含任何质量分布的
- 那些被质量分布完全填满的
- 那些被质量分布部分填充的
质心 - 方法 1
当从参考点到质量分布的距离相对于分布的角度范围较大时，并且当质量分布（或任何几个分布）没有参考的几何包围时，这将很有效。
然后，您可以通过对每个多边形的贡献求和来找到分布的质心R，
其中M是分布的总质量，r i是到第i个多边形几何中心的空间向量，而m i是密度乘以包含质量的多边形部分（即1.00表示完全填充的多边形，0.00对于完全空的多边形）。随着您增加采样大小（网格点的数量），质心的近似值将接近解析解。一旦你有了质心，计算产生的引力场就很简单了：你只需在上面的R点方程上放置一个质量为M的点质量。并使用
为了演示，这里是使用多边形操作的shapely库在 Python 中在二维中实现所描述的方法：
import numpy as np
import matplotlib.pyplot as plt
import shapely.geometry as geom
def centerOfMass(r, density = 1.0, n = 100):
theta = np.linspace(0, np.pi*2, len(r))
xy = np.stack([np.cos(theta)*r, np.sin(theta)*r], 1)
mass_dist = geom.Polygon(xy)
x, y = mass_dist.exterior.xy
# Create the grid and populate with polygons
gx, gy = np.meshgrid(np.linspace(min(x), max(x), n), np.linspace(min(y),
max(y), n))
polygons = [geom.Polygon([[gx[i,j], gy[i,j]],
[gx[i,j+1], gy[i,j+1]],
[gx[i+1,j+1],gy[i+1,j+1]],
[gx[i+1,j], gy[i+1,j]],
[gx[i,j], gy[i,j]]])
for i in range(gx.shape[0]-1) for j in range(gx.shape[1]-1)]
# Calculate center of mass
R = np.zeros(2)
M = 0
for p in polygons:
m = (p.intersection(mass_dist).area / p.area) * density
M += m
R += m * np.array([p.centroid.x, p.centroid.y])
return geom.Point(R / M), M
density = 1.0 # kg/m^2
G = 6.67408e-11 # m^3/kgs^2
theta = np.linspace(0, np.pi*2, 100)
r = np.cos(theta*2+np.pi)+5+np.sin(theta)+np.cos(theta*3+np.pi/6)
R, M = centerOfMass(r, density)
m = geom.Point(20, 0)
r_1 = m.distance(R)
m_1 = 5.0 # kg
F = G * (m_1 * M) / r_1**2
rhat = np.array([R.x - m.x, R.y - m.y])
rhat /= (rhat[0]**2 + rhat[1]**2)**0.5
# Draw the mass distribution and force vector, etc
plt.figure(figsize=(12, 6))
plt.axis('off')
plt.plot(np.cos(theta)*r, np.sin(theta)*r, color='k', lw=0.5, linestyle='-')
plt.scatter(m.x, m.y, s=20, color='k')
plt.text(m.x, m.y-1, r'$m$', ha='center')
plt.text(1, -1, r'$M$', ha='center')
plt.quiver([m.x], [m.y], [rhat[0]], [rhat[1]], width=0.004,
scale=0.25, scale_units='xy')
plt.text(m.x - 5, m.y + 1, r'$F = {:.5e}$'.format(F))
plt.scatter(R.x, R.y, color='k')
plt.text(R.x, R.y+0.5, 'Center of Mass', va='bottom', ha='center')
plt.gca().set_aspect('equal')
plt.show()
这种方法有点矫枉过正：在大多数情况下，找到质心和多边形面积乘以质心和总质量的密度就足够了。然而，它甚至适用于不均匀的质量分布——这就是我使用它进行演示的原因。
现场求和 - 方法 2
在许多情况下，这种方法也是矫枉过正，尤其是与第一种方法相比，但它会在任何分布下（在经典范围内）提供最佳近似值。
这里的想法是将质量分布的每一块对一个点质量的影响相加，以确定净引力（基于引力场可以独立添加的前提）
class pointMass:
def __init__(self, mass, x, y):
self.mass = mass
self.x = x
self.y = y
density = 1.0 # kg/m^2
G = 6.67408e-11 # m^3/kgs^2
def netForce(r, m1, density = 1.0, n = 100):
theta = np.linspace(0, np.pi*2, len(r))
xy = np.stack([np.cos(theta)*r, np.sin(theta)*r], 1)
# Create a shapely polygon for the mass distribution
mass_dist = geom.Polygon(xy)
x, y = mass_dist.exterior.xy
# Create the grid and populate with polygons
gx, gy = np.meshgrid(np.linspace(min(x), max(x), n), np.linspace(min(y),
max(y), n))
polygons = [geom.Polygon([[gx[i,j], gy[i,j]],
[gx[i,j+1], gy[i,j+1]],
[gx[i+1,j+1],gy[i+1,j+1]],
[gx[i+1,j], gy[i+1,j]],
[gx[i,j], gy[i,j]]])
for i in range(gx.shape[0]-1) for j in range(gx.shape[1]-1)]
g = np.zeros(2)
for p in polygons:
m2 = (p.intersection(mass_dist).area / p.area) * density
rhat = np.array([p.centroid.x - m1.x, p.centroid.y - m1.y])
rhat /= (rhat[0]**2 + rhat[1]**2)**0.5
g += m1.mass * m2 / p.centroid.distance(geom.Point(m1.x, m1.y))**2 * rhat
g *= G
return g
theta = np.linspace(0, np.pi*2, 100)
r = np.cos(theta*2+np.pi)+5+np.sin(theta)+np.cos(theta*3+np.pi/6)
m = pointMass(5.0, 20.0, 0.0)
g = netForce(r, m)
plt.figure(figsize=(12, 6))
plt.axis('off')
plt.plot(np.cos(theta)*r, np.sin(theta)*r, color='k', lw=0.5, linestyle='-')
plt.scatter(m.x, m.y, s=20, color='k')
plt.text(m.x, m.y-1, r'$m$', ha='center')
plt.text(1, -1, r'$M$', ha='center')
ghat = g / (g[0]**2 + g[1]**2)**0.5
plt.quiver([m.x], [m.y], [ghat[0]], [ghat[1]], width=0.004,
scale=0.25, scale_units='xy')
plt.text(m.x - 5, m.y + 1, r'$F = ({:0.3e}, {:0.3e})$'.format(g[0], g[1]))
plt.gca().set_aspect('equal')
plt.show()
对于所使用的相对简单的测试用例，其中的结果与第一种方法非常接近：
但是，虽然存在第一种方法无法正常工作的情况，但不会出现第二种方法会失败的情况（在经典体制中），因此建议采用这种方法。
1这在极端情况下确实会崩溃，例如超过黑洞的事件视界，或者当r接近普朗克长度时，但这些情况不是本问题的主题。
2在密度不均匀的情况下，这会变得更加复杂，并且在质量分布无法用符号描述的情况下没有简单的解析解。
3应该注意的是，积分实际上是这样做的；找到质心。
4对于质点内质量分布牛顿壳牌定理，或现场总和必须使用。
5在天文学中，这被称为重心，物体总是围绕系统的重心运行——而不是任何给定物体的质心。
6在某些情况下，使用牛顿壳定理就足够了，但是这些情况与分布几何无关。