Luogu4207 月下柠檬树

题面

简化题意

一棵合法的柠檬树如图，其中除最顶层外，每一层都是一个圆台，最顶层是一个圆锥：

给定 $h_i,r_i,\alpha$ ，求其在月光下投影的面积，保留两位小数。

$n\leqslant 500,0.3<\alpha<\frac{2}{\pi}$ 。

题解

单独的圆锥

这同时也是 subtask1 这个部分分的条件。考虑一个圆锥的投影：

根据投影知识，不难得到， $AJ=\frac{AD}{\tan\alpha}$ 。

我们可以直接得到 $J$ 的坐标，投影的面积可以用 $\triangle JKM$ 的面积和扇形 $KCG$ 的面积的和的两倍表示。所以唯一的问题就是：已知 $J,C$ 坐标， $\odot C$ 的半径，求 $S_{\triangle JKC},\angle KCG$ 。因为 $JC,r$ 已知，不难得到 $\angle KCJ$ 和它的三角函数，进而解决问题，不多赘述。

单独的圆台

考虑单独的一个圆台怎么做（geogebra 画不出圆台，只能用两个圆简陋表示一下）：

显然， $E,F$ 坐标都能求出来，且 $\odot E$ 的半径 $r_1$ ， $\odot F$ 的半径 $r_2$ 都知道。为了方便，画出俯视图：

设两条外公切线交于点 $A$ 。有：

\sin\angle EAK=\frac{FG}{AF}=\frac{EK}{AF+EF}

$FG,EK,EF$ 为已知，可以解出 $AF$ 的值，然后 $A$ 的坐标和 $\angle EAK$ 的三角函数就出来了，于是能求出公切线的解析式，于是 $G,K$ 坐标都可以求出来了。

还是只考虑它的一半，不难发现就是直角梯形 $GKEF$ 的面积加上两个扇形。都可以由已知量简单推导得到。

我们联合！

考虑这两个东西加在一起会产生什么效应：

题目瞬间感觉毒瘤了不少。我们还是只考虑一半的情况，发现这个东西是可以分割的，于是作这些切点关于 $y$ 轴的垂线：

至于梯形部分，和上面圆台的解决方法是一样的。唯一的三角形和圆锥的解决方法类似。因此重点考虑“斜边”是一段弧的东西。即图中类似于封闭图形 $DHUN$ 的部分。

法一

所求面积为：

S_{\text{扇形} FAC}-S_{\triangle AEC}-S_{\text{扇形} FAB}+S_{\triangle ADB}

$C$ 在 $A$ 下方时同理。

法二

因为圆一定在 $y$ 轴上，所以圆的解析式为 $r^2=x^2+(y-a)^2$ 。

为了方便我们把 $x,y$ 反过来，让它落在 $x$ 轴上，解析式变为 $r^2=(x-a)^2+y^2$ 。先考虑 $a=0$ 的情况，后续再把 $x$ 用 $x-a$ 代进去。

y=\sqrt{r^2-x^2}

求它的一个不定积分：

\int \sqrt{r^2-x^2}\, {\rm d}x

我们知道 $-r\leqslant x\leqslant r$ ，所以 $x$ 可以表示为 $r\sin\varphi$ 的形式：

\int \sqrt{r^2-r^2\sin^2\varphi} \, {\rm d}r\sin\varphi

然后我们发现根号里面的东西就是 $r^2\cos^2\varphi$ ：

\int r\cos\varphi \, {\rm d}r\sin\varphi = \int r^2\cos\varphi(\sin\varphi)'\, {\rm d}\varphi

$(\sin x)'=\cos x$ ， $r$ 为常数，可以提出来

r^2\int \cos^2 \varphi \, {\rm d}\varphi

诱导公式变一下，再提一波常数：

\frac{r^2}{2}\int 1+\cos 2\varphi \, {\rm d}\varphi

和的积分等于积分的和，考虑第二项怎么积分，看成 $f(x)=\cos x,g(x)=2x$ 的复合函数即可：

\frac{r^2}{2}(\varphi+\frac{\sin 2\varphi}{2})

然后把 $x$ 代回去，二倍角公式展开三角函数， $\varphi=\arcsin \frac{x}{r}$ ：

