【计算几何】【NOI2005】月下柠檬树

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

为了便于计算,假设月亮的光线是平行光,且地面是水平的,在计算时忽略树干所产生的影子。李哲当然会算了,但是他希望你也来练练手。
【输入格式】
从文件 lemon.in 中读入数据。
文件的第 1 行包含一个整数 n 和一个实数 alpha,表示柠檬树的层数和月亮的光线与地面夹角(单位为弧度)。
第 2 行包含 n+1 个实数 h0,h1,h2,...,hn,表示树离地的高度和每层的高度。
第 3 行包含 n 个实数 r1,r2,...,rn,表示柠檬树每层下底面的圆的半径。
上述输入文件中的数据,同一行相邻的两个数之间用一个空格分隔。
输入的所有实数的小数点后可能包含 1 至 10 位有效数字。
【输出格式】
将你的结果输出到文件 lemon.out 中。
输出 1 个实数,表示树影的面积。四舍五入保留两位小数。
【输入样例】
2 0.7853981633
10.0 10.00 10.00
4.00 5.00
【输出样例】
171.97
【数据范围】
1≤n≤500,0.3<alpha<π/2,0<hi≤100,0<ri≤100。
10%的数据中,n=1。
30%的数据中,n≤2。
60%的数据中,n≤20。
100%的数据中,n≤500。

这道题用simpson公式可以做，只是一开始把题理解错误，把需要计算的区域弄错了所以“误差”很大。
将整个立体图形投影到平面，则可以得到一系列圆（平行投影且被投影图形与投影面平行不改变图形的形状和大小）和等腰梯形。
（样例如图。）

而这些等腰梯形都是与圆相切的（最开始以为等腰梯形的底边就是各个圆的直径，过不了样例，还以为是simpson公式的误差太大……）
代码：

/*****************************\
 * @prob: NOI2005 lemon      *
 * @auth: Wang Junji         *
 * @stat: Accepted.          *
 * @date: June. 6th, 2012    *
 * @memo: simpson自适应公式   *
\*****************************/
#include <cstdio>
#include <cstdlib>
#include <algorithm>
#include <cstring>
#include <string>
#include <cmath>

const int maxN = 510;
const double zero = 1e-10, PI = atan(1.l) * 4, INF = 1e198;
int cnt_c, cnt_t, n;
double x[maxN], r[maxN], _x1[maxN], _x2[maxN], _y1[maxN], _y2[maxN], alpha;

inline void gmax(double &a, double b) {if (b - a > zero) a = b; return;}
inline void gmin(double &a, double b) {if (b - a < -zero) a = b; return;}
template <typename _Tp> inline _Tp sqr(const _Tp &x) {return x * x;}

inline double simpson(double L, double R, double fL, double fM, double fR)
{return (R - L) * (fL + 4 * fM + fR) / 6;}

inline double f(double X)
{
    double len = 0;
    for (int i = 0; i < n; ++i)
    {
        if (fabs(x[i] - X) < r[i])
            gmax(len, sqrt(sqr(r[i]) - sqr(X - x[i])));
        if (x[i + 1] - x[i] - fabs(r[i] - r[i + 1]) > zero
            && _x1[i] - X < -zero && X - _x2[i] < -zero)
            gmax(len, (_y2[i] - _y1[i]) / (_x2[i] - _x1[i])
                    * (X - _x1[i]) + _y1[i]);
    }
    return len;
}

double calc(double L, double Mid, double R,
            double fL, double fM, double fR, double pre)
{
    double LM = (L + Mid) / 2, RM = (Mid + R) / 2,
        fLM = f(LM), fRM = f(RM),
        sL = simpson(L, Mid, fL, fLM, fM),
        sR = simpson(Mid, R, fM, fRM, fR);
    if (fabs(sL + sR - pre) < zero) return pre;
    else return calc(L, LM, Mid, fL, fLM, fM, sL) +
                calc(Mid, RM, R, fM, fRM, fR, sR);
}

inline double get_area(double L, double R)
{
    double Mid = (L + R) / 2, fL = f(L), fR = f(R), fM = f(Mid),
        pre = simpson(L, R, fL, fM, fR);
    return calc(L, Mid, R, fL, fM, fR, pre);
}

int main()
{
    freopen("lemon.in", "r", stdin);
    freopen("lemon.out", "w", stdout);
    scanf("%d%lf", &n, &alpha); alpha = 1 / tan(alpha);
    for (int i = 0; i < n + 1; ++i)
        scanf("%lf", x + i), x[i] *= alpha;
    for (int i = 0; i < n; ++i) scanf("%lf", r + i);
    for (int i = 1; i < n + 1; ++i) x[i] += x[i - 1];
    double L = x[n], R = x[n];
    for (int i = 0; i < n; ++i)
    {
        gmin(L, x[i] - r[i]), gmax(R, x[i] + r[i]);
        if (x[i + 1] - x[i] - fabs(r[i] - r[i + 1]) > zero)
        {
            double tmp = (r[i + 1] - r[i]) / (x[i + 1] - x[i]);
            _x1[i] = x[i] - r[i] * tmp;
            _y1[i] = sqrt(sqr(r[i]) - sqr(_x1[i] - x[i]));
            _x2[i] = x[i + 1] - r[i + 1] * tmp;
            _y2[i] = sqrt(sqr(r[i + 1]) - sqr(_x2[i] - x[i + 1]));
        }
    }
    printf("%.2lf\n", 2 * get_area(L, R)); return 0;
}