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;
}