1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
|
题意:一束光线射向三维空间中的很多个镜面球体,遇到球面发生反射,经过有
限次后光线不会再与球面相遇问最后光线照射到的球面上的坐标
思路:遇到球面发生反射,求出起点关于法线对称点,光线的起点改为照射到的
球面的坐标,射向方向为新起点照向对称点最先遇到的球面上的坐标即为反射点
需先求一次关于时间t的二次方程,t < 0 舍去
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#include<set>
#include<map>
#include<algorithm>
const int maxn = 1e2;
const double pi = acos(-1.0);
const double eps = 1e-5;
using namespace std;
//加法误差
double add(double a, double b) {
if
(fabs(a + b) < eps)
return
0;
return
a + b;
}
//向量运算
struct P {
double x, y, z;
P() {}
P(double x, double y, double z) : x(x), y(y), z(z) {}
P operator + (P p) {
return
P(add(x, p.x), add(y, p.y), add(z, p.z));
}
P operator - (P p) {
return
P(add(x, -p.x), add(y, -p.y), add(z, -p.z));
}
P operator * (double d) {
return
P(x * d, y * d, z * d);
}
//点乘
double dot(P p) {
return
add(add(x * p.x, y * p.y), z * p.z);
}
//面积
double area(P p) {
double dx = add(y * p.z, -z * p.y);
double dy = -add(x * p.z, -z * p.x);
double dz = add(x * p.y, -y * p.x);
double s = add(add(dx * dx, dy * dy), dz * dz);
return
sqrt(s);
}
//向量模长
double lenth() {
return
sqrt(add(add(x * x, y * y), z * z));
}
//单位向量
P unit() {
return
P(x, y, z) * (1 / lenth());
}
//向量p1-p在向量p1-p2上的投影模长
double prj_len(P p1, P p2) {
P p(x, y, z);
return
(p - p1).dot((p2 - p1).unit());
}
//在p1-p2上的投影(垂足)
P prj(P p1, P p2) {
P p(x, y, z);
return
p1 + (p2 - p1).unit() * p.prj_len(p1, p2);
}
//关于p1-p2的对称点
P sym(P p1, P p2) {
P p(x, y, z);
P e = p.prj(p1, p2);
return
e * 2 - p;
}
};
struct dis {
P p;
double d;
int ind;
dis() {}
dis(P p, double d, int ind) : p(p), d(d), ind(ind) {}
};
bool operator < (dis a, dis b) {
return
a.d < b.d;
}
int main() {
int n, num;
double r[maxn];
P st, a[maxn], s;
dis df[maxn];
while
(scanf(
"%d"
, &n) && n) {
int m = -1;
scanf(
"%lf %lf %lf"
, &s.x, &s.y, &s.z);
st.x = st.y = st.z = 0;
for
(int i = 0; i < n; i++) {
scanf(
"%lf %lf %lf %lf"
, &a[i].x, &a[i].y, &a[i].z, &r[i]);
}
while
(1) {
num = 0;
P next = st + s, toy;
for
(int i = 0; i < n; i++) {
if
(i == m)
continue
;
P e = a[i].prj(st, next);
P ss = e - a[i];
double d = ss.lenth();
if
(d > r[i] || !add(r[i], -d))
continue
;
//相切或相离
double dx = s.x, dy = s.y, dz = s.z;
double Dx = add(st.x, -a[i].x), Dy = add(st.y, -a[i].y), Dz = add(st.z, -a[i].z);
// A*t^2 + B*t + C = 0
//printf("dxyz = %lf %lf %lf\n", dx, dy, dz);
//printf("Dxyz = %lf %lf %lf\n", Dx, Dy, Dz);
double A = add(add(dx * dx, dy * dy), dz * dz);
double B = 2 * add(add(dx * Dx, dy * Dy), dz * Dz);
double C = add(add(Dx * Dx, Dy * Dy), add(Dz * Dz, -r[i] * r[i]));
double lf = add(B * B, -4 * A * C);
//printf("Number.%d dt = %lf\n\n", i, sqrt(lf));
if
(lf <= 0)
continue
;
double dd = sqrt(lf);
double t1 = add(-B, dd) / 2 / A, t2 = add(-B, -dd) / 2 / A;
P p1 = P(add(t1 * dx, st.x), add(t1 * dy, st.y), add(t1 * dz, st.z));
double d1 = (p1 - st).lenth();
P p2 = P(add(t2 * dx, st.x), add(t2 * dy, st.y), add(t2 * dz, st.z));
double d2 = (p2 - st).lenth();
if
(s.dot(p1 - st) > 0) df[num++] = dis(p1, d1, i);
if
(s.dot(p2 - st) > 0) df[num++] = dis(p2, d2, i);
}
if
(!num)
break
;
sort(df, df + num);
m = df[0].ind;
P e = st.sym(a[df[0].ind], df[0].p);
st = df[0].p;
s = e - st;
P p2 = e;
}
printf(
"%.4lf %.4lf %.4lf\n"
, st.x, st.y, st.z);
}
return
0;
}
|
POJ3944 Spherical Mirrors
最新推荐文章于 2020-03-06 15:47:33 发布