[SGU110 Dungeon]【计算几何】【模拟】【球面反射】

【题目大意】
给定n个球，与一条射线，然后让你模拟这个射线在这些球之间弹，输出弹到得球的标号。如果弹到<=10个就有多少输出多少。否则，输出前10个+ 'etc.'

【算法分析】
我们用三维向量来解决。
1、首先根据点到球心距离等于半径 + 我们的运动向量列出一个方程，解这个一元二次方程可以得到我们什么时候碰到这个球。
2、然后碰到以后呢。就把射线起点，球上的碰触点，球心这三个点，搞成一个平面，然后再在上面以向量：碰触点->球心作为法线，考虑镜面反射即可。
然后这里我再次用了一个解方程+向量伸缩的方法得到了他反射以后的向量。

关于第二步的反射部分，我们所知道的条件有：入射向量，法线向量。我的做法具体就是先构造出向量tmp=发现向量-入射向量。并使得|tmp|=|入射向量|。然后就可以通过一些几何知识来搞了。

【其它】WA了几遍，因为我解一元二次方程的时候判断无解写得是delta<1e-6，然后如果1e-6<=delta<0的话，我就先将它置为0再进行后面的计算。
我改成delta<0就为无解就AC了。T_T。。。精度什么的，最讨厌了。

【CODE】
#include #include #include #include #include #define dis(A,B) sqrt( (A.x-B.x)*(A.x-B.x)+(A.y-B.y)*(A.y-B.y)+(A.z-B.z)*(A.z-B.z) )
#define Sqr(x) ((x)*(x))
#define max(x,y) ((x)>(y)?(x):(y))
#define min(x,y) ((x)<(y)?(x):(y))
using namespace std;
const int N=1005;
const double eps=1e-6;
int n;
double res1,res2;
struct Circle{double x,y,z,r;}a[N];
struct Point{double x,y,z;}p,fx;

void init(){
     scanf("%d",&n);
     int i,j;
     for (i=1;i<=n;i++)
       cin >> a[i].x >> a[i].y >> a[i].z >> a[i].r;
     cin >> p.x >> p.y >> p.z >> fx.x >> fx.y >> fx.z;
     fx.x-=p.x;
     fx.y-=p.y;
     fx.z-=p.z;
}

Point Get_Point(double rate){
      Point res;
      res.x=fx.x*rate+p.x;
      res.y=fx.y*rate+p.y;
      res.z=fx.z*rate+p.z;
      return res;
}

double Get_Meet_Rate(Circle Cir){
       double a,b,c,delta,x1,x2;
       Point tmp;
       tmp.x=Cir.x-p.x;
       tmp.y=Cir.y-p.y;
       tmp.z=Cir.z-p.z;
       a=Sqr(fx.x)+Sqr(fx.y)+Sqr(fx.z);
       b=-2*(fx.x*tmp.x + fx.y*tmp.y + fx.z*tmp.z);
       c=-Sqr(Cir.r)+Sqr(tmp.x)+Sqr(tmp.y)+Sqr(tmp.z);
       delta=b*b-4*a*c;
       if (delta<0) return 1e22;
       delta=sqrt(delta);
       x1=(-b+delta)/(2*a);
       x2=(-b-delta)/(2*a);
       if (x1<=eps) x1=1e22;
       if (x2<=eps) x2=1e22;
       res1=x1; res2=x2;
       return min(x1,x2);
}

Point opposite(Point tmp){
      Point res;
      res.x=-tmp.x;
      res.y=-tmp.y;
      res.z=-tmp.z;
      return res;
}

void Reflect(int i,double rate){
     Point cross=Get_Point(rate);
     Get_Meet_Rate(a[i]);
     if ( fabs(res1-res2) < eps ){
       p=cross;
       return;
     }
     Point normal_line,tmp;
     normal_line.x=(a[i].x-cross.x)/2;
     normal_line.y=(a[i].y-cross.y)/2;
     normal_line.z=(a[i].z-cross.z)/2;
     double t=Sqr(normal_line.x)+Sqr(normal_line.y)+Sqr(normal_line.z);
     t/=(normal_line.x*fx.x + normal_line.y*fx.y + normal_line.z*fx.z);
     fx.x*=t; fx.y*=t; fx.z*=t;
     normal_line.x*=2;
     normal_line.y*=2;
     normal_line.z*=2;
    
     tmp.x=normal_line.x-fx.x;
     tmp.y=normal_line.y-fx.y;
     tmp.z=normal_line.z-fx.z;
    
     fx=opposite(tmp);
     p=cross;
    
     bool flag=false;
     if (fx.x>10) flag=true;
     if (fx.y>10) flag=true;
     if (fx.z>10) flag=true;
     if (!flag){
       fx.x*=10;
       fx.y*=10;
       fx.z*=10;
     }
}

void work(){
     int i,best,cnt=0;
     double rate,tmp;
     while (1){
           rate=1e20;
           best=0;
           for (i=1;i<=n;i++){
               tmp=Get_Meet_Rate(a[i]);
               if (tmp                 rate=tmp;
                 best=i;
               }
           }
           if (best==0) break;
           Reflect(best,rate);
           cnt++;
           if (cnt>10) break;
           if (cnt!=1) cout << " ";
           cout << best;
     }
     if (cnt>10) cout << " etc." << endl;
            else cout << endl;
}

int main(){
    init();
    work();
    return 0;
}