\frac{r^2}{2}\left(\arcsin \frac{x}{r}+\frac{x}{r}\times\sqrt{1-\frac{x^2}{r^2}}\right)

然后化简一下就可以了：

\frac{r^2}{2}\arcsin \frac{x}{r} + \frac{x\sqrt{r^2-x^2}}{2}

记得把 $x-a$ （ $a$ 为圆心横坐标）代到 $x$ 中而不是直接用 $x$ ：

\frac{r^2}{2}\arcsin\frac{x-a}{r}+\frac{(x-a)\sqrt{r^2-(x-a)^2}}{2}

于是求那个东西的面积的时候，假设弧的端点横坐标分别为 $x_1,x_2$ ，将其分别代入然后减一下即可。

细节

这一段我一直没有讨论全，看了 iuyi 大佬的题解后豁然开朗%%%。

一些细节问题不拿出来单独讨论会直接出错，例如两个圆之间没有公切线（包含关系）的情况：

针对这个图，我们添加一个判断：

把最上面和最下面的点也视为分段点。
如果一个段内同时有多个圆弧或线段，取它们面积的最大值为这一段的面积。

但是我们又发现上面的图有些极端，不一定存在一个大圆包含了所有其它的圆，也不一定是圆心坐标大的点包含圆心坐标小的点。最为保险的还是添加两条：

圆和圆、线段之间的交点、线段和线段的交点，也记作分段的点。
如果一个圆包含下一个圆（按圆心坐标大小枚举时），不要计算公切线。

圆和圆的交点、线段和线段的交点比较简单，圆和线段的交点直接联立解析式都能做。

还有就是注意精度误差产生的问题，有些地方不加限制条件会爆 WA（特别是代码中用注释提到的那一段）。

于是这个题就做出来了，剩下的都交给代码实现了 qwq。

代码

有借鉴于 @iuyi 的题解代码，仅供参考。

#include<bits/stdc++.h>
#define double long double
using namespace std;
const int max_n=502;
const double inf=1000000000.0;

inline double Rt(double a,double b,double c){
    if(b==inf) swap(a,b);
    if(c==inf)
        return sqrt(a*a+b*b);
    return sqrt(c*c-b*b);
}

int n,cnt;
double alpha,h[max_n],r[max_n];
double dn[max_n],up[max_n];

struct dot{
    double x,y;
}p[max_n*max_n*10];

inline dot maked(double a,double b){
    dot res;
    res.x=a,res.y=b;
    return res;
}

double sep[max_n*max_n];

inline void Circ_Circ(){ // 圆和圆的交点
    for(register int i=1;i<=n;++i)
        for(register int j=i+1;j<=n;++j)
            if(dn[i]<dn[j] && up[i]>dn[j] && up[i]<up[j])
                sep[++cnt]=h[i]-((r[j]*r[j]-r[i]*r[i])/(h[j]-h[i])-h[j]+h[i])/2;
}

struct line{
    dot a,b;
    inline bool valid(){
        if(a.x==inf || a.y==inf) return 0;
        if(b.x==inf || b.y==inf) return 0;
        return 1;
    }
    inline double slope(){
        if(a.x==b.x) return inf;
        return (a.y-b.y)/(a.x-b.x);
    }
}l[max_n];

inline line makel(dot a,dot b){
    line res;
    res.a=a,res.b=b;
    return res;
}

inline void Line_Line(){ // 线段和线段的交点，联立解析式
    for(register int i=1;i<=n;++i)
        for(register int j=i+1;j<=n;++j)
            if(l[i].valid() && l[j].valid()){
                double k1=l[i].slope(),k2=l[j].slope();
                if(k1==inf || k2==inf || k1==k2)
                    continue;
                double b1=l[i].a.y-k1*l[i].a.x,b2=l[j].a.y-k2*l[j].a.x;
                double X=(b2-b1)/(k1-k2);
                if(X>=l[i].a.x && X<=l[i].b.x && X>=l[j].a.x && X<=l[j].b.x)
                    sep[++cnt]=X;
            }
}

