真是一道让人头疼的几何模拟题。。。wa了n多遍。。。注意圆台部分上底半径的求解,否则圆台体积会求错。。。横放时圆台部分体积要微元求解。。。最终s的确定要靠二分逼近。。。
#include<stdio.h>
#include<stdlib.h> #include<string.h> #include<math.h> #define pi 3.1415926535898 #define N 1000 double k,hb,db,hn,dn,h,s,v,V,r1,r2,r,delta; double area(double r,double s) //横放时弓形底面积 { double t,s1,s2,ang; //if(s>r) return (pi*(r*r)-area(r,2*r-s)); t=db/2-s; ang=t<r?acos(t/r):0; s1=ang*r*r; s2=r*cos(ang)*r*sin(ang); //s2=t*sqrt(r*r-t*t); return s1-s2; } double volume(double s) //横放液面高度为s时的体积 { if(s*2>db) return V-volume(db-s); double vb,vn,vc; vb=area(db/2,s)*hb; vn=area(dn/2,s)*hn; vc=area(db/2,s)+area(dn/2,s); int i; for(i=1;i<N;i+=2) vc+=4*area((db+(dn-db)*i/N)/2,s); for(i=2;i<N;i+=2) vc+=2*area((db+(dn-db)*i/N)/2,s); vc*=(h-hb-hn)/3/N; return vb+vc+vn; } int main() { double i,d; while(scanf("%lf%lf%lf%lf%lf%lf",&k,&hb,&db,&hn,&dn,&h),!(k==0&&hb==0&&db==0&&hn==0&&dn==0&&h==0)){ r1=db/2; r2=dn/2; //V=pi*r1*r1*hb+pi*(r1*r1+r1*r2+r2*r2)*(h-hb-hn)/3+pi*r2*r2*hn; //瓶子总体积 V = 2 * volume(r1); if(k<=hb){ v=pi*r1*r1*k; //液体体积 } else if(k<(h-hn)){ d=r2*(h-hn-hb)/(r1-r2); //wa了两天,就因为此处d和r的推导有错误 r=(d+h-k-hn)/(d+h-hb-hn)*r1; v=pi*r1*r1*hb+pi*(r1*r1+r1*r+r*r)*(k-hb)/3; } else{ v=pi*r1*r1*hb; v+=pi*(r1*r1+r1*r2+r2*r2)*(h-hb-hn)/3; v+=pi*r2*r2*(hn+k-h); } s=db/2; for(i=s/2;i>.00001;i/=2){ //二分逼近 if(volume(s)>v) s-=i; else s+=i; } printf("%.2lf\n",s); } system("pause"); return 0; }