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\n",tmp.x,tmp.y,tmp.z);}
}