POJ3944 Spherical Mirrors-优快云博客

本文链接：https://blog.youkuaiyun.com/hnust_Derker/article/details/52433470

本文介绍了一种在三维空间中模拟光线与多个镜面球体交互的方法。通过求解光线与球体表面的交点，并利用向量运算来确定光线反射路径，最终找到光线最后一次接触的球面坐标。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

题意：一束光线射向三维空间中的很多个镜面球体，遇到球面发生反射，经过有

限次后光线不会再与球面相遇问最后光线照射到的球面上的坐标

思路：遇到球面发生反射，求出起点关于法线对称点，光线的起点改为照射到的

球面的坐标，射向方向为新起点照向对称点最先遇到的球面上的坐标即为反射点

需先求一次关于时间t的二次方程，t < 0 舍去

#include<cstdio>

#include<cstring>

#include<cmath>

#include<vector>

#include<set>

#include<map>

#include<algorithm>

const int maxn = 1e2;

const double pi = acos(-1.0);

const double eps = 1e-5;

using namespace std;

//加法误差

double add(double a, double b) {

    if(fabs(a + b) < eps) return 0;

    return a + b;

}

//向量运算

struct P {

    double x, y, z;

    P() {}

    P(double x, double y, double z) : x(x), y(y), z(z) {}

    P operator + (P p) {

        return P(add(x, p.x), add(y, p.y), add(z, p.z));

    }

    P operator - (P p) {

        return P(add(x, -p.x), add(y, -p.y), add(z, -p.z));

    }

    P operator * (double d) {

        return P(x * d, y * d, z * d);

    }

    //点乘

    double dot(P p) {

        return add(add(x * p.x, y * p.y), z * p.z);

    }

    //面积

    double area(P p) {

        double dx =  add(y * p.z, -z * p.y);

        double dy = -add(x * p.z, -z * p.x);

        double dz =  add(x * p.y, -y * p.x);

        double s = add(add(dx * dx, dy * dy), dz * dz);

        return sqrt(s);

    }

    //向量模长

    double lenth() {

        return sqrt(add(add(x * x, y * y), z * z));

    }

    //单位向量

    P unit() {

        return P(x, y, z) * (1 / lenth());

    }

    //向量p1-p在向量p1-p2上的投影模长

    double prj_len(P p1, P p2) {

        P p(x, y, z);

        return (p - p1).dot((p2 - p1).unit());

    }

    //在p1-p2上的投影(垂足)

    P prj(P p1, P p2) {

        P p(x, y, z);

        return p1 + (p2 - p1).unit() * p.prj_len(p1, p2);

    }

    //关于p1-p2的对称点

    P sym(P p1, P p2) {

        P p(x, y, z);

        P e = p.prj(p1, p2);

        return e * 2 - p;

    }

};

struct dis {

    P p;

    double d;

    int ind;

    dis() {}

    dis(P p, double d, int ind) : p(p), d(d), ind(ind) {}

};

bool operator < (dis a, dis b) {

    return a.d < b.d;

}

int main() {

    int n, num;

    double r[maxn];

    P st, a[maxn], s;

    dis df[maxn];

    while(scanf("%d", &n) && n) {

        int m = -1;

        scanf("%lf %lf %lf", &s.x, &s.y, &s.z);

        st.x = st.y = st.z = 0;

        for(int i = 0; i < n; i++) {

            scanf("%lf %lf %lf %lf", &a[i].x, &a[i].y, &a[i].z, &r[i]);

        }

        while(1) {

            num = 0;

            P next = st + s, toy;

            for(int i = 0; i < n; i++) {

                if(i == m) continue;

                P e = a[i].prj(st, next);

                P ss = e - a[i];

                double d = ss.lenth();

                if(d > r[i] || !add(r[i], -d)) continue; //相切或相离

                double dx = s.x, dy = s.y, dz = s.z;

                double Dx = add(st.x, -a[i].x), Dy = add(st.y, -a[i].y), Dz = add(st.z, -a[i].z);

                // A*t^2 + B*t + C = 0

                //printf("dxyz = %lf %lf %lf\n", dx, dy, dz);

                //printf("Dxyz = %lf %lf %lf\n", Dx, Dy, Dz);

                double A = add(add(dx * dx, dy * dy), dz * dz);

                double B = 2 * add(add(dx * Dx, dy * Dy), dz * Dz);

                double C = add(add(Dx * Dx, Dy * Dy), add(Dz * Dz, -r[i] * r[i]));

                double lf = add(B * B,  -4 * A * C);

                //printf("Number.%d dt = %lf\n\n", i, sqrt(lf));

                if(lf <= 0) continue;

                double dd = sqrt(lf);

                double t1 = add(-B, dd) / 2 / A, t2 = add(-B, -dd) / 2 / A;

                P p1 = P(add(t1 * dx, st.x), add(t1 * dy, st.y), add(t1 * dz, st.z));

                double d1 = (p1 - st).lenth();

                P p2 = P(add(t2 * dx, st.x), add(t2 * dy, st.y), add(t2 * dz, st.z));

                double d2 = (p2 - st).lenth();

                if(s.dot(p1 - st) > 0) df[num++] = dis(p1, d1, i);

                if(s.dot(p2 - st) > 0) df[num++] = dis(p2, d2, i);

            }

            if(!num) break;

            sort(df, df + num);

            m = df[0].ind;

            P e = st.sym(a[df[0].ind], df[0].p);

            st = df[0].p;

            s = e - st;

            P p2 = e;

        }

        printf("%.4lf %.4lf %.4lf\n", st.x, st.y, st.z);

    }

    return 0;

}