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?

[align=left]Input[/align] The first line is the number of test cases. For each test case, a line contains a float-point number d.
[align=left]Output[/align] For each test case, output a line containing the area of the surface rounded to 5 decimal places.
[align=left]Sample Input[/align]

4
0
1
2
0.424413182

[align=left]Sample Output[/align]

0.00000
4.44288
3.14159
3.51241

题目大意：

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

解题思路：

二分底面，然后积分求面积和体积。

由高中学过的几何知识可以知道：水面相当于对于一个很长的圆柱体倾斜的用刀切开，那么这个切面就是一个完整的椭圆，当然如果不倾斜则得到特殊的椭圆——圆，如果水面经过杯底，那么水面就是一个缺少一部分的椭圆，所以我们需要分开讨论水面经过杯底和不经过杯底两种情况。

那么这两种情况的d的临界值是多少呢？ 可以发现对于水刚到杯底的时候，有水和无水的部分各占一半，所以分界点d=1；

1、水面不经过杯底（d≥1）



这种情况如上图所示，h+(2−h)2=d, 所以h=2∗d−2,那么可以求出水面这个完整椭圆的长半径a=22+(2−h)22−−−−−−−√,而椭圆的短半经是b=1，所以水面面积为S=PI∗a∗b.

2、水面经过杯底（d<1）



对于上图中 mid 越大则水的体积越大，那么我们可以根据体积二分mid 求出 mid 真实长，最后根据真实的 mid 求出水面的面积，可以知道二分范围为(0,2)。每次我们需要根据当前 mid 求出水的体积，因为水体不规则所以必须积分求水的体积。

积分：我们根据水的高度积分，如上图所示利用 y 积分求体积，那么我们需要根据 y 求出每个水体截面的长 t （类似于杯底的 mid），求相似三角形 tmid=(2−y)2 得t=(2−y)∗mid2 ; 然后根据水截面长 t 求出当前水体截面的面积 S；

可以知道水体截面为一个扇形减去一个三角形组成，面积及体积如下图所示：





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



如上图所示利用二分求出的 mid 得 len=22+mid2−−−−−−−−√，设水面与杯底的一个交点为 (x,h)，h=1−(1−mid)2−−−−−−−−−−−−√ (由杯底水面交线所在圆很容易求出 h ),那么 len=a−x 且x2a2+h2b2=1,解上面两个方程 求得a=len1+flag∗1−h2√（flag表示正负），x=a−len 可以发现 mid<1 时水面为小半个椭圆，这时 x>0,所以 flag<0, 反之mid<1时，水面为大半个椭圆，这时 x<0,flag>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;
}