一道计算几何题,求给定图形面积 题目对联系Simpson的新手很友好,但要注意一些细节 保留两位小数嘛..求近似解就阔以啦

关键词:NOI  计算几何 相关题目:[NOI2005]月下柠檬树

参考:自适应辛普森积分

Description

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

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

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

李哲将整棵柠檬树分成了n 层,由下向上依次将层编号为1,2,…,n。从第1到n-1 层,每层都是一个圆台型,第n层(最上面一层)是圆锥型。对于圆台型,其上下底面都是水平的圆。对于相邻的两个圆台,上层的下底面和下层的上底面重合。第n 层(最上面一层)圆锥的底面就是第n-1层圆台的上底面。所有的底面的圆心(包括树顶)处在同一条与地面垂直的直线上。李哲知道每一层的高度为h1,h2,…,hn,第1层圆台的下底面距地面的高度为h0,以及每层的下底面的圆的半径r1,r2,…,rn。李哲用熟知的方法测出了月亮的光线与地面的夹角为alpha。 月下柠檬树 为了便于计算,假设月亮的光线是平行光,且地面是水平的,在计算时忽略树干所产生的影子。

李哲当然会算了,但是他希望你也来练练手。

Input&Output

Input 第1行包含一个整数n和一个实数alpha,表示柠檬树的层数和月亮的光线与地面夹角(单位为弧度)。 第2行包含n+1个实数h0,h1,h2,…,hn,表示树离地的高度和每层的高度。 第3行包含n个实数r1,r2,…,rn,表示柠檬树每层下底面的圆的半径。

上述输入文件中的数据,同一行相邻的两个数之间用一个空格分隔。 输入的所有实数的小数点后可能包含1至10位有效数字。

1≤n≤500,0.3 < alpha < π/2,0 < hi ≤ 100,0 < ri ≤ 100

Output 输出1个实数,表示树影的面积。四舍五入保留两位小数。

Sample Input

2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output

171.97

Solution

我们求得是类似于下面这样的图形: moon

AC代码

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
#include <cstdio>
#include <cmath>
#include <algorithm>
using namespace std;
#define re register int
#define rep(i,a,b) for(re i=(a);i<=(b);++i)
const int N = 520;
const double eps = 1e-6;
inline double sqr(double x) {return x*x;}

int n, skd=0, tag[N];
double deg, h, ld=1e10, rd=-1e10;
struct Circle {
double x, r;
} t[N];
struct Point { double x, y; };
struct Segment {
Point s, t;
double k, b;
Segment(Point _s = (Point){0,0}, Point _t = (Point){0,0}) {
s = _s, t = _t;
if (s.x>t.x) swap(s,t);
k = (t.y-s.y)/(t.x-s.x);
b = s.y - k*s.x;
};
inline double f(double x) {return k*x+b;}
} l[N];

inline double F(double x) {
double res = 0;
rep(i,1,n)
if(fabs(x - t[i].x) < t[i].r)
res = max(res, sqrt(sqr(t[i].r) - sqr(x-t[i].x)));
rep(i,1,skd)
if(l[i].s.x <= x && l[i].t.x >= x)
res = max(res, l[i].f(x));
return res;
}

inline double Simpson(double l, double r) {
return (r-l)/6.0*(F(l) + 4.0*F((l+r)/2) + F(r));
}

inline double SOL(double l, double r, double ans) {
double mid = (l+r)/2, lans = Simpson(l, mid), rans = Simpson(mid, r);
if(fabs(lans + rans - ans) < eps) return ans;
return SOL(l, mid, lans) + SOL(mid, r, rans);
}

inline void init() {
double sina, cosa;
rep(i,1,n-1) {
double d = t[i+1].x-t[i].x;
if(d - fabs(t[i].r - t[i+1].r) <0.0) continue;
sina = (t[i].r-t[i+1].r) / d;
cosa = sqrt(1.0 - sqr(sina));
l[++skd] = (Segment) {
(Point) {t[i].x + t[i].r*sina, t[i].r*cosa},
(Point) {t[i+1].x + t[i+1].r*sina, t[i+1].r*cosa}
};
} //公切线
rep(i,1,n)
ld = min(ld, t[i].x-t[i].r), rd = max(rd, t[i].x+t[i].r);
}

int main() {
scanf("%d%lf", &n, &deg);
deg = 1.0/tan(deg);
rep(i,1,n+1)
scanf("%lf", &t[i].x), h+=t[i].x, t[i].x = h*deg; //投影转换
rep(i,1,n)
scanf("%lf", &t[i].r);
t[++n].r = 0;
init();
// rep(i,1,skd)
// printf("%lf %lf * %lf %lf\n", l[i].s.x, l[i].s.y, l[i].t.x, l[i].t.y);
printf("%.2lf\n", 2.0*SOL(ld, rd, Simpson(ld, rd)));
}