sgu 110 射线关于球的反射光线
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 #include <iostream> 2 #include <fstream> 3 #include <sstream> 4 #include <cstdlib> 5 #include <cstdio> 6 #include <cmath> 7 #include <string> 8 #include <cstring> 9 #include <algorithm> 10 #include <queue> 11 #include <stack> 12 #include <vector> 13 #include <set> 14 #include <map> 15 #include <list> 16 #include <iomanip> 17 #include <cctype> 18 #include <cassert> 19 #include <bitset> 20 #include <ctime> 21 22 using namespace std; 23 24 #define pau system("pause") 25 #define ll long long 26 #define pii pair<int, int> 27 #define pb push_back 28 #define mp make_pair 29 #define mp make_pair 30 #define pli pair<ll, int> 31 #define pil pair<int, ll> 32 #define clr(a, x) memset(a, x, sizeof(a)) 33 34 const double pi = acos(-1.0); 35 const int INF = 0x3f3f3f3f; 36 const int MOD = 1e9 + 7; 37 const double EPS = 1e-9; 38 39 /* 40 #include <ext/pb_ds/assoc_container.hpp> 41 #include <ext/pb_ds/tree_policy.hpp> 42 using namespace __gnu_pbds; 43 #define TREE tree<pli, null_type, greater<pli>, rb_tree_tag, tree_order_statistics_node_update> 44 TREE T; 45 */ 46 47 inline double p2(double x) {return x * x;} 48 struct Point { 49 double x, y, z; 50 Point () {} 51 Point (double x, double y, double z) : x(x), y(y), z(z) {} 52 Point operator - (const Point &p) const { 53 return Point(x - p.x, y - p.y, z - p.z); 54 } 55 Point operator + (const Point &p) const { 56 return Point(x + p.x, y + p.y, z + p.z); 57 } 58 Point operator * (const double &k) const { 59 return Point(k * x, k * y, k * z); 60 } 61 double operator | (const Point &p) const { 62 return x * p.x + y * p.y + z * p.z; 63 } 64 void input() { 65 scanf("%lf%lf%lf", &x, &y, &z); 66 } 67 void output() { 68 printf("x = %.6f, y = %.6f, z = %.6f\n", x, y, z); 69 } 70 double len2() { 71 return p2(x) + p2(y) + p2(z); 72 } 73 double len() { 74 return sqrt(len2()); 75 } 76 } s, e; 77 struct Line { 78 Point s, e; 79 Line () {} 80 Line (Point s, Point e) : s(s), e(e) {} 81 void output() { 82 printf("l.s = "); 83 s.output(); 84 printf("l.e = "); 85 e.output(); 86 } 87 }; 88 struct Ball { 89 double x, y, z, r; 90 void input(int i) { 91 scanf("%lf%lf%lf%lf", &x, &y, &z, &r); 92 } 93 } b[55]; 94 double solve(double a, double b, double c) { 95 if (b > 0) return 1e18; 96 double delta = b * b - 4 * a * c; 97 if (delta < -EPS) return 1e18; 98 if (delta < EPS) return -b / (2 * a); 99 return (-b - sqrt(delta)) / (2 * a); 100 } 101 Point reflect(Line l1, Line l2) { 102 double k = (l2.e - l2.s) | (l1.s - l1.e); 103 k = k * 2.0 / (l2.e - l2.s).len2(); 104 return (l2.e - l2.s) * k + (l1.e - l1.s); 105 } 106 int n; 107 int main() { 108 scanf("%d", &n); 109 for (int i = 1; i <= n; ++i) { 110 b[i].input(i); 111 } 112 s.input(), e.input(); 113 int pre_id = 0; 114 vector<int> ans; 115 for (int i = 1; i <= 10; ++i) { 116 double lambda = 1e18; 117 Point v = e - s; 118 int tid = pre_id; 119 for (int j = 1; j <= n; ++j) { 120 if (j == pre_id) continue; 121 Point B = s - Point(b[j].x, b[j].y, b[j].z); 122 double res = solve(v.len2(), 2 * (v.x * B.x + v.y * B.y + v.z * B.z), B.len2() - b[j].r * b[j].r); 123 if (res < lambda) { 124 lambda = res; 125 tid = j; 126 } 127 //cout << "i = " << i << ' ' << res << endl; 128 } 129 if (tid == pre_id) break; 130 ans.pb(tid); 131 pre_id = tid; 132 Point pp = reflect(Line(s, e), Line(Point(b[tid].x, b[tid].y, b[tid].z), s + (e - s) * lambda)); 133 //pp.output(); 134 s = s + (e - s) * lambda; 135 //s.output(); 136 e = s + pp; 137 } 138 for (int i = 0; i < ans.size(); ++i) { 139 printf("%d ", ans[i]); 140 } 141 if (ans.size() == 10) { 142 Point v = e - s; 143 int tid = pre_id; 144 double lambda = 1e18; 145 for (int i = 1; i <= n; ++i) { 146 if (i == pre_id) continue; 147 Point B = s - Point(b[i].x, b[i].y, b[i].z); 148 double res = solve(v.len2(), (v.x * B.x + v.y * B.y + v.z * B.z) * 2.0, B.len2() - b[i].r * b[i].r); 149 if (res < lambda) { 150 lambda = res; 151 tid = i; 152 } 153 } 154 if (tid != pre_id) puts("etc."); 155 } 156 return 0; 157 } 158 /* 159 * 2 160 0 0 0 1 161 1 3 0 1 162 -1 2 0 1 -1 0 163 */
