AOJ 1289 Spherical Mirrors

問題

3次元空間上にn個の球があり、それぞれの半径および中心の座標が与えられる。
今、空間の原点(0, 0, 0)から、与えられた向き(u, v, w)の方向へレーザーを発射する。

レーザーは、球面にぶつかると、入射角と反射角が等しくなるように反射する。
レーザーが最後に跳ね返る点の座標を求めよ。

制約条件

n≦100
レーザーが反射する回数は5回以下

方針

図を描いて、がんばって実装する。
3次元ベクトルの内積、三次元の直線と点の距離、反射などが必要。

ソースコード

P operator+(P a, const  P &b) { rep(i, a.size()) a[i] += b[i]; return a; }
P operator-(P a, const  P &b) { rep(i, a.size()) a[i] -= b[i]; return a; }
P operator*(double d, P p) { rep(i, p.size()) p[i] *= d; return p; }

double dot(const P &a, const P &b){
	double res = 0;
	rep(i, a.size()) res += a[i] * b[i];
	return res;
}
double norm(const P &a) { return dot(a, a); }
double abs(const P &a) { return sqrt(dot(a, a)); }
P projection(const L &l, const P &p){
	double t = dot(p - l.first, l.first - l.second) / norm(l.first - l.second);
	return l.first + t * (l.first - l.second);
}
P reflection(const L &l, const P &p) {
  return p + 2.0 * (projection(l, p) - p);
}
double distanceLP(const L &l, const P &p) {
  return abs(p - projection(l, p));
}

int n;
double u, v, w;
double x[100], y[100], z[100], r[100];

int main(){
  while(cin >> n, n){
		P last(3), dir(3);
		vector<P> ps;
		cin >> dir[0] >> dir[1] >> dir[2];
		rep(i, n){
			cin >> x[i] >> y[i] >> z[i] >> r[i];
			P p;
			p.pb(x[i]); p.pb(y[i]); p.pb(z[i]);
			ps.pb(p);
		}
  	while(1){
  		double mn = INF;
  		int mni = -1;
  		rep(i, n){
  			const L l = mp(last, dir);
  			if(distanceLP(l, ps[i]) >= r[i]) continue;
  			double lo = 0, mid;
  			double hi = dot(ps[i] - l.first, l.second - l.first) / norm(l.second - l.first);
  			rep(it, 100){
  				mid = (lo + hi) / 2;
  				if(abs(l.first + mid * (l.second - l.first) - ps[i]) < r[i]) hi = mid;
  				else lo = mid;
  			}
  			if(lo < EPS) continue;
  			if(mn > mid) mn = mid, mni = i;
			}
			if(mni < 0) break;
			P next = last + mn * (dir - last);
			P ndir = reflection(mp(ps[mni], next), last);
			last = next;
			dir = ndir;
  	}
  	printf("%.4f %.4f %.4f\n", last[0], last[1], last[2]);
	}
  return 0;
}