SGU 110 Dungeon（计算几何）

Description
给出一束光线的坐标和入射方向，再给出n个球的球心坐标和半径（保证不重合不相交），这束光纤会在n个球之间反射（满足反射定律），输出反射的序列（即球的编号），注意相切也算反射，如果反射次数超过10次则只需要输出10次，后面输出etc.表示超过10次
Input
第一行一整数n(<=50)表示球的个数，之后n行每行四个数表示每个球的球心坐标和半径，最后一行六个数表示光线的入射起点和入射方向
Output
输出光线在球上的反射序列
Sample Input
2
0 0 2 1
0 0 -2 1
0 0 0 0 0 100
Sample Output
1 2 1 2 1 2 1 2 1 2 etc.
Solution
计算几何，每次从光线与各个球的交点中找到一个离光线位置最近的一个表示光线在这个球上发生反射，之后更新光线的位置和方向即可
Code

#include<cstdio>
#include<iostream>
#include<cmath>
using namespace std;
#define maxn 55
#define eps 1e-10
struct node
{
    double x,y,z,r;
    node operator +(const node &b)const
    {
        node c;
        c.x=x+b.x,c.y=y+b.y,c.z=z+b.z;
        return c;
    }
    node operator -(const node &b)const 
    {
        node c;
        c.x=x-b.x,c.y=y-b.y,c.z=z-b.z;
        return c;
    }
    node operator *(const double n)const    
    {
        node c;
        c.x=x*n,c.y=y*n,c.z=z*n;
        return c;
    }
    double dis(node b)
    {
        double i=x-b.x,j=y-b.y,k=z-b.z;
        return sqrt(i*i+j*j+k*k);
    }
    double len()
    {
        return sqrt(x*x+y*y+z*z);
    }
    double multiply(node b)
    {
        double i=y*b.z-z*b.y,j=-(x*b.z-z*b.x),k=x*b.y-y*b.x;
        return sqrt(i*i+j*j+k*k);
    }
    void unit()
    {
        double l=len();
        x/=l,y/=l,z/=l;
    }
    void out()
    {
        printf("%.2lf %.2lf %.2lf\n",x,y,z);
    }
};
int n,res,pos,ans[maxn];
node a,b,s,p[maxn];
int reflect()
{
    double m=1e100,temp=pos;
    node ta,ts;
    for(int i=0;i<n;i++)    
        if(i!=pos)
        {
            double d=s.multiply((p[i]-a));
            double D=a.dis(p[i]);
            if(d-p[i].r<eps)
            {
                double l=sqrt(D*D-d*d)-sqrt(abs(p[i].r*p[i].r-d*d));
                node c,ss;
                c=a+s*l;
                if(abs(c.dis(p[i])-p[i].r)<eps&&m-l>eps)
                {
                    m=l,temp=i,ta=c;
                    double angle=2.0*sqrt(abs(p[i].r*p[i].r-d*d))/p[i].r;
                    ss=c-p[i],ss.unit(),ss=ss*angle;
                    ts=ss+s,ts.unit();
                }
            }
        }
    if(temp==pos)return 0;
    a=ta,s=ts,pos=temp,ans[res++]=pos+1;
    return 1;
}
int main()
{
    while(~scanf("%d",&n))
    {
        res=0;
        for(int i=0;i<n;i++)
            scanf("%lf%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z,&p[i].r);
        scanf("%lf%lf%lf%lf%lf%lf",&a.x,&a.y,&a.z,&b.x,&b.y,&b.z);
        s=b-a,s.unit(),pos=-1;
        while(res<11)
        {
            if(!reflect())break;
        }
        if(res>10)
        {
            for(int i=0;i<10;i++)printf("%d ",ans[i]);
            puts("etc.");
        }
        else
            for(int i=0;i<res;i++)
                printf("%d%c",ans[i],i==res-1?'\n':' ');
    }
    return 0;
}