SGU110 计算几何 Computational Geometry

题意：三维空间N个球体，给一个光束的方向，求它在这些球体间的反射情况(球面反射)。

Problem: N spheres are in a 3D-space. Give you a direction of a light and write a program to output the reflections(the light rays reflect from the surface of the spheres according the ordinary law).

题解：如下图所示，每次反射分为两步。第一步利用图中关系求出tt，进而求出反射点。第二步利用投影定理求出反射光线的方向向量。

Solutions:  As shown in the figure, each reflection can be considered as two steps. Firstly, calculate the value of tt according to the figure, then make the reflect point. Secondly,  calculate the direction vector of the reflected ray by the projection theorem.

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <algorithm>
using namespace std;

#define P(x) ((x)*(x))

const int N=50+10;
const double oo=(1LL<<60);
const double eps=1e-8;

int n,cnt;
double r[N],tt;
struct Point {
    Point() {}
    Point(double x,double y,double z): x(x),y(y),z(z) {}
    double x,y,z;
}c[N],s,e,dir;

double operator *(Point A,Point B) {return A.x*B.x+A.y*B.y+A.z*B.z;}
Point operator *(double B,Point A) {return Point(A.x*B,A.y*B,A.z*B);}
Point operator +(Point A,Point B) {return Point(A.x+B.x,A.y+B.y,A.z+B.z);}
Point operator -(Point A,Point B) {return Point(A.x-B.x,A.y-B.y,A.z-B.z);}

inline int next() {
    int w = 0;
    double a,b,cc,delta,gen;
    tt = oo;
    for(int i=1;i<=n;i++) {
        a = dir*dir;
        b = 2.0*(s-c[i])*dir;
        cc = P(s-c[i]) - P(r[i]);
        if ((delta = P(b)-4*a*cc)>-eps) {
           if (fabs(delta)<eps) delta = 0;
           else delta = sqrt(delta);
           
           gen = (-b+delta)/a*0.5;
           if (gen>eps && gen<tt) tt=gen,w=i;
           
           gen = (-b-delta)/a*0.5;
           if (gen>eps && gen<tt) tt=gen,w=i;
        }
    }
    return w;
}
int main() {
    scanf("%d",&n);
    cnt = 0;
    for (int i=1;i<=n;i++)
        scanf("%lf%lf%lf%lf",&c[i].x,&c[i].y,&c[i].z,&r[i]);
    scanf("%lf%lf%lf%lf%lf%lf",&s.x,&s.y,&s.z,&e.x,&e.y,&e.z);
    dir = e-s;
    while (1) {
          int t = next();
          if (t==0) break;
          printf("%d ",t);
          cnt++;
          s = s+tt*dir;
          Point R = s-c[t];
          dir = dir-2.0*(R*dir)/(R*R)*R;
          if (cnt==10) break;
    }
    if (cnt==10 && next()) printf("etc.");
    printf("\n");
    return 0;
}