HDU 2912

直线关于球的多次反射，求最后一次反射点

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>

using namespace std;
const double inf=1e10;
const double eps=1e-8;
struct point {
    double x,y,z;
//    point (double _x,double _y,double _z){ x=_x; y=_y; z=_z; };
};
struct sphe {
    point cent;
    double r;
};
struct vect {
    point st,des;
};
sphe cir[110];
 vect livc;
int n;

point operator -(const point &u,const point &v){
    point ret;
    ret.x=u.x-v.x; ret.y=u.y-v.y; ret.z=u.z-v.z;
    return ret;
}

double dot(point x,point y){
    return x.x*y.x+x.y*y.y+x.z*y.z;
}

point xmulti(point u,point v){
    point ret;
    ret.x=(u.y*v.z-v.y*u.z);
    ret.y=(u.z*v.x-u.x*v.z);
    ret.z=(u.x*v.y-u.y*v.x);
    return ret;
}

double dis(point x,point y){
    return sqrt((x.x-y.x)*(x.x-y.x)+(x.y-y.y)*(x.y-y.y)+(x.z-y.z)*(x.z-y.z));
}

double vlen(point x){
    return sqrt(x.x*x.x+x.y*x.y+x.z*x.z);
}

point  construct(){
    point crop;
    crop.x=crop.y=crop.z=0;
    double stoc=inf; point tmpcrop; point foot,tmpfoot; bool flag; point tmp; int k;
    while(true){
        flag=false; stoc=inf;
        for(int i=0;i<n;i++){
            if(dot(livc.des-livc.st,cir[i].cent-livc.st)>=-eps){//判断圆是否与直线同向，通过点积判方向 
                double D=vlen(xmulti(livc.des-livc.st,cir[i].cent-livc.st))/dis(livc.st,livc.des);  
            //    cout<<D<<' '<<i<<endl;
                if(D-cir[i].r<=eps){       //半径小于D,相交 
                    flag=true;
                //    cout<<"YES"<<endl;
                    double u=dot(cir[i].cent-livc.st,livc.des-livc.st)/(dis(livc.st,livc.des)*dis(livc.st,livc.des));
                    	//计算垂足。可通过向量的比例所得方程，联合垂直点积为0的方程解得 
                    tmpfoot=livc.st;
                    tmpfoot.x+=u*(livc.des.x-livc.st.x);
                    tmpfoot.y+=u*(livc.des.y-livc.st.y);
                    tmpfoot.z+=u*(livc.des.z-livc.st.z);
                //    cout<<tmpfoot.x<<' '<<tmpfoot.y<<' '<<tmpfoot.z<<' '<<endl;
                    u=sqrt((cir[i].r*cir[i].r-D*D))/dis(livc.st,livc.des);  //计算交点。垂足到圆上交点方向与直线反方向相同
												//通过两者距离比计算出向量的转化 
                    tmpcrop=tmpfoot;
                    tmp=livc.st-livc.des;
                    tmpcrop.x+=tmp.x*u;
                    tmpcrop.y+=tmp.y*u;
                    tmpcrop.z+=tmp.z*u;
                    D=dis(tmpcrop,livc.st);
                //    cout<<D<<endl;
                    if(D<stoc){     //若与多个圆相交，选取较近的一个 
                        stoc=D; crop=tmpcrop;
                        k=i;
                    }
                }
            }
        }
        if(!flag) return crop;
        double tu=dot(livc.st-cir[k].cent,crop-cir[k].cent)/(dis(crop,cir[k].cent)*dis(crop,cir[k].cent));
        tmpfoot=cir[k].cent;     //计算反射线。直线st点关于交点与球心的直线 对称点作为反射线的des点 
        tmpfoot.x+=tu*(crop.x-cir[k].cent.x);
        tmpfoot.y+=tu*(crop.y-cir[k].cent.y);
        tmpfoot.z+=tu*(crop.z-cir[k].cent.z);    //知直线st点到反射线des点的方向与st点到关于对称线垂足方向相同且为两倍 
        livc.des.x=((tmpfoot.x-livc.st.x)*2+livc.st.x);  //通过这样可以求对称点 
        livc.des.y=((tmpfoot.y-livc.st.y)*2+livc.st.y);
        livc.des.z=((tmpfoot.z-livc.st.z)*2+livc.st.z);
        livc.st=crop;
    //    cout<<livc.des.x<<' '<<livc.des.x<<' '<<livc.des.x<<endl;
    }
}

int main(){
    point tmp; double r;
    while(scanf("%d",&n),n){
        livc.st.x=livc.st.y=livc.st.z=0;
        scanf("%lf%lf%lf",&tmp.x,&tmp.y,&tmp.z);
        livc.des=tmp;
        for(int i=0;i<n;i++){
            scanf("%lf%lf%lf%lf",&cir[i].cent.x,&cir[i].cent.y,&cir[i].cent.z,&cir[i].r);
        }
        tmp=construct();
        printf("%.4lf %.4lf %.4lf
",tmp.x,tmp.y,tmp.z);
    }
}