inline void Circ_Line(){ // 圆和线段的交点，联立解析式
    for(register int i=1;i<=n;++i)
        for(register int j=1;j<=n;++j)
            if(j!=i && j!=i-1 && l[j].valid()){
                double k=l[j].slope(),b=l[j].a.y-k*l[j].a.x;
                double A=k*k+1,B=2*k*b-2*h[i],C=b*b+h[i]*h[i]-r[i]*r[i],delta=B*B-4*A*C;
                if(delta<0)
                    continue;
                double x1=(-B+sqrt(delta))/(A*2),x2=(-B-sqrt(delta))/(A*2);
                if(x1>=l[j].a.x && x1<=l[j].b.x)
                    sep[++cnt]=x1;
                if(x2>=l[j].a.x && x2<=l[j].b.x)
                    sep[++cnt]=x2;
            }
}

double ans=0;

inline double S_Line(double x1,double x2,int i){ // 线段在区间内的面积
    x1=max(x1,l[i].a.x),
    x2=min(x2,l[i].b.x);
    if(x1>=x2)
        return 0;
    double k=l[i].slope(),b=l[i].a.y-k*l[i].a.x;
    if(k==inf)
        return 0;
    return (k*x1+b+k*x2+b)*(x2-x1)/2;
}

inline double Calc(double x,double a,double r){ // 积分
    x-=a;
    return r*r/2*asin(max((double)-1,min((double)1,x/r)))+x*sqrt(max((double)0,r*r-x*x))/2;
    // 如果没有限制条件会因为精度误差 WA
}

inline double S_Circ(double x1,double x2,int i){ // 圆在区间内的面积
    x1=max(x1,dn[i]),
    x2=min(x2,up[i]);
    if(x1>=x2)
        return 0;
    return Calc(x2,h[i],r[i])-Calc(x1,h[i],r[i]);
}

signed main(){
//  freopen(".in","r",stdin),
//  freopen(".out","w",stdout);
    scanf("%d%Lf",&n,&alpha);
    alpha=1.0/tan(alpha);
    for(register int i=1;i<=n+1;++i){
        scanf("%Lf",&h[i]);
        h[i]=h[i]*alpha+h[i-1];
    }
    sep[++cnt]=inf,
    sep[++cnt]=-inf;
    for(register int i=1;i<=n+1;++i){
        if(i<=n)
            scanf("%Lf",&r[i]);
        dn[i]=h[i]-r[i],
        up[i]=h[i]+r[i];
        sep[1]=min(sep[1],dn[i]),
        sep[2]=max(sep[2],up[i]);
    }
    Circ_Circ();
    for(register int i=1;i<=n;++i){ // 求出所有公切线
        if(up[i]>=up[i+1] || dn[i]>=dn[i+1]){
            l[i]=makel(maked(inf,inf),maked(inf,inf));
            continue;
        }
        if(r[i]==r[i+1]){
            sep[++cnt]=h[i],
            sep[++cnt]=h[i+1];
            l[i]=makel(maked(h[i],r[i]),maked(h[i+1],r[i+1]));
        }
        else if(r[i+1]==0){
            double delta=r[i]*r[i]/(h[i+1]-h[i]);
            sep[++cnt]=h[i]+delta,
            sep[++cnt]=h[i+1];
            l[i]=makel(maked(h[i]+delta,Rt(delta,inf,r[i])),maked(h[i+1],0));
        }
        else{
            double AF=r[i+1]*(h[i+1]-h[i])/(r[i]-r[i+1]);
            double Gx=h[i+1]+r[i+1]*r[i+1]/AF,Kx=h[i]+r[i]*r[i+1]/AF;
            sep[++cnt]=Kx,
            sep[++cnt]=Gx;
            l[i]=makel(maked(Kx,Rt(r[i]*r[i+1]/AF,inf,r[i])),maked(Gx,Rt(r[i+1]*r[i+1]/AF,inf,r[i+1])));
        }
    }
    Line_Line(),
    Circ_Line();
    sort(sep+1,sep+1+cnt),
    cnt=unique(sep+1,sep+1+cnt)-sep-1;
    for(register int i=1;i<cnt;++i) if(sep[i]!=sep[i+1]){
        double mx=0;
        for(register int j=1;j<=n;++j){
            mx=max(mx,S_Line(sep[i],sep[i+1],j)),
            mx=max(mx,S_Circ(sep[i],sep[i+1],j));
        }
        ans+=mx;
    }
    printf("%.2Lf",ans*2);
    return 0;
}