sgu 110 Dungeon 三维计算几何

题意：宇宙探险家们在一次任务中发现了M行星上的一个巨大的地牢，地牢的一个大厅中充满了闪闪发光的小球。探险家们发现光线

在小球表面反射时遵循光的镜面反射规律（入射角=反射角）。根据古老的传说，如果将光线以恰当的顺序在小球之间反射，那么通

道一间隐藏着古老的珍贵知识的房间的门就会开启。放心，你的任务不是去猜这个恰当的顺序是什么。你的任务简单得多，你被告知

小球的位置和半径，以及入射光线的光源位置和方向，然后你要找出光线在小球上反射的序列。

思路：设初始点为s，方向向量为dir，那么如果与圆有交点，必然为交点p可以表示为s+dir*k，其中k>0 。那么我们可以暴力枚举所

有的圆，如何判断与该圆相交呢？ 假如相交，设起点s到圆心的向量为vec，起点s到交点p的向量为dir*k,那么两个向量相减的绝对值

一定为该圆的半径，于是列出一元二次方程，解出k(k>0)。那么在所有的解中取最小的k，不妨设为kmin，那么交点p=dir*kmin+s。

那么如何求出反射的方向向量呢？设交点p所在的圆为i,从该圆圆心到起点s的向量设为vec1，从该圆圆心到交点p的向量为vec2，

我们对vec1、vec2求点积再除以vec2的长度就可以得到从圆心到起点s投影到vec2上点（该点即为s垂直于向量vec2的垂足）的向

量，不妨设为vec3，那么vec3+圆心坐标就是垂足的坐标，设为p1点。那么可以得到s相对与p1对称的点，不妨设点为t，那么反射向

量就是从交点p指向t的向量。那么循环10次。。每次循环之后如果找到点就用交点更新起点，用反向向量更新方向向量，找不到就直

接break。有一点需要注意，对于上一次循环找到的圆在这一次扫的时候要屏蔽，不然永远落在该圆上。详见代码：

// file name: sgu110.cpp //
// author: kereo //
// create time:  2014年10月16日 星期四 19时13分02秒 //
//***********************************//
#include<iostream>
#include<cstdio>
#include<cstring>
#include<queue>
#include<set>
#include<map>
#include<vector>
#include<stack>
#include<cmath>
#include<string>
#include<algorithm>
using namespace std;
typedef long long ll;
const int MAXN=100+100;
const double eps=1e-8;
const int inf=0x3fffffff;
const int mod=1000000000+7;
#define L(x) (x<<1)
#define R(x) (x<<1|1)
int n;
struct point{
	double x,y,z,r;
	point(double x=0,double y=0,double z=0,double r=0) : x(x),y(y),z(z),r(r){ }
}p[MAXN];
typedef point vec;
vec operator + (vec a,vec b){
	return vec(a.x+b.x,a.y+b.y,a.z+b.z);
}
vec operator - (vec a,vec b){
	return vec(a.x-b.x,a.y-b.y,a.z-b.z);
}
vec operator * (vec a,double k){
	return vec(a.x*k,a.y*k,a.z*k);
}
vec operator / (vec a,double k){
	return vec(a.x/k,a.y/k,a.z/k);
}
int dcmp(double x){
	if(fabs(x)<eps) return 0;
	return x < 0 ? -1 : 1;
}
double cross(vec a,vec b){
	return a.x*b.x+a.y*b.y+a.z*b.z;
}
bool judge(point s,vec dir,point cir,double &k){
	double a=dir.x*dir.x+dir.y*dir.y+dir.z*dir.z;
	double b=2*(dir.x*(s.x-cir.x)+dir.y*(s.y-cir.y)+dir.z*(s.z-cir.z));
	double c=(s.x-cir.x)*(s.x-cir.x)+(s.y-cir.y)*(s.y-cir.y)+(s.z-cir.z)*(s.z-cir.z)-cir.r*cir.r;
	double deta=b*b-4*a*c;
	if(dcmp(deta)<0) return false;
	double x1=(-b-sqrt(deta))/(2*a);
	double x2=(-b+sqrt(deta))/(2*a);
	if(dcmp(x2)<0)
		return false;
	k=x2;
	if(dcmp(x1)>0)
	k=x1;
	return true;
}
vec reflect(point s,point p,point cir){
	point tmp=cir+((p-cir)*cross(s-cir,p-cir)/cross(p-cir,p-cir));
	return point(2*tmp.x-s.x,2*tmp.y-s.y,2*tmp.z-s.z)-p;
}
int main()
{
	while(~scanf("%d",&n)){
		for(int i=0;i<n;i++)
			scanf("%lf%lf%lf%lf",&p[i].x,&p[i].y,&p[i].z,&p[i].r);
		point s; vec dir;
		scanf("%lf%lf%lf",&s.x,&s.y,&s.z);
		scanf("%lf%lf%lf",&dir.x,&dir.y,&dir.z);
		dir.x-=s.x; dir.y-=s.y; dir.z-=s.z;
		int pre=-1; //记录之前撞的
		for(int i=0;i<11;i++){
			int now=-1;
			double Min=1e50;
		   	for(int j=0;j<n;j++) if(j!=pre){
				double t;
				if(judge(s,dir,p[j],t) && Min>t)
				   now=j,Min=t;	
			}
			if(now == -1){
				printf("\n");
				break;
			}
			if(i == 10){
				printf(" etc.\n");
				break;
			}
			else{
				if(i)
					printf(" ");
					point tmp=s+dir*Min;
					dir=reflect(s,tmp,p[now]);
					s=tmp;
				printf("%d",now+1);
			}
			pre=now;	
		}
	}
	return 0;
}