1502: [NOI2005]月下柠檬树

1502: [NOI2005]月下柠檬树

Time Limit: 5 Sec  Memory Limit: 64 MB
Submit: 1077  Solved: 600
[Submit][Status][Discuss]

Description

Input

文件的第1行包含一个整数n和一个实数alpha，表示柠檬树的层数和月亮的光线与地面夹角(单位为弧度)。第2行包含n+1个实数h0,h1,h2,…,hn，表示树离地的高度和每层的高度。第3行包含n个实数r1,r2,…,rn，表示柠檬树每层下底面的圆的半径。上述输入文件中的数据，同一行相邻的两个数之间用一个空格分隔。输入的所有实数的小数点后可能包含1至10位有效数字。

Output

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

Sample Input

2 0.7853981633
10.0 10.00 10.00
4.00 5.00

Sample Output

171.97

HINT

1≤n≤500，0.3

Source

[Submit][Status][Discuss]



首先要解决的大概是。。投影的形状如何？

很好想象的是，一个圆形的投影还是一个圆形，而且大小形状完全一样

而一个圆台可以当作无限个半径依次缩小的圆堆叠而成，

所以圆台的投影，是无限个离得很近，半径逐渐减小的圆，只取上下底面的投影圆作关键圆标出，大致如下

一个圆台的投影大概是两个圆以及公切线夹在它们中间的那一段

所以整个图形的投影是一段连续的图像，可以抽象成一些圆和一些公切线组成的图形

这就可以当作一个连续的函数了求函数图像下方的面积，可以使用微积分，那就用辛普森法则把它积出来

对于在区间[a,b]的函数图像，它的积分可以近似地下列方法替换

用一条过图像起止点及中间点的一条抛物线作曲边，则这个曲边梯形的面积即视作积分，有公式

不过如果区间过长，这样的式子显然有极大误差

所以每次求完当前区间的积分后，看看左右子区间积分是否和当前区间相似

如果相似就直接返回值了，否则继续细分下去

f(x)的求法，可以暴力枚举每个圆和每条公切线，看在哪个图像上取最值

公切线的求法可以简单地用相似三角形配合三角函数做

#include<iostream>#include<cstdio>#include<cstring>#include<cmath>using namespace std; const int maxn = 505;typedef double DB;const DB EPS = 1E-7;const DB INF = 1E16; struct Point{    DB x,y; Point(){}    Point(DB x,DB y): x(x),y(y){}    Point operator * (const DB &t) {return Point(x * t,y * t);}    Point operator + (const Point &B) {return Point(x + B.x,y + B.y);}};typedef Point Vector; struct Circle{    Point o; DB r;}c[maxn]; struct Segment{    Point l,r; DB k,b; Segment(){}    Segment(Point l,Point r): l(l),r(r)    {        k = (r.y - l.y) / (r.x - l.x);        b = l.y - k * l.x;    }}s[maxn]; int n,cnt;DB alpha,sum; int cmp(DB x,DB y){    if (fabs(x - y) <= EPS) return 0;    return x < y ? -1 : 1;} DB Dot(Vector v1,Vector v2) {return v1.x * v2.x + v1.y * v2.y;}DB Length(Vector v) {return sqrt(Dot(v,v));}DB Cal(DB c,DB a) {return sqrt(c * c - a * a);} DB F(DB x){    DB ret = 0;    for (int i = 1; i < n; i++)    {        DB d = fabs(c[i].o.x - x);        if (cmp(d,c[i].r) <= 0) ret = max(ret,Cal(c[i].r,d));    }    for (int i = 1; i <= cnt; i++)        if (cmp(s[i].l.x,x) <= 0 && cmp(x,s[i].r.x) <= 0)            ret = max(ret,s[i].k * x + s[i].b);    return ret;} DB Calc(DB L,DB R) {    DB A = (R - L) / 6.00;    DB B = (F(L) + 4.00 * F((L + R) / 2.00) + F(R));    return A * B;} DB Simpson(DB L,DB R,DB now){    DB mid = (L + R) / 2.00;    DB A = Calc(L,mid),B = Calc(mid,R);    if (!cmp(A + B,now)) return now;    return Simpson(L,mid,A) + Simpson(mid,R,B);} DB Work(){    DB L = c[n].o.x - c[n].r,R = c[n].o.x + c[n].r;    for (int i = 1; i < n; i++)    {        L = min(L,c[i].o.x - c[i].r); R = max(R,c[i].o.x + c[i].r);        if (cmp(c[i+1].o.x - c[i+1].r,c[i].o.x - c[i].r) <= 0) continue;        if (cmp(c[i+1].o.x + c[i+1].r,c[i].o.x + c[i].r) <= 0) continue;        if (!cmp(c[i].r,c[i+1].r))        {            Point l = c[i].o + Point(0,c[i].r);            Point r = c[i+1].o + Point(0,c[i+1].r);            s[++cnt] = Segment(l,r); continue;        }        DB d = c[i+1].o.x - c[i].o.x,r1 = c[i].r,r2 = c[i+1].r;        if (r1 > r2) swap(r1,r2); DB delta = r1 * d / (r2 - r1);        DB theta = asin(r2 / (delta + d)),len; Vector v; Point p,l,r;        if (c[i].r < c[i+1].r) p = Point(c[i].o.x - delta,0),v = Point(1,tan(theta));        else p = Point(c[i+1].o.x + delta,0),v = Point(-1,tan(theta));        len = Length(v);        l = p + v * (Cal(delta,r1) / len);        r = p + v * (Cal(delta + d,r2) / len);         if (l.x > r.x) swap(l,r); s[++cnt] = Segment(l,r);    }    return Simpson(L,R,Calc(L,R));} int main(){    #ifdef DMC        freopen("DMC.txt","r",stdin);    #endif         cin >> n >> alpha; ++n;    for (int i = 1; i <= n; i++)    {        DB h; scanf("%lf",&h);        sum += h; c[i].o = Point(sum / tan(alpha),0);    }    for (int i = 1; i < n; i++) scanf("%lf",&c[i].r); c[n].r = 0;    printf("%.2lf\n",Work() * 2.00);    return 0;}