sgu110：Dungeon

110. Dungeon time limit per test: 0.25 sec.  memory limit per test: 4096 KB The mission of space explorers found on planet M the vast dungeon. One of the dungeon halls is fill with the bright spheres. The explorers find out that the light rays reflect from the surface of the spheres according the ordinary law (the incidence angle is equal to the reflectance angle, the incidence ray, the reflected ray and the perpendicular to the sphere surface lay in the one plane). The ancient legend says that if the light ray will reflect from the spheres in the proper order, than the door to the room with very precious ancient knowledge will open. You are not to guess the right sequence; your task is much simpler. You are given the positions and the radii of the spheres, the place where the laser shot was made and the direction of light propagation. And you must find out the sequence in which the light will be reflected from the spheres. Input The first line of input contains the single integer n (1≤n≤50) - the amount of the spheres. The next n lines contain the coordinates and the radii of the spheres xi, yi, zi, ri (the integer numbers less or equal to 10000 by absolute value). The last line contains 6 real numbers - the coordinates of two points. The first one gives the coordinates of the place of laser shot, and the second gives the direction in which it was made (the second point is the point on the ray). The starting point of the ray lies strictly outside of any sphere. Output Your program must output the sequence of sphere numbers (spheres are numbers from 1 as they was given in input), from which the light ray was reflected. If the ray will reflect more the 10 times, than you must output first 10, then a space and the word 'etc.' (without quotes). Notice: if the light ray goes at a tangent to the sphere you must assume that the ray was reflected by the sphere. Sample Input 1 1 0 0 2 1 0 0 0 0 0 1 Sample Output 1 1 Sample Input 2 2 0 0 2 1 0 0 -2 1 0 0 0 0 0 100 Sample Output 2 1 2 1 2 1 2 1 2 1 2 etc.

一道计算几何题，想清楚了还是不难的。

对于一条光线射到球面上，光源为(x1,y1,z1)，方向向量为t。

设交点为(x0,y0,z0)，必有方程(x0,y0,z0)=(x1,y1,z1)+k*t...①；

又焦点在圆上，圆心为(a1,b1,c1)，必有|(a1,b1,c1)-(x0,y0,z0)|=r...②；

联立①②，可以解出交点(x0,y0,z0)。

连结交点和圆心求出法向量(x0-a1,y0-b1,z0-c1)，做出方向向量t关于法向量的对称向量即可。

方法有很多，这里可以通过做投影的方法来求出对称向量。

具体做法参考http://www.cnblogs.com/graphics/archive/2010/08/03/1791626.html

代码如下：

#include <cstdio>
#include <cmath>
using namespace std;
const int MAXN = 55;
const double INF = 10000005;

int N = 0;
struct sphere
{
  double x, y, z, r;	
}ball[MAXN] = {0, 0, 0, 0};
struct vector
{
  double x, y, z;	
};
typedef vector point;
vector drt = {0, 0, 0}, drtl = {0, 0, 0}, drtn = {0, 0, 0};
point st = {0, 0, 0}; 

point operator + (const point &A, const vector &B)
{
  point C = {A.x+B.x, A.y+B.y, A.z+B.z};
  return C;	
}
point operator - (const point &A, const vector &B)
{
  point C = {A.x-B.x, A.y-B.y, A.z-B.z};
  return C;	
}
vector operator * (const double &k, const vector &B)
{
  vector C = {k*B.x, k*B.y, k*B.z};
  return C;	
}

double dot(const vector &A, const vector &B)
{
  return A.x*B.x+A.y*B.y+A.z*B.z;	
}

int main()
{
  scanf("%d", &N);
  for(int i = 1; i <= N; ++i)
    scanf("%lf%lf%lf%lf", &ball[i].x, &ball[i].y, &ball[i].z, &ball[i].r);
  scanf("%lf%lf%lf%lf%lf%lf", &st.x, &st.y, &st.z, &drt.x, &drt.y, &drt.z);
  drt.x -= st.x;drt.y -= st.y;drt.z -= st.z;
  
  int last = 0, neww = 0;
  for(int i = 1; i <= 11; ++i)
  {
    double k = INF, tmpk = 0;
    double x, y, z, _x, _y, _z, a, b, c;
    double shade;
    neww = 0;
	for(int j = 1; j <= N; ++j)
	{
	  if(last == j) continue;
	  x = st.x-ball[j].x, y = st.y-ball[j].y, z = st.z-ball[j].z;
	  _x = drt.x, _y = drt.y, _z = drt.z;
	  a = _x*_x+_y*_y+_z*_z;
	  b = 2*(x*_x+y*_y+z*_z);
	  c = x*x+y*y+z*z-ball[j].r*ball[j].r;
	  if(b*b < 4*a*c) continue;
	  tmpk = (-b-sqrt(b*b-4*a*c))/a*0.5;
	  if(tmpk > 0)
	  {
	    if(tmpk < k)
	    {
		  k = tmpk;
		  neww = j;	
	    }
	  }
	  else
	  {
	    tmpk = (-b+sqrt(b*b-4*a*c))/a*0.5;
	    if(tmpk > 0 && tmpk < k)
	    {
	      k = tmpk;
		  neww = j;	
	    }
	  }
    }
    if(neww == 0)
      break;	
    if(i <= 10) printf("%d ", neww);
    else 
    {
      printf("etc.");
	  break;	
    }
    last = neww;
    if(ball[last].r == 0) continue;
    st = st+k*drt;
    drt = (-1)*drt;
    drtl.x = st.x-ball[last].x, drtl.y = st.y-ball[last].y, drtl.z = st.z-ball[last].z;
    drtn = dot(drt, drtl)/(drtl.x*drtl.x+drtl.y*drtl.y+drtl.z*drtl.z)*drtl;
    drt = 2*drtn-drt;
  }
  printf("\n");
  return 0;	
}