HDU 5954 Do not pour out(二分+积分)

传送门

You have got a cylindrical cup. Its bottom diameter is 2 units and its height is 2 units as well.
The height of liquid level in the cup is d (0 ≤ d ≤ 2). When you incline the cup to the maximal angle such that the liquid inside has not been poured out, what is the area of the surface of the liquid?

 

Input

The first line is the number of test cases. For each test case, a line contains a float-point number d.

 

Output

For each test case, output a line containing the area of the surface rounded to 5 decimal places.

 

Sample Input

  
  

   
   4
0
1
2
0.424413182
  
  

 

Sample Output

  
  

   
   0.00000
4.44288
3.14159
3.51241
  
  

题目大意：

有一个圆柱瓶子，高为 $2$， 底面直径为 $2$，现在瓶子里有高为 $d(0\le d\le 2)$ 的水，现在把瓶子倾斜到最大使水不露出来，求水面的面积。

解题思路：

二分底面，然后积分求面积和体积。
由高中学过的几何知识可以知道：水面相当于对于一个很长的圆柱体倾斜的用刀切开，那么这个切面就是一个完整的椭圆，当然如果不倾斜则得到特殊的椭圆——圆，如果水面经过杯底，那么水面就是一个缺少一部分的椭圆，所以我们需要分开讨论水面经过杯底和不经过杯底两种情况。
那么这两种情况的d的临界值是多少呢？ 可以发现对于水刚到杯底的时候，有水和无水的部分各占一半，所以分界点d=1；

1、水面不经过杯底$（d\ge 1）$

这种情况如上图所示，$h+\frac {(2-h)}2=d$, 所以$h=2*d-2$,那么可以求出水面这个完整椭圆的长半径$a=\sqrt\frac {2^2+(2-h)^2}2$,而椭圆的短半经是$b=1$，所以水面面积为$S=PI*a*b$.

2、水面经过杯底$（d\lt 1）$

对于上图中 $mid$ 越大则水的体积越大，那么我们可以根据体积二分$mid$ 求出 $mid$ 真实长，最后根据真实的 $mid$ 求出水面的面积，可以知道二分范围为$(0,2)$。每次我们需要根据当前 $mid$ 求出水的体积，因为水体不规则所以必须积分求水的体积。
积分：我们根据水的高度积分，如上图所示利用 $y$ 积分求体积，那么我们需要根据 $y$ 求出每个水体截面的长 $t$ （类似于杯底的 $mid$），求相似三角形 $\frac t {mid}=\frac {(2-y)} 2$ 得$t=\frac{(2-y)*mid}2$ ; 然后根据水截面长 $t$ 求出当前水体截面的面积 $S$；
可以知道水体截面为一个扇形减去一个三角形组成，面积及体积如下图所示：

求出真实的 $mid$ 以后，那么就可以求出水面的面积了。

如上图所示利用二分求出的 $mid$ 得 $len=\sqrt{2^2+mid^2}$，设水面与杯底的一个交点为 $(x,h)$，$h=\sqrt{1-(1-mid)^2}$ (由杯底水面交线所在圆很容易求出 h ),那么 $len=a-x$ 且$\frac {x^2}{a^2}+\frac {h^2}{b^2}=1$,解上面两个方程 求得$a=\frac {len}{1+flag*\sqrt{1-h^2}}（flag表示正负），x=a-len$ 可以发现 $mid\lt 1$ 时水面为小半个椭圆，这时 $x\gt 0$,所以 $flag\lt 0$, 反之$mid\lt 1$时，水面为大半个椭圆，这时 $x\lt 0,flag\gt 0$

现在椭圆方程已求出，$X$ 范围 $(x,a)$ （不是完整的椭圆），所以需要积分，如下图所述：

代码：

#include <iostream>
#include <string.h>
#include <string>
#include <algorithm>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <map>
using namespace std;
typedef long long LL;
const int MAXN = 1e6+5;
const double PI = acos(-1);
const double eps = 1e-8;
const LL MOD = 1e9+7;

double cal2(double a,double x)
{
    double sum=x+sin(2*x)/2;
    return sum*a*0.5;
}
double area(double a,double x)
{
    double sum=cal2(a,0.5*PI);
    sum=sum-cal2(a,asin(x/a));
    return sum*2;
}
double cal(double x)
{
    double ans=sin(x)-sin(x)*sin(x)*sin(x)/3-x*cos(x);
    return ans;
}
double getV(double mid)
{
    double V=cal(0)-cal(acos(1-mid));
    V=V*(-2)/mid;
    return V;
}
int main()
{
    int T; cin>>T;
    while(T--)
    {
        double d; scanf("%lf",&d);
        if(d>1)
        {
            double h=2*d-2;
            double a=sqrt(4+(2-h)*(2-h))/2;
            printf("%.5f\n",PI*a);
            continue;
        }
        double l=0,r=2.0;
        for(int i=0;i<100;i++)
        {
            double mid=(l+r)/2;
            double V=getV(mid);
            if(V-PI*d>eps) r=mid;
            else l=mid;
        }
        double mid=l;
        int flag=1;
        if(mid<1) flag=-1;
        double len=sqrt(mid*mid+4);
        double h=sqrt(1-(1-mid)*(1-mid));
        double a=len/(1+flag*sqrt(1-h*h));
        double x=a-len;
        double res=area(a,x);
        printf("%.5f\n",res);
    }
    return 0;
}