NOI2005 月下柠檬树

[NOI2005] 月下柠檬树

Description

李哲非常非常喜欢柠檬树,特别是在静静的夜晚,当天空中有一弯明月温柔地照亮地面上的景物时,他必会悠闲地坐在他亲手植下的那棵柠檬树旁,独自思索着人生的哲理。

李哲是一个喜爱思考的孩子,当他看到在月光的照射下柠檬树投在地面上的 影子是如此的清晰,马上想到了一个问题:树影的面积是多大呢?

李哲知道,直接测量面积是很难的,他想用几何的方法算,因为他对这棵柠 檬树的形状了解得非常清楚,而且想好了简化的方法。

李哲将整棵柠檬树分成了 $n$ 层,由下向上依次将层编号为 $1,2,…,n$。从第 $1$ 到 $n-1$ 层,每层都是一个圆台型,第 $n$ 层(最上面一层)是圆锥型。对于圆台型, 其上下底面都是水平的圆。对于相邻的两个圆台,上层的下底面和下层的上底面 重合。第 $n$ 层(最上面一层)圆锥的底面就是第 $n-1$ 层圆台的上底面。所有的底面的圆心(包括树顶)处在同一条与地面垂直的直线上。李哲知道每一层的高度为 $h_1,h_2,…,h_n$,第 $1$ 层圆台的下底面距地面的高度为 $h_0$,以及每层的下底面的圆的半径 $r_1,r_2,…,r_n$。李哲用熟知的方法测出了月亮的光线与地面的夹角为 $alpha$。

uInxQf.png

为了便于计算,假设月亮的光线是平行光,且地面是水平的,在计算时忽略树干所产生的影子。李哲当然会算了,但是他希望你也来练练手。

Solution

自适应辛普森

乍一看十分麻烦。。。

由于是平行光,投影一下就变成了下图。

首先,第一个 $h_0$ 可以直接忽略(题目说了忽略树干影子)。其次,圆的半径不会改变,只有高度与圆台投影形成的梯形尺寸需要注意(最后一个圆锥也可看成圆台,三角形看成梯形)。

只要我们得出了各个圆以及梯形的坐标与大小,不必考虑重合与否,直接套上辛普森,问题便迎刃而解。

我们发现,高度为 $H$ 的物体的长度会变为$H \cdot \displaystyle \frac 1 {tan \ \alpha}$,这可以帮助我们确定横坐标。

比较麻烦的是梯形。我们需要知道它的左右端点以及上下底长。它的腰是由相邻两圆的公切线构成的(知道即可,我们并不需要储存这两条线),两底则可以由左右端点距圆心距离和半径通过勾股定理得出。因此我们先要得出左右端点的位置。

观察下图,发现这两个端点离圆心很近,我们将 $AE,CH$ 分别相连,分别过$E,H$ 向 $x$ 轴作垂线,过$A$ 点向$CH$ 作垂线,应用简单的相似三角形知识即可求出左右端点的坐标。这里图就不给出了。

a

至此,我们已经得出了所有圆和梯形的位置与大小。套上辛普森板子即可。

还有,在计算 $F$ 值时梯形可能会麻烦些,但是这只是初中几何题,这里不再赘述。

Code

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
#include <cmath>
#include <cstdio>
#include <algorithm>
#define LL long long
#define sqr(x) ((x) * (x))

const double eps = 1e-8, pi = acos(-1);

inline int Sgn(double x) {
return x < -eps ? -1 : x > eps;
}

struct Circle {
double x, y, r;
Circle() {}
Circle(double _x, double _y, double _r) : x(_x), y(_y), r(_r) {}
};

struct Trape {
double l, r, al, ar;
Trape() {}
Trape(double _l, double _r, double _al, double _ar) : l(_l), r(_r), al(_al), ar(_ar) {}
};

const int N = 1e3;
int n;
double Alpha, L, R, h[N + 5], r[N + 5];
Circle a[N + 5];
Trape t[N + 5];

inline double F(double pos) {
double res = 0;
for (int i = 1; i <= n; ++i)
if (Sgn(a[i].r - fabs(pos - a[i].x)) > 0) {
double tmp = sqrt(sqr(a[i].r) - sqr(fabs(pos - a[i].x)));
res = std::max(res, tmp * 2);
}
for (int i = 1; i <= n; ++i) {
if (Sgn(pos - t[i].l) > 0 && Sgn(t[i].r - pos) > 0) {
double d = pos - t[i].l, len = t[i].r - t[i].l;
double Len = len * t[i].al / (t[i].al - t[i].ar);
double tmp = (Len - d) / Len * t[i].al;
res = std::max(res, tmp * 2);
}
}
return res;
}

inline double calc(double l, double r) {
return (r - l) * (F(l) + F(r) + 4 * F((l + r) / 2)) / 6;
}

double Simpson(double l, double r, double Eps, double ans) {
double Mid = (l + r) / 2, SL = calc(l, Mid), SR = calc(Mid, r);
if (fabs(SL + SR - ans) <= Eps * 15) return SL + SR + (SL + SR - ans) / 15;
return Simpson(l, Mid, Eps / 2, SL) + Simpson(Mid, r, Eps / 2, SR);
}

int main() {
scanf("%d%lf", &n, &Alpha);
for (int i = 0; i <= n; ++i) scanf("%lf", &h[i]);
for (int i = 1; i <= n; ++i) scanf("%lf", &r[i]);
double H = 0, Tan = 1.0 / tan(Alpha), L = 0, R = 0;
for (int i = 1; i <= n; ++i) {
a[i] = Circle(Tan * H, 0, r[i]);
double len = Tan * h[i], ll = (r[i + 1] - r[i]) * r[i] / len, rr = (r[i + 1] - r[i]) * r[i + 1] / len;
t[i] = Trape(Tan * H - ll, Tan * (H + h[i]) - rr, sqrt(sqr(r[i]) - sqr(ll)), sqrt(sqr(r[i + 1]) - sqr(rr)));
H += h[i];
L = std::min(L, a[i].x - a[i].r);
R = std::max(R, a[i].x + a[i].r);
}
R = std::max(R, Tan * H);
printf("%.2lf\n", Simpson(L, R, 5e-4, calc(L, R)));
return 0;
}