POJ3944 Spherical Mirrors

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140

题意：一束光线射向三维空间中的很多个镜面球体，遇到球面发生反射，经过有 限次后光线不会再与球面相遇问最后光线照射到的球面上的坐标 思路：遇到球面发生反射，求出起点关于法线对称点，光线的起点改为照射到的 球面的坐标，射向方向为新起点照向对称点最先遇到的球面上的坐标即为反射点 需先求一次关于时间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; }