CCPC-Wannafly Winter Camp Day2 H

 

题目描述

 

在一片小行星带里有 nn 颗小行星，它们在万有引力的作用下绕着一颗行星旋转。在这一刻时，它们之间不存在碰撞的情况。一位清洁工奉命前来清理这颗行星，Ta 会动用某种先进技术使这颗行星顷刻间从宇宙中消失，任何距离这颗行星的中心在一定范围内的事物都会在一瞬间被清除。假设这些天体都是完整的球体，你能计算出清除的区域里有多少体积的事物原本属于这些小行星吗？

注意，这些天体在此刻满足两两不存在交集的条件。

下图是对样例的解释，其中清理区域是标记为红色的球体内部，而小行星则被依次标记为橙色、蓝色和绿色。

输入描述

 

输入包含多组测试数据。第一行包含一个整数 TT，表示测试数据的组数。随后的内容是各组测试数据。对于每组测试数据：

第一行包含一个整数 nn。

接下来的 nn 行里，每行包含四个整数 x, y, zx,y,z 和 rr，表示有一颗中心位于 (x, y, z)(x,y,z)、半径为 rr 的小行星。

最后一行包含四个整数 x', y', z'x′,y′,z′ 和 r'r′，表示行星的中心位于 (x', y', z')(x′,y′,z′)，而清洁工的清理半径为 r'r′（一个大于该行星半径的值）。

• 1 \leq T \leq 60001≤T≤6000
• 1 \leq n \leq 1001≤n≤100
• -10^3 \leq x, y, z, x', y', z' \leq 10^3−103≤x,y,z,x′,y′,z′≤103
• 1 \leq r, r' \leq 10^31≤r,r′≤103

输出描述

 

对于每组测试数据，输出一行Case #x: y，其中x是测试数据的编号（从 11 开始编号），y是这组数据的答案，要求相对误差或绝对误差不超过 10^{-6}10−6。

严格来讲，如果你的答案是 aa，而标准答案是 bb，那么当 \frac{|a - b|}{\max{1, |b|}} \leq 10^{-6}max1,∣b∣∣a−b∣​≤10−6 时你的答案会被认为是正确的。

样例输入 1 

1
3
5 5 5 2
-6 -7 6 1
6 -5 0 3
1 -1 0 10

样例输出 1 

Case #1: 142.76246874761383764962

 

两个球的体积交，模版套一下。 

#include<bits/stdc++.h>
typedef long long int ll;

const int maxn=105;
const double PI=acos(-1.0);
double check(int x,int y,int z,int x1,int y1,int z1)
{
	return (double)sqrt((x-x1)*(x-x1)+(y-y1)*(y-y1)+(z-z1)*(z-z1));
}
typedef struct point {
    double x,y,z,r;
    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));
}
double SphereInterVS(point a, point b) {
	double v;
    double d =dist(a,b);//球心距
    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;//相交部分体积
    return v;
}

struct node
{
	int x,y,z,r;
}a[maxn];

double get_v(int r)
{
	return (double)4*PI*r*r*r/3.0;
}
int main()
{
	int t;scanf("%d",&t);
	for(int cc=1;cc<=t;cc++)
	{
		int n;scanf("%d",&n);
		for(int i=1;i<=n;i++)
		{
			scanf("%d%d%d%d",&a[i].x,&a[i].y,&a[i].z,&a[i].r);
		}
		
		int X,Y,Z,R;scanf("%d%d%d%d",&X,&Y,&Z,&R);
		double ans=0;
		for(int i=1;i<=n;i++)
		{
			if(check(a[i].x,a[i].y,a[i].z,X,Y,Z)<=R-a[i].r)  ans+=get_v(a[i].r);
			else if((check(a[i].x,a[i].y,a[i].z,X,Y,Z)>R-a[i].r)&&(check(a[i].x,a[i].y,a[i].z,X,Y,Z)<R+a[i].r))
			{
			    point c;point b;
			    c.x=(double)a[i].x,c.y=(double)a[i].y,c.z=(double)a[i].z,c.r=(double)(a[i].r);
				b.x=(double)X,b.y=(double)Y,b.z=(double)Z,b.r=(double)R;
				ans+=SphereInterVS(b,c);
			}
		}
		printf("Case #%d: %lf\n",cc,ans);
	}
	return 0;
}