sgu110计算几何

题意：

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.

题解（向量法）：

1、计算直线与球的交点：dir*k(k>0)表示光线、vec2表示起点到目标球的球心的向量，vec2-dir*k的绝对值等于目标球的半径R，解一元二次方程得到k的值，取较小的非负数解，(x+k*(dx-x), y+k*(dy-y), z+k*(dz-z))即为交点坐标

2、计算反射光线：设交点所在球为i，交点为p，vec1表示球心到起点的向量，s表示光线起点，vec2表示球心到交点的向量，dc表示vec1与vec2的点乘，d表示vec2的长度，vec1在vec2方向上的投影坐标为 v((px - x[i]) / d * dc + x[i], (py - y[i]) / d * dc + y[i], (pz - z[i]) / d * dc + z[i])，v*2-s-p即为反射光线的方向向量

3、循环十次，每次循环之后更新光线向量，如果与任何球都没有交点则直接退出，循环十次后若还与球有交点则输出etc.

附上代码：

/*|vec1-dir*k|=r(x-x1*k)2+(y-y1*k)2+(z-z1*k)2=r2(x1^2+y1^2+z1^2)k^2-(2x0x1+2y0y1+2z0z1)k=r^2-x^2-y^2-z^2*/#include <cstdio>#include <cmath>#define ZERO 1e-10using namespace std;struct ball{double x, y, z, r;}b[55];double x, y, z, xd, yd, zd;int n, last;void init(){scanf("%d", &n);for (int i = 1; i <= n; i++) scanf("%lf%lf%lf%lf", &b[i].x, &b[i].y, &b[i].z, &b[i].r);scanf("%lf%lf%lf%lf%lf%lf", &x, &y, &z, &xd, &yd, &zd);}double js(double x1, double y1, double z1, double dirx, double diry, double dirz, double r){double a = dirx * dirx + diry * diry + dirz * dirz;double b = -(2 * x1 * dirx + 2 * y1 * diry + 2 * z1 * dirz);double c = x1 * x1 + y1 * y1 + z1 * z1 - r * r;if (a >= -ZERO && a <= ZERO){double k1 = - c / b;return k1;}double drt = b * b - 4 * a * c;if (drt < -ZERO) return -1;double gh = sqrt((double)drt);double k1 = (-b + gh) / (2 * a);double k2 = (-b - gh) / (2 * a);if (k1 < -ZERO) return -1;if (k2 < -ZERO) return k1;return k2;}int main(){init();for (int i = 1; i <= 11; i++){double k = 1000000000;int get = 0;for (int j = 1; j <= n; j++){if (j == last) continue;double kk = js(b[j].x - x, b[j].y - y, b[j].z - z, xd - x, yd - y, zd - z, b[j].r);if (kk !=-1 && kk < k){k = kk;get = j;}}if (get == 0) return 0;if (i == 11) break;printf("%d ", get);last = get;double x1, y1, z1, xd1, yd1, zd1, vx, vy, vz, dc, d;struct ball g = b[get];x1 = x + k * (xd - x);y1 = y + k * (yd - y);z1 = z + k * (zd - z);dc = (x1 - g.x) * (x - g.x) + (y1 - g.y) * (y - g.y) + (z1 - g.z) * (z - g.z);d = sqrt((x1 - g.x) * (x1 - g.x) + (y1 - g.y) * (y1 - g.y) + (z1 - g.z) * (z1 - g.z));dc /= d;vx = (x1 - g.x) / d * dc + g.x;vy = (y1 - g.y) / d * dc + g.y;vz = (z1 - g.z) / d * dc + g.z;xd1 = vx * 2 - x;yd1 = vy * 2 - y;zd1 = vz * 2 - z;x = x1;y = y1;z = z1;xd = xd1;yd = yd1;zd = zd1;}printf("etc.\n");return 0;}

写得比较丑不要在意，注意精度问题！！！