首先,定义一些东西
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;//剩余部分表面积
}