求两个球体相交部分体积，公式及推导，C++实现

公式

两球相交分为两种情况：
第一种 大球完全包围了小球
第二种 两球相交但无完全包含关系

第一种情况的答案是显而易见的，直接求出小球体积即可
下面给出第二种情况的求解方法

$$V=\pi h_1^2(R-\frac{1}{3}{h}_1)+\pi h_2^2(r-\frac{1}{3}{h}_2)$$
图中角$\alpha$与角$\beta$位置标反了（多谢评论区小伙伴的提醒）

推导

我们把相交部分面积分为两部分，一部分是大球刨出去的体积，一部分是小球抛出去的体积，并且这两部分面积计算公式相同。
接下来我们证明其中大球刨去体积的公式。
先做一些准备工作
$$\cos \alpha =\frac{R^2+d^2-r^2}{2Rd}$$
$$h_1=R-Rcos\alpha$$
为了应用黎曼积分求体积，首先要得到大球被抛去部分的圆面面积公式。
由于刨去部分圆面的半径随着圆面的圆心从E点到F点的距离变小而变小，为了研究这两者间的关系，我们设圆面的圆心到F点的距离为 $x$ 。
得到圆面半径的表达式为 $\sqrt{R^2-(Rcos\alpha +x{)}^2})$ ，从而圆面面积表达式为 $\pi (R^2-(Rcos\alpha +x{)}^2)$
于是得到大球抛去体积 $V_1$ 的积分公式为
$$V_1={\int }_0^{h_1}\pi (R^2-(Rcos\alpha +x{)}^2)dx$$
经过一系列牛顿-莱布尼茨公式之后化简得到
$$V_1=\pi (-\frac{1}{3}{h}_1^2-R{h}_1^{2}cos\alpha +R^2{h}_{1}si{n}^2\alpha )$$
再把 $cos\alpha$ 代入化简后得到
$$V_1=\pi h_1^2(R-\frac{1}{3}{h}_1)$$
同理可得小球被刨去的体积
$$V_2=\pi h_2^2(r-\frac{1}{3}{h}_2)$$
总体积
$$V=V_1+V_2=\pi h_1^2(R-\frac{1}{3}{h}_1)+\pi h_2^2(r-\frac{1}{3}{h}_2)$$

代码

话不多说，上代码

#include<iostream>
#include<cmath>
#define Pi 3.14

using namespace std;

class Sphere	//属性分别为:(x,y,z)三个坐标以及半径
{
private:
	double m_x;
	double m_y;
	double m_z;
	double m_r;
public:
	Sphere():m_x(0),m_y(0), m_z(0), m_r(0)
	{
	}
	Sphere(double x,double y,double z,double r) :m_x(x), m_y(y), m_z(z), m_r(r)
	{
	}

	double getx();
	double gety();
	double getz();
	double getr();
	void modi_x(double x);
	void modi_y(double y);
	void modi_z(double z);
	void modi_r(double r);
	bool intersect(Sphere& o);
	double vol_ints(Sphere& o);
	double distance(Sphere& o);
};
//获得属性
double Sphere::getx()
{
	return m_x;
}
double Sphere::gety()
{
	return m_y;
}
double Sphere::getz()
{
	return m_z;
}
double Sphere::getr()
{
	return m_r;
}
//修改属性
void Sphere::modi_x(double x)
{
	m_x = x;
}
void Sphere::modi_y(double y)
{
	m_y = y;
}
void Sphere::modi_z(double z)
{
	m_z = z;
}
void Sphere::modi_r(double r)
{
	m_r = r;
}
//判断两个球是否相交
bool Sphere::intersect(Sphere& o)
{
	if (m_x == o.m_x && m_y == o.m_y && m_z == o.m_z)
		return true;
	else
	{
		double dis = distance(o);
		if ((m_r + o.m_r) * (m_r + o.m_r) > dis)
			return true;
		else
			return false;
	}
}
//计算两球圆心之间的距离
double Sphere::distance(Sphere& o)
{
	return (m_x - o.m_x) * (m_x - o.m_x) + (m_y - o.m_y) * (m_y - o.m_y) + (m_z - o.m_z) * (m_z - o.m_z);
}
//计算两球相交部分体积
double Sphere::vol_ints(Sphere& o)
{
	if (!intersect(o))
		return 0;
	else
	{
		double dis = distance(o);
		double max = m_r > o.m_r ? m_r : o.m_r;
		double min = m_r > o.m_r ? o.m_r : m_r;
		if (sqrt(dis) + min <= max)
			return 4 / 3 * Pi * pow(min, 3);
		else
		{
			double cos_a = (max * max + dis * dis - min * min) / (2 * max * dis);
			double cos_b = (min * min + dis * dis - max * max) / (2 * min * dis);
			double h_a = max * (1 - cos_a);
			double h_b = min * (1 - cos_b);
			return (Pi / 3) * (3 * max - h_a) * h_a * h_a + (Pi / 3) * (3 * min - h_b) * h_b * h_b;
		}
	}
}

第一次写文章，发现错误的话请各位多多指教