AOJ 1289,PKU 3944 Spherical Mirrors
問題概要
球上の鏡がN個与えられる。(0,0,0)から(u,v,w)の方向にレーザー光線を発射したとき最後に反射する点の座標を求めよ
解法
シミュレート。
球にぶつかった後のレーザー光線の方向ベクトルは二分探索で求めることができる。
実装(C++)
#include <iostream> #include <algorithm> #include <valarray> #include <vector> #include <cmath> using namespace std; struct Circle{ double x,y,z,r; }; struct Point{ double x,y,z; Point(){} Point(double x,double y,double z):x(x),y(y),z(z){} Point(const Circle &c){ x=c.x;y=c.y;z=c.z; } }; Point operator-(const Point &a,const Point &b){ return Point(a.x-b.x,a.y-b.y,a.z-b.z); } Point operator+(const Point &a,const Point &b){ return Point(a.x+b.x,a.y+b.y,a.z+b.z); } Point operator*(const Point &a,const double &b){ return Point(a.x*b,a.y*b,a.z*b); } double dot(const Point &a,const Point &b){ return a.x*b.x+a.y*b.y+a.z*b.z; } const double EPS=1e-7; Point direction; Point now; Circle c[100]; int N; //交点をすべて求める vector<pair<Point,int> > cross_all(){ vector<pair<Point,int> > res; for(int i=0;i<N;i++){ Point t=now-c[i]; double A=dot(direction,direction); double B=2*dot(direction,t); double C=dot(t,t)-c[i].r*c[i].r; double D=B*B-4*A*C; if(D>=0){ double x=(-B+sqrt(D))/2/A; if(x>EPS)res.push_back(make_pair(now+direction*x,i)); x=(-B-sqrt(D))/2/A; if(x>EPS)res.push_back(make_pair(now+direction*x,i)); } } return res; } int main() { for(;cin>>N,N;){ cin>>direction.x>>direction.y>>direction.z; now=Point(0,0,0); for(int i=0;i<N;i++){ cin>>c[i].x>>c[i].y>>c[i].z>>c[i].r; } while(1){ vector<pair<Point,int> > cp=cross_all(); Point next; int next_index=-1; double distance=99999999; for(int i=0;i<(int)cp.size();i++){ double d=dot(cp[i].first-now,cp[i].first-now); if(distance>d){ distance=d; next_index=cp[i].second; next=cp[i].first; } } if(next_index==-1)break; //反射方向を考える //二分探索によって求める double low=0,high=10000; Point u=next-c[next_index]; double C=dot(now-next,now-next); for(int i=0;i<64;i++){ double mid=(low+high)*.5; double B=dot(u*mid+next-now,u*mid+next-now); double A=dot(u*mid,u*mid); if(A+B>C){ high=mid; }else{ low=mid; } } Point Q=u*low+next; direction=now+(Q-now)*2-next; now=next; } cout<<now.x<<" "<<now.y<<" "<<now.z<<endl; } return 0; }