样例输入:
1
1 0 0
3 0 0
2 0 0
4 0 0
3 3
样例输出:
0.262
简化题意:给你三维空间中的四个点A,B,C,D,然后给你一个K1和K2,一个点P1到A的距离不小于P1到B的距离的K1倍,一个点P2到C的距离不小于P2到D的距离的K2倍,求P1点所形成的轨迹与P2点所形成的轨迹相交的体积。
以P1点的轨迹为例:设P1(x,y,z),则有
(x-xA)*(x-xA)+(y-yA)*(y-yA)+(z-zA)*(z-zA)>=k*k*((x-xB)*(x-xB)+(y-yB)*(y-yB)+(z-zB)*(z-zB))
经过对这个等式化简不难发现这是一个球,化成标准形式很容易得到圆心坐标和半径,然后我们只需要待入计算几何的模板即可求得两球相交体积
下面是代码:
#include<iostream>
#include<cstdio>
#include<string>
#include<ctime>
#include<cmath>
#include<cstring>
#include<algorithm>
#include<stack>
#include<climits>
#include<queue>
#include<map>
#include<set>
#include<sstream>
#include<cassert>
#include<bitset>
#include<list>
#include<unordered_map>
#define double long double
#define lowbit(x) (x&-x)
using namespace std;
typedef long long LL;
typedef unsigned long long ull;
template<typename T>
inline void read(T &x)
{
T f=1;x=0;
char ch=getchar();
while(0==isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
while(0!=isdigit(ch)) x=(x<<1)+(x<<3)+ch-'0',ch=getchar();
x*=f;
}
template<typename T>
inline void write(T x)
{
if(x<0){x=~(x-1);putchar('-');}
if(x>9)write(x/10);
putchar(x%10+'0');
}
const int inf=0x3f3f3f3f;
const int N=1e6+100;
const double pi=acos(-1);
double pow2(double x){return x*x;}
double pow3(double x){return x*x*x;}
double dis(double x1,double y1,double z1,double x2,double y2,double z2)
{
return pow2(x1-x2)+pow2(y1-y2)+pow2(z1-z2);
}
double cos(double a,double b,double c){return (b*b+c*c-a*a)/(2*b*c);}
double cap(double r,double h){return pi*(r*3-h)*h*h/3;}
//2球体积交
double sphere_intersect(double x1,double y1,double z1,double r1,double x2,double y2,double z2,double r2)
{
double d=dis(x1,y1,z1,x2,y2,z2);
//相离
if(d>=pow2(r1+r2))return 0;
//包含
if(d<=pow2(r1-r2))return pow3(min(r1,r2))*4*pi/3;
//相交
double h1=r1-r1*cos(r2,r1,sqrt(d)),h2=r2-r2*cos(r1,r2,sqrt(d));
return cap(r1,h1)+cap(r2,h2);
}
//2球体积并
double sphere_union(double x1,double y1,double z1,double r1,double x2,double y2,double z2,double r2)
{
double d=dis(x1,y1,z1,x2,y2,z2);
//相离
if(d>=pow2(r1+r2))return (pow3(r1)+pow3(r2))*4*pi/3;
//包含
if(d<=pow2(r1-r2))return pow3(max(r1,r2))*4*pi/3;
//相交
double h1=r1+r1*cos(r2,r1,sqrt(d)),h2=r2+r2*cos(r1,r2,sqrt(d));
return cap(r1,h1)+cap(r2,h2);
}
void solve()
{
double x[5],y[5],z[5];
double x1, y1, z1, r1, x2, y2, z2, r2;
for(int i=1;i<=4;i++) cin>>x[i]>>y[i]>>z[i];
double k1,k2;
cin>>k1>>k2;
x1 = ((k1*k1*x[2]-x[1])/(k1*k1-1));
y1 = ((k1*k1*y[2]-y[1])/(k1*k1-1));
z1 = ((k1*k1*z[2]-z[1])/(k1*k1-1));
r1 = (x1*x1+y1*y1+z1*z1-(k1*k1*(x[2]*x[2]+y[2]*y[2]+z[2]*z[2]))/(k1*k1-1)+(x[1]*x[1]+y[1]*y[1]+z[1]*z[1])/(k1*k1-1));
r1 = sqrt(r1);
k1 = k2;
x[1] = x[3];y[1] = y[3];z[1] = z[3];
x[2] = x[4];z[2] = z[4];y[2] = y[4];
x2 = ((k1*k1*x[2]-x[1])/(k1*k1-1));
y2 = ((k1*k1*y[2]-y[1])/(k1*k1-1));
z2 = ((k1*k1*z[2]-z[1])/(k1*k1-1));
r2 = (x2*x2+y2*y2+z2*z2-(k1*k1*(x[2]*x[2]+y[2]*y[2]+z[2]*z[2])/(k1*k1-1))+(x[1]*x[1]+y[1]*y[1]+z[1]*z[1])/(k1*k1-1));
r2 = sqrt(r2);
double v = sphere_intersect(x1, y1, z1, r1, x2, y2, z2, r2);
printf("%Lf\n", v);
}
int main(){
int _;
scanf("%d",&_);
while(_--) solve();
return 0;
}