hdu 5954 -- Do not pour out(积分+二分)

题目链接

 

Problem Description

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的水，求将水杯倾斜到最大角度（即水到杯口时，水不能流出）时，水面的面积？

 

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

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

 

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

         

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

 

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

 

      

         上图即是水面经过杯底的样子（我画的不是水杯倾斜，直接让水倾斜了，凑合看吧）。

         

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

 

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

 

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

 

       

         

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

         

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

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

     

       面积终于求出来了，太惊喜了！

       那么现在就该上代码了：

 

#include <iostream>
#include <algorithm>
#include <cstdio>
#include <cmath>
using namespace std;
const double PI=acos(-1);
const double eps=1e-12;

double cal2(double x)
{
    double sum=x+sin(2*x)/2;
    return sum;
}
double area(double a,double x)
{
    double sum=cal2(0.5*PI);
    sum=sum-cal2(asin(x/a));
    return sum*a;
}
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(acos(1))-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.0,r=2.0;
        for(int i=0;i<50;i++)
        {
            double mid=(l+r)/2;
            //double du=acos(1-mid);
            //double S=du-sin(du)*(1-mid);
            double V=getV(mid);
            if(fabs(V-PI*d)<eps) break;
            if(V>PI*d) r=mid;
            else l=mid;
        }
        double mid=l;
        if(mid==0.0) mid=eps;
        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;
}