球相交模板

模板"/>

球相交模板

#define PI acos(-1.0)
#define ll long long int
const int maxn = 1e2 + 5;
using namespace std;
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;}
};
double dist(point p1, point p2) {       //返回平面上两点距离return sqrt((p1 - p2)*(p1 - p2));
}
struct sphere {   //球double r, area, tj; // 半径, 表面积, 体积等.point centre;
};
sphere s,a[maxn];
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;//剩余部分表面积
}
void solve() {int n; cin >> n;double x, y, z, r;for(int i = 1; i <= n; i++) {scanf("%lf%lf%lf%lf",&x,&y,&z,&a[i].r); //其他球a[i].centre = {x,y,z};}scanf("%lf%lf%lf%lf",&x,&y,&z,&r);s.r = r;s.centre = {x,y,z}; //中心球double ans = 0, v = 0;for(int i = 1; i <= n; i++) {double ss, dis = dist(s.centre, a[i].centre);if(dis >= s.r + a[i].r)continue;  //在外部if(dis + min(s.r, a[i].r) <= max(s.r, a[i].r))  //在内部{ans += 4.0 / 3.0 * PI * min(s.r,a[i].r) * min(s.r,a[i].r) * min(s.r,a[i].r);continue;}SphereInterVS(s, a[i], v, ss); //相交部分ans += v;}printf("%.8f\n" ,ans);
}