bzoj1502: [NOI2005]月下柠檬树

题目描述：

李哲非常非常喜欢柠檬树，特别是在静静的夜晚，当天空中有一弯明月温柔地照亮地面上的景物时，他必会悠闲地坐在他亲手植下的那棵柠檬树旁，独自思索着人生的哲理。李哲是一个喜爱思考的孩子，当他看到在月光的照射下柠檬树投在地面上的影子是如此的清晰，马上想到了一个问题：树影的面积是多大呢？李哲知道，直接测量面积是很难的，他想用几何的方法算，因为他对这棵柠檬树的形状了解得非常清楚，而且想好了简化的方法。李哲将整棵柠檬树分成了 $n$  层，由下向上依次将层编号为 $1,2,…,n$ 。从第 $1$ 到 $n-1$  层，每层都是一个圆台型，第 $n$  层(最上面一层)是圆锥型。对于圆台型，其上下底面都是水平的圆。对于相邻的两个圆台，上层的下底面和下层的上底面重合。第 $n$  层(最上面一层)圆锥的底面就是第 $n-1$ 层圆台的上底面。所有的底面的圆心(包括树顶)处在同一条与地面垂直的直线上。李哲知道每一层的高度为 $h1,h2,…,hn,$ 第 $1$ 层圆台的下底面距地面的高度为 $h0$ ，以及每层的下底面的圆的半径 $r1,r2,…,rn$ 。李哲用熟知的方法测出了月亮的光线与地面的夹角为 $alpha$ 。为了便于计算，假设月亮的光线是平行光，且地面是水平的，在计算时忽略树干所产生的影子。李哲当然会算了，但是他希望你也来练练手。

算法标签：辛普森积分

思路：

模拟光线，找出地面投影出的每一个圆，和每一条圆间的线段。用辛普森积分求面积。

以下代码：

#include<bits/stdc++.h>
#define il inline
#define db double
#define _(d) while(d(isdigit(ch=getchar())))
using namespace std;
const int N=505;
const db pi=acos(-1);
int n;
db alp,h[N];
struct point{
    db x,y;
    point(){};
    point(db a,db b){x=a;y=b;}
}line[N][2];
struct Circle{
    point o;db r;
    point getpoint(db ang){
        return point(o.x+r*cos(ang),o.y+r*sin(ang));
    }
}c[N];
il int read(){
    int x,f=1;char ch;
    _(!)ch=='-'?f=-1:f;x=ch^48;
    _()x=(x<<1)+(x<<3)+(ch^48);
    return f*x;
}
il void getline(Circle a,Circle b,point &p,point &q){
    if(a.r<b.r)swap(a,b);
    if((a.r-b.r)>=fabs(a.o.x-b.o.x))return;
    db ang=acos((a.r-b.r)/fabs(a.o.x-b.o.x));
    if(a.o.x<=b.o.x)p=a.getpoint(ang),q=b.getpoint(ang);
    else q=a.getpoint(pi-ang),p=b.getpoint(pi-ang);
}
db pow(db x){
    return x*x;
}
il db f(db x){
    db res=0;
    for(int i=0;i<n;i++)
        if(line[i][0].x<=x&&line[i][1].x>=x){
            res=max(res,(x-line[i][0].x)*(line[i][1].y-line[i][0].y)/(line[i][1].x-line[i][0].x)+line[i][0].y);
        }
    for(int i=0;i<=n;i++){
        if(pow(c[i].r)-pow(c[i].o.x-x)>=0){
            res=max(res,sqrt(pow(c[i].r)-pow(c[i].o.x-x)));
        }
    }
    return res;
}
il db simpson(db a,db b){
    db c=(a+b)/2.0;
    return (f(a)+4*f(c)+f(b))*(b-a)/6;
}
il db asr(db a,db b,db eps,db A){
    db c=(a+b)/2.0;
    db l=simpson(a,c),r=simpson(c,b);
    if(fabs(l+r-A)<=15*eps)return l+r+(l+r-A)/15.0;
    return asr(a,c,eps/2,l)+asr(c,b,eps/2,r);
}
int main()
{
    n=read();scanf("%lf",&alp);
    db sum=0;
    for(int i=0;i<=n;i++){
        scanf("%lf",&h[i]);
        sum+=h[i];h[i]=sum/tan(alp);
    }
    db r,mn=1e9,mx=0;
    for(int i=0;i<=n;i++){
        if(i!=n)scanf("%lf",&r);
        else r=0;
        mn=min(mn,h[i]-r);
        mx=max(mx,h[i]+r);
        c[i]=(Circle){point(h[i],0),r};
    }
    for(int i=0;i<n;i++){
        getline(c[i],c[i+1],line[i][0],line[i][1]);
    }
    printf("%.2lf
",2*asr(mn,mx,1e-7,simpson(mn,mx)));
    return 0;
}