[NOI2005]月下柠檬树[计算几何（simpson）]

1502: [NOI2005]月下柠檬树

Time Limit: 5 Sec  Memory Limit: 64 MB
Submit: 1169  Solved: 626
[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

求一棵树(圆锥加圆台组成)在平面上的投影的面积。

给定投影角度(0.3 < alpha <= pi/2)。

先来想想圆的投影是什么样子

还是他自己。

再想圆锥投影是什么样子

一个点加一个圆，并且有这个点与该圆的两条切线(该点在圆内部时没有切线)

再想圆台

两个圆，加上两个圆的外公切线组成的一坨图形。

不妨随意画一个。

好难画- -!

大概就转化成这个样子了。

观察这个图形…

轴对称啊- -!

首先AC长是圆心距，可求。

AI长是半径差，可求。

所以CI可求。

连接FC，观察△FAC

2*S△FAC=FG*AC=CI*AF

AF为半径，已知。

所以FG可求。

于是AG可求。

A点坐标已知，所以F点坐标已知。

E点，直接相似即可。

#include<cmath>

#include<cstdio>

#include<algorithm>

#define pf(x) ((x)*(x))

using namespace std;

const int N=+;

const double eps=1e-;

typedef pair<double,double> point;

typedef pair<double,double> circle;

struct line{

    point s,t;

    double k,b;

    line(){}

    line(point _s,point _t){

        s=_s;t=_t;

        k=(s.second-t.second)/(s.first-t.first);

        b=s.second-s.first*k;

    }

    const double f(const double x){

        return k*x+b;

    }

};

int n,n1;double alpha,H[N];

point p;line L[N];circle C[N];

double lb=2e9,rb;

double sina,cosa,tana;

inline void add(const circle &a,const circle &b){

    n1++;

    sina=(a.second-b.second)/(b.first-a.first);

    cosa=sqrt(-pf(sina));

    tana=sina/cosa;

    L[n1].s=make_pair(a.first+a.second*sina,a.second*cosa);

    L[n1].t=make_pair(b.first+b.second*sina,b.second*cosa);

    L[n1].k=-tana;

    L[n1].b=L[n1].s.second-L[n1].s.first*L[n1].k;

}

inline const double F(const double x){

    double re=;

    for(int i=;i<=n1;i++) if(x>=L[i].s.first&&x<=L[i].t.first) re=max(re,L[i].f(x));

    for(int i=;i<=n;i++) if(x>=C[i].first-C[i].second&&x<=C[i].first+C[i].second) re=max(re,sqrt(pf(C[i].second)-pf(x-C[i].first)));

    return re;

}

inline const double simpson(const double l,const double r){

    double mid=(l+r)/;

    return (F(l)+F(r)+*F(mid))*(r-l)/;

}

inline double asr(double l,double r,double eps,double last){

    double mid=(l+r)/;

    double L=simpson(l,mid),R=simpson(mid,r);

    if(fabs(L+R-last)<=*eps) return L+R+(L+R-last)/;

    return asr(l,mid,eps/,L)+asr(mid,r,eps/,R);

}

inline int cmp(const double x){

    if(fabs(x)<eps) return ;

    return x>?:-;

}

int main(){

    scanf("%d%lf",&n,&alpha);

    for(int i=;i<=n+;i++) scanf("%lf",&H[i]),H[i]+=H[i-];

    for(int i=;i<=n;i++) scanf("%lf",&C[i].second);

    double ta=tan(alpha);

    p=make_pair(H[n+]/ta,);rb=max(rb,p.first);

    double x,r,l,h;

    C[n].first=H[n]/ta;

    x=C[n].first;r=C[n].second;

    lb=min(lb,x-r);

    rb=max(rb,x+r);

    if(x+r<p.first){

        l=pf(r)/(p.first-x);// 射影定理

        h=sqrt(pf(r)-pf(l));

        L[++n1]=line(make_pair(x+l,h),p);

    }

    for(int i=n-;i;i--){

        C[i].first=H[i]/ta;

        x=C[i].first;r=C[i].second;

        lb=min(lb,x-r);

        rb=max(rb,x+r);

        if(cmp(C[i+].first-x-fabs(C[i+].second-r))>)//内含

            add(C[i],C[i+]);

    }

    printf("%.2lf\n",*asr(lb,rb,eps,simpson(lb,rb)));

    return ;

}