sgu110



SGU110 Dungeon
题目大意
空间探测器在星球M上发现了巨大的地牢，地牢被明亮的球充满，探测器发现光线能按自然规律被球表面反射（入射角等于反射角，入射光线、反射光线、法线在同一平面）。古老的传说说如果光按一定顺序被球表面反射，房间的门就会打开。
你不需要去猜这个顺序；你的任务更简单一些。你会知道球的位置和半径、激光发射的位置及光传播的方向。你要找出光被球反射的顺序。
输入

n(1<=n<=50)球的数量 下面n行读入球坐标，半径xi, yi, zi, ri(integer范围内)。
最后一行包含6个数。
前三个是激光发射的坐标(发射点严格的在任何球体外)。
后三个告诉你激光发射的方向(这个坐标在光线上)。
输出

光被球反射的顺序(球在输入中从1开始依次编号)。
如果球被反射10次以上，输出前十次，然后用一个空格和'etc.'隔开（不含引号）。
注意:如果光的轨迹是某个球的切线,认为光被球反射。
样例输入1

1
0 0 2 1
0 0 0 0 0 1

样例输出1

1

样例输入2

2
0 0 2 1
0 0 -2 1
0 0 0 0 0 100

样例输出2

1 2 1 2 1 2 1 2 1 2 etc.

第一道吓到我的题目。
首先想到的是解析式求解交点，不过扩展到三维的时候成功挂掉。。。
更改思路使用向量，就可以成功AC。

对于点(xs,ys,zs)和光线向量(xl,yl,zl)，在这条光线上的点(x,y,z)满足：
x=xs+k*xl;
y=ys+k*yl;
z=zs+k*zl;
其中，k∈R

对于球(x0,y0,z0)和半径r，在这个球上的点(x,y,z)满足：
(x0-x)^2+(y0-y)^2+(z0-z)^2=r^2;

联立上述四个式子，可以得到一个关于k的一元二次方程。
若该方程无解，则证明两者不相交。

在得到交点之后，问题就是反射了(请先行了解向量加减法的几何意义和光的反射定律)
法线向量v1为从球心指向交点的向量。
入射的反向量v2为从交点指向出发点的向量。
我们可以利用公式求出v1在v2上的投影向量v3(投影公式参考网友博客：向量投影)。
然后将(v3-v1)*2+v1得到反射向量(无法理解的话可以在二维平面上模拟画图帮助理解)。
反射出发点变为交点。

注意事项：1.一元二次方程的解可能有两个，我们要取的是离出发点进的那个交点。
          2.光线是一条射线而不是直线，所以一元二次方程的解应该>0。
          3.在变换向量时入射向量不等于光线向量(长度不同)，不能混用。
          4.开始时记得计算光线向量(输入所给并非光线向量，而是这上面的一个点)。

下面附上我的代码：

#include <stdio.h>
#include <stdlib.h>
#include <stdbool.h>
#include <math.h>
struct b
{
  double x,y,z;
  double r;
  bool pass;
}ball[51],now,line,next,zero;
int n;
double dis(struct b a1,struct b a2)
{
  return (sqrt((a1.x-a2.x)*(a1.x-a2.x)+(a1.y-a2.y)*(a1.y-a2.y)+(a1.z-a2.z)*(a1.z-a2.z)));
}
struct b equation(double a,double b,double c) //解方程ak^2+bk+c=0和对应的交点
{
  double k1,k2;
  struct b x1,x2;
  x1.pass=false;
  k1=(-b+sqrt(b*b-4*a*c))/(2*a);
  x1.x=k1*line.x+now.x;
  x1.y=k1*line.y+now.y;
  x1.z=k1*line.z+now.z;
  k2=(-b-sqrt(b*b-4*a*c))/(2*a);
  x2.x=k2*line.x+now.x;
  x2.y=k2*line.y+now.y;
  x2.z=k2*line.z+now.z;
  if ((k1==k2 || dis(now,x1)<dis(now,x2)) && k1>0)
    {
    x1.pass=true;
    return x1;
    }
  else if (k2>0)
    {
    x2.pass=true;
    return x2;
    }
  return x1;
}
struct b find(int k) //判断交点
{
  double a,b,c;
  struct b re;
  if (k==2)
    k=2;
  re.pass=false;
  a=line.x*line.x+line.y*line.y+line.z*line.z;
  b=2*line.x*(now.x-ball[k].x)+2*line.y*(now.y-ball[k].y)+2*line.z*(now.z-ball[k].z);
  c=(now.x-ball[k].x)*(now.x-ball[k].x)+(now.y-ball[k].y)*(now.y-ball[k].y)+(now.z-ball[k].z)*(now.z-ball[k].z)-ball[k].r*ball[k].r;
  if (a!=0 && b*b>=4*a*c)
    re=equation(a,b,c);
  else if (a==0 && b!=0 && (-c)/b>0)
    {
	re.x=(-c)/b*line.x+now.x;
    re.y=(-c)/b*line.y+now.y;
    re.z=(-c)/b*line.z+now.z;
    re.pass=true;
	}
  return re;
}
double point(struct b a1,struct b a2)
{
  return a1.x*a2.x+a1.y*a2.y+a1.z*a2.z;
}
struct b makev(struct b a1,struct b a2)
{
  struct b re;
  re.x=a2.x-a1.x;
  re.y=a2.y-a1.y;
  re.z=a2.z-a1.z;
  return re;
}
struct b change(struct b a1,double i)
{
  struct b re;
  re.x=a1.x*i;
  re.y=a1.y*i;
  re.z=a1.z*i;
  return re;
}
struct b add(struct b a1,struct b a2)
{
  struct b re;
  re.x=a1.x+a2.x;
  re.y=a1.y+a2.y;
  re.z=a1.z+a2.z;
  return re;
}
void turn(int k)
{
  struct b v1,v2,v3,v4,v5;
  double i;
  v1=makev(next,now);  //从交点向起始点指一个向量
  v2=makev(ball[k],next);  //从球心向交点指一个向量
  i=point(v1,v2)/(dis(zero,v2)*dis(zero,v2));
  v3=change(v2,i);  //以上两步求出投影向量v3
  v4=makev(v1,v3);  //从v1向v3指一个向量(实际就是减法)
  v5=change(v4,2);  //将上面的向量延伸
  line=add(v1,v5);  //计算反射向量
  return ;
}
void work()
{
  int t,i,ballk,lastball=0;
  struct b k;
  t=1;
  while (t<=11)
    {
	next.pass=false;
	for (i=1;i<=n;i++)
	  if (i!=lastball) //因为起点在球上，特判
	    {
	    k=find(i); //对第i个球进行判断
	    if (k.pass && (!next.pass || dis(k,now)<dis(next,now))) //判断交点更新
	      {
		  next=k;
		  next.pass=true;
		  ballk=i;
		  }
        }
    if (!next.pass)
      break;
    turn(ballk); //光线反射
    now=next;  //更改入射点
    if (t<=10)
      if (t==1)
        printf("%d",ballk);
      else
        printf(" %d",ballk);
    else
      printf(" etc.");
    t++;
    lastball=ballk; //记录上一个到达的球
	}
  printf("\n");
  return ;
}
void init()
{
  int i;
  scanf("%d",&n);
  for (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",&now.x,&now.y,&now.z,&line.x,&line.y,&line.z);
  line.x-=now.x;
  line.y-=now.y;
  line.z-=now.z; //计算向量
  work();
  return ;
}
int main()
{
  init();
  return 0;
}