兩圓相交到兩球相交

首先，定義一些東西

const double PI = acos(-1.0);
typedef struct point {
	double x,y;
	point() {
	}
	point(double a, double b) {
		x = a;
		y = b;
	}
	point operator -(const point &b)const {		//返回減去后的新點
		return point(x - b.x, y - b.y);
	}
	point operator +(const point &b)const {		//返回加上后的新點
		return point(x + b.x, y + b.y);
	}
	point operator *(const double &k)const {	//返回相乘后的新點
		return point(x * k, y * k);
	}
	point operator /(const double &k)const {	//返回相除后的新點
		return point(x / k, y / k);
	}
	double operator ^(const point &b)const {	//叉乘
		return x*b.y - y*b.x;
	}
	double operator *(const point &b)const {	//點乘
		return x*b.x + y*b.y;
	}
}point;
typedef struct circle {//圓
	double r;
	point centre;
}circle;

double dist(point p1, point p2) {		//返回平面上兩點距離
	return sqrt((p1 - p2)*(p1 - p2));
}

• 兩圓相交
  兩圓關系，可以根據圓心距離和半徑的關系來判斷，現在只考慮相交的情況，即圓心距\(L\)在兩圓半徑之和\(|r_1+r_2|\)及兩圓半徑之差\(|r_1-r_2|\)之間。
  如上圖所示，已知\(r_1,r_2,L\)那就可以得到很多東西。
  根據勾股定理，可以得到
  \(r_1^2-h_1^2=l_1^2\)
  \(r_2^2-h_1^2=l_2^2\)
  \(L=l_1+l_2\)
  聯立推出
  \(h_1=\sqrt{r_1^2-(\frac{L^2+r_1^2-r_2^2}{2L})^2}\)
  \(h_1=\sqrt{r_2^2-(\frac{L^2+r_2^2-r_1^2}{2L})^2}\)
  \(h_1=h_2\)
• 有了這個就可以算兩圓相交的弧長了，但是需要先算出角度

用余弦公式可以算出\(angle_a\)和\(angle_b\)
\(cos(angle_a)=\frac{r_1^2+L^2-r_2^2}{2Lr_1}\)
\(cos(angle_b)=\frac{r_2^2+L^2-r_1^2}{2Lr_2}\)
但是實際相交的弧長所對應的圓心角是上述所求角的兩倍，所以乘以2。
最后利用弧長公式即可計算兩圓相交部分的弧長。
弧長=\(Pi\)*對應圓心角(弧度制)

void CircleInterLen(circle a, circle b, double &la, double &lb) {
	double d = dist(a.centre, b.centre);//圓心距
	double t = (d*d + a.r*a.r - b.r*b.r) / (2.0 * d);//
	double h = sqrt((a.r*a.r) - (t*t)) * 2;//h1=h2
	double angle_a = 2 * acos((a.r*a.r + d*d - b.r*b.r) / (2.0 * a.r*d));  
	//余弦公式計算r1對應圓心角，弧度
	double angle_b = 2 * acos((b.r*b.r + d*d - a.r*a.r) / (2.0 * b.r*d));  
	//余弦公式計算r2對應圓心角，弧度
	la = angle_a*a.r;//r1所對應的相交弧長
	lb = angle_b*b.r;//r2所對應的相交弧長
	//double rest_la = 2.0 * PI * a.r - la;//r1圓剩余部分弧長
	//double rest_lb = 2.0 * PI * b.r - lb;//r2圓剩余部分弧長
}

剩下部分的弧長只需要用原來的弧長相減就行。

• 再看相交部分的面積:

扇形面積公式:\(S=\frac{lr}{2}\)(l為弧長，r為半徑)=\(\frac{nr^2}{2}\)(n為圓心角，弧度制)
相交部分面積包含扇形減去三角形的面積

