两球体积交

#include<bits/stdc++.h>
using namespace std;
#define LL long long
const int N = 2e6;
const long long  inf = 1e9;
const long long mod = 998244353;
const double PI = acos(-1);
struct Point
{
    long  double x, y, z;
}p[6];
struct Ball
{
    Point o;
    long double r;
} o[4];
int n;
double r1, r2;
long double a, b, c, d, w;
int k, k2;
long double getdis(Point a, Point b)
{
    return sqrt((b.x - a.x) * (b.x - a.x) + (b.y - a.y) * (b.y - a.y) + (b.z - a.z) * (b.z - a.z));
}
long double get_Ball_area(struct Ball a, struct Ball b)
{
    long double ans = 0;
    long double dis = getdis(a.o, b.o);
    long double  r1 = a.r, r2 = b.r;
    Point  o1 = a.o, o2 = b.o;
    if (dis >= r1 + r2) return 0;
    if (dis + r1 <= r2) return 1.0 * 4 / 3 * PI * r1 * r1 * r1;
    if (dis + r2 <= r1) return 1.0 * 4 / 3 * PI * r2 * r2 * r2;
    long double cosa = (r1 * r1 + dis * dis - r2 * r2) / (2.0 * r1 * dis);
    long double d1 = r1 - r1 * cosa;
    ans += PI * (r1 * d1 * d1 - 1 / 3.0 * d1 * d1 * d1);//求冠面积公式
    long double cosb = (r2 * r2 + dis * dis - r1 * r1) / (2.0 * r2 * dis);
    long double d2 = r2 - r2 * cosb;
    ans += PI * (r2 * d2 * d2 - 1 / 3.0 * d2 * d2 * d2);
    return ans;
}
int main()
{
    int T;
    cin >> T;
    while (T--)
    {
        for (int i = 1; i <= 4; i++)
            cin >> p[i].x >> p[i].y >> p[i].z;
           cin >> k >> k2;
        a = k * k - 1;
        b = k * k * p[2].x - p[1].x;
        c = k * k * p[2].y - p[1].y;
        d = k * k * p[2].z - p[1].z;
        w = k * k*(p[2].x * p[2].x + p[2].y * p[2].y + p[2].z * p[2].z) - p[1].x * p[1].x - p[1].y * p[1].y - p[1].z * p[1].z;
        o[1].o = { b / a,c / a,d / a };
        o[1].r = sqrt((b * b + c * c + d * d - w * a) / (a * a));
        a = k2 * k2 - 1;
        b = k2 * k2 * p[4].x - p[3].x;
        c = k2 * k2 * p[4].y - p[3].y;
        d = k2 * k2 * p[4].z - p[3].z;
        w = k2 * k2*(p[4].x * p[4].x + p[4].y * p[4].y + p[4].z * p[4].z) - p[3].x * p[3].x - p[3].y * p[3].y - p[3].z * p[3].z;
        o[2].o = { b / a,c / a,d / a };
        o[2].r = sqrt((b * b + c * c + d * d - w * a) / (a * a));
        long double ans = get_Ball_area(o[1], o[2]);
        cout<<ans<<endl;
    }
    return 0;
}