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 */

View Code

 

转载于:https://www.cnblogs.com/BIGTOM/p/9735616.html