我怎样才能算出任何二维多边形的引力？
How can I work out the gravitational force of any 2d polygon?
（为简单起见，在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
这是一个可行的解决方案吗？如果是这样，我将如何高效快速地打包圈子？如果不行，有没有更好的解决办法？
编辑：请注意我想如何找到物体在特定点的力，因此需要计算圆和观察者之间的角度（和向量求和）。它不同于计算所施加的总力。
背景
有些解释会有些跑题，但我认为有必要帮助澄清评论中提出的一些事情，因为其中大部分内容有点违反直觉。
这种引力相互作用的解释依赖于point masses的概念。假设你有两个质点，它们在一个孤立的系统中，彼此相隔一定距离，r1，质量为 m1和m2，
m1创建的gravitational field由
给出
其中G是universal gravitational constant，r是到m[=172=的距离]1和r̂是沿m1[连线的单位方向=205=] 和 m2.
这个场施加在m2上的引力由
给出
注意 - 重要的是，对于 any 距离为 any 的两个质点来说，这是正确的.1
引力相互作用的场性质允许我们使用 superposition 来计算由于多重相互作用而产生的净引力。考虑如果我们在前面的场景中添加另一个质量，m3，
那么质量 m2 上的引力就是每个 产生的场的引力之和=]其他质量，
和ri,j = rj,i。这适用于任何间隔的任何数量的质量。这也意味着，如果您更喜欢这种形式主义，则可以通过 vector sum 聚合质量集合创建的场。
现在考虑一下，如果我们有大量的质点，M，聚集在一起形成一个密度均匀的连续刚体。然后我们想计算单个空间不同点质量 m 上的引力，由于总质量 M:
然后我们可以考虑不同大小的质量区域（或体积），而不是考虑质点，并对这些区域（或体积）对质点的影响进行积分或求和。在二维情况下，引力的大小则为
其中σ是总质量的密度。2这相当于将每个微分质量引起的重力矢量场相加, σdxdy。这种等效性非常重要，因为它意味着对于质量分布外足够远的任何点质量，质量分布引起的引力几乎正好 与位于质量分布 center of mass 的质量 M 的点质量相同。3 4
这意味着，对于非常好的近似，在计算由于任何质量分布引起的引力场时，质量分布可以用质量分布中心的等效质量点质量代替。这适用于任何数量的空间上不同的质量分布，无论这些分布是否构成刚体。此外，这意味着您甚至可以将 组 分布聚合到系统质心的单个点质量中。5 只要参考点足够远.
但是，为了找到由于任何点处的质量分布引起的点质量上的引力，对于形状和分离不可知方式的任何质量分布，我们有通过对质量分布的每个部分的贡献求和来计算该点的引力场。6
回到问题
当然，对于任意多边形或多面体，解析解可能非常困难，因此使用求和要简单得多，算法方法也将类似地使用求和。
从算法上讲，这里最简单的方法实际上不是 geometric packing（使用 circles/spheres 或 squares/cubes）。使用打包并非不可能，但从数学上讲，这种方法存在重大挑战 - 最好采用依赖于更简单数学的方法。一种这样的方法是定义一个包含质量分布空间范围的网格，然后创建简单的（square/cubic 或 rectangular/cuboidic）多边形或以网格点为顶点的多面体。这将创建三种多边形或多面体：
- 不包含任何质量分布的那些
- 那些被质量分布完全填满的
- 那些被质量分布部分填充的
质心 - 方法 1
当参考点到质量分布的距离相对于分布的 angular 范围较大，并且没有质量分布对参考点的几何包围时（或任何几个发行版）。
然后您可以通过对每个多边形的贡献求和来找到分布的质心 R，
其中M是分布的总质量，ri是到第ith多边形几何中心的空间向量，mi 是密度乘以包含质量的多边形部分（即 1.00 对于完全填充的多边形和 0.00 用于完全空的多边形）。当您增加采样大小（网格点的数量）时，质心的近似值将接近解析解。一旦你有了质心，计算产生的引力场就很简单了：你只需在 R 点放置一个质量 M 的点质量并使用 equation from above.
为了演示，这里是 Python 中使用 shapely 库进行多边形操作的二维描述方法的实现：
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 接近 Planck length，但这些情况不是本问题的主题。
2在密度不均匀的情况下，这会变得更加复杂，并且在质量分布无法用符号描述的情况下，没有平凡的解析解。
3可能应该注意到，这实际上是 integral 正在做的事情；找到质心。
4对于在质量分布牛顿Shell Theorem内的点质量，或者必须使用场求和。
5在天文学中这被称为 barycenter，物体 总是 围绕系统的质心运行 - 而不是任何给定物体的质心。
6在某些情况下，使用牛顿 Shell Theorem 就足够了，但是这些情况与分布几何无关。
（为简单起见，在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
这是一个可行的解决方案吗？如果是这样，我将如何高效快速地打包圈子？如果不行，有没有更好的解决办法？
编辑：请注意我想如何找到物体在特定点的力，因此需要计算圆和观察者之间的角度（和向量求和）。它不同于计算所施加的总力。
背景
有些解释会有些跑题，但我认为有必要帮助澄清评论中提出的一些事情，因为其中大部分内容有点违反直觉。
这种引力相互作用的解释依赖于point masses的概念。假设你有两个质点，它们在一个孤立的系统中，彼此相隔一定距离，r1，质量为 m1和m2，
m1创建的gravitational field由
给出其中G是universal gravitational constant，r是到m[=172=的距离]1和r̂是沿m1[连线的单位方向=205=] 和 m2.
这个场施加在m2上的引力由
给出注意 - 重要的是，对于 any 距离为 any 的两个质点来说，这是正确的.1
引力相互作用的场性质允许我们使用 superposition 来计算由于多重相互作用而产生的净引力。考虑如果我们在前面的场景中添加另一个质量，m3，
那么质量 m2 上的引力就是每个 产生的场的引力之和=]其他质量，
和ri,j = rj,i。这适用于任何间隔的任何数量的质量。这也意味着，如果您更喜欢这种形式主义，则可以通过 vector sum 聚合质量集合创建的场。
现在考虑一下，如果我们有大量的质点，M，聚集在一起形成一个密度均匀的连续刚体。然后我们想计算单个空间不同点质量 m 上的引力，由于总质量 M:
然后我们可以考虑不同大小的质量区域（或体积），而不是考虑质点，并对这些区域（或体积）对质点的影响进行积分或求和。在二维情况下，引力的大小则为
其中σ是总质量的密度。2这相当于将每个微分质量引起的重力矢量场相加, σdxdy。这种等效性非常重要，因为它意味着对于质量分布外足够远的任何点质量，质量分布引起的引力几乎正好 与位于质量分布 center of mass 的质量 M 的点质量相同。3 4
这意味着，对于非常好的近似，在计算由于任何质量分布引起的引力场时，质量分布可以用质量分布中心的等效质量点质量代替。这适用于任何数量的空间上不同的质量分布，无论这些分布是否构成刚体。此外，这意味着您甚至可以将 组 分布聚合到系统质心的单个点质量中。5 只要参考点足够远.
但是，为了找到由于任何点处的质量分布引起的点质量上的引力，对于形状和分离不可知方式的任何质量分布，我们有通过对质量分布的每个部分的贡献求和来计算该点的引力场。6
回到问题
当然，对于任意多边形或多面体，解析解可能非常困难，因此使用求和要简单得多，算法方法也将类似地使用求和。
从算法上讲，这里最简单的方法实际上不是 geometric packing（使用 circles/spheres 或 squares/cubes）。使用打包并非不可能，但从数学上讲，这种方法存在重大挑战 - 最好采用依赖于更简单数学的方法。一种这样的方法是定义一个包含质量分布空间范围的网格，然后创建简单的（square/cubic 或 rectangular/cuboidic）多边形或以网格点为顶点的多面体。这将创建三种多边形或多面体：
- 不包含任何质量分布的那些
- 那些被质量分布完全填满的
- 那些被质量分布部分填充的
质心 - 方法 1
当参考点到质量分布的距离相对于分布的 angular 范围较大，并且没有质量分布对参考点的几何包围时（或任何几个发行版）。
然后您可以通过对每个多边形的贡献求和来找到分布的质心 R，
其中M是分布的总质量，ri是到第ith多边形几何中心的空间向量，mi 是密度乘以包含质量的多边形部分（即 1.00 对于完全填充的多边形和 0.00 用于完全空的多边形）。当您增加采样大小（网格点的数量）时，质心的近似值将接近解析解。一旦你有了质心，计算产生的引力场就很简单了：你只需在 R 点放置一个质量 M 的点质量并使用 equation from above.
为了演示，这里是 Python 中使用 shapely 库进行多边形操作的二维描述方法的实现：
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 接近 Planck length，但这些情况不是本问题的主题。
2在密度不均匀的情况下，这会变得更加复杂，并且在质量分布无法用符号描述的情况下，没有平凡的解析解。
3可能应该注意到，这实际上是 integral 正在做的事情；找到质心。
4对于在质量分布牛顿Shell Theorem内的点质量，或者必须使用场求和。
5在天文学中这被称为 barycenter，物体 总是 围绕系统的质心运行 - 而不是任何给定物体的质心。
6在某些情况下，使用牛顿 Shell Theorem 就足够了，但是这些情况与分布几何无关。