void CircleInterArea(circle a, circle b,double &s1,double &s2) {
	double d = dist(a.centre, b.centre);//圓心距
	double t = (d*d + a.r*a.r - b.r*b.r) / (2.0 * d);//
	double h = sqrt((a.r*a.r) - (t*t)) * 2;//h1=h2
	double angle_a = 2 * acos((a.r*a.r + d*d - b.r*b.r) / (2.0 * a.r*d));  
	//余弦公式計算r1對應圓心角，弧度
	double angle_b = 2 * acos((b.r*b.r + d*d - a.r*a.r) / (2.0 * b.r*d));  
	//余弦公式計算r2對應圓心角，弧度
	double la = angle_a*a.r;//r1所對應的相交弧長
	double lb = angle_b*b.r;//r2所對應的相交弧長
	s1 = la*a.r / 2.0 - a.r*a.r*sin(angle_a) / 2.0;	//相交部分r1圓的面積
	s2 = lb*b.r / 2.0 - b.r*b.r*sin(angle_b) / 2.0;	//相交部分r2圓的面積
	double rest_s1 = PI*a.r*a.r - s1 - s2;//r1圓剩余部分面積，不含相交部分面積
	double rest_s2 = PI*b.r*b.r - s1 - s2;//r2圓剩余部分面積，不含相交部分面積
}

下面考慮兩球相交：

相交部分如下:

但實際上，如果將其投影至平面，還是剛才的樣子

從上可以知道，相交部分體積是兩個球缺的和。
球冠面積\(S=2πrh\)
球缺的體積公式為\(V=πh^2(r-\frac{h}{3})\)
\(h\)球冠高
\(r\)球半徑

兩球心距離=\(L\) 半徑分別為\(r_1 r_2\)
\(|r_1- r_2|< L< |r_1+ r_2|\)
兩球相交的截面為平面，相交線為圓半徑為\(r_3\)
截面到球心的距離分別為\(l_1 l_2\)
\(l_1+l_2=L\)
\(L\)直線過相交圓心並垂直相交圓直徑

\(r_1^2=r_3^2+l_1^2\)
\(r_2^2=r_3^2+l_2^2\)

\(r_1^2-r_2^2=l_1^2-l_2^2\)
\(r_1^2-r_2^2=(l_1+l_2)(l_1-l_2)\)
\(r_1^2-r2^2=(2l_1-L)L\)
\(l_1=[(r_1^2-r_2^2)/L+L]/2\)
\(l_2=L-l_1\)

\(x_1=r_1-l_1\)
\(x_2=r_2-l_2\)

typedef struct point {
	double x,y,z;
	point() {

	}
	point(double a, double b,double c) {
		x = a;
		y = b;
		z = c;
	}
	point operator -(const point &b)const {		//返回減去后的新點
		return point(x - b.x, y - b.y,z-b.z);
	}
	point operator +(const point &b)const {		//返回加上后的新點
		return point(x + b.x, y + b.y,z+b.z);
	}
	//數乘計算
	point operator *(const double &k)const {	//返回相乘后的新點
		return point(x * k, y * k,z*k);
	}
	point operator /(const double &k)const {	//返回相除后的新點
		return point(x / k, y / k,z/k);
	}
	double operator *(const point &b)const {	//點乘
		return x*b.x + y*b.y+z*b.z;
	}
}point;
double dist(point p1, point p2) {		//返回平面上兩點距離
	return sqrt((p1 - p2)*(p1 - p2));
}
typedef struct sphere {//球
	double r;
	point centre;
}sphere;
void SphereInterVS(sphere a, sphere b,double &v,double &s) {
	double d = dist(a.centre, b.centre);//球心距
	double t = (d*d + a.r*a.r - b.r*b.r) / (2.0 * d);//
	double h = sqrt((a.r*a.r) - (t*t)) * 2;//h1=h2，球冠的高
	double angle_a = 2 * acos((a.r*a.r + d*d - b.r*b.r) / (2.0 * a.r*d));  //余弦公式計算r1對應圓心角，弧度
	double angle_b = 2 * acos((b.r*b.r + d*d - a.r*a.r) / (2.0 * b.r*d));  //余弦公式計算r2對應圓心角，弧度
	double l1 = ((a.r*a.r - b.r*b.r) / d + d) / 2;
	double l2 = d - l1;
	double x1 = a.r - l1, x2 = b.r - l2;//分別為兩個球缺的高度
	double v1 = PI*x1*x1*(a.r - x1 / 3);//相交部分r1圓所對應的球缺部分體積
	double v2 = PI*x2*x2*(b.r - x2 / 3);//相交部分r2圓所對應的球缺部分體積
	 v = v1 + v2;//相交部分體積
	double s1 = PI*a.r*x1;  //r1對應球冠表面積
	double s2 = PI*a.r*x2;	//r2對應球冠表面積
	 s = 4 * PI*(a.r*a.r + b.r*b.r) - s1 - s2;//剩余部分表面積
}