// ［このプログラムの目的］卵の表面積を計算，2009年12月08日（火）

// 曲線の方程式：(x*x+y*y)**2=a*(x**3)+(a-b)*x*(y**2),
// b/a=0.8 の場合が最も実際の卵に近い。
// なお、一般的条件は b<=a で、b=0 のときは円となる。

// file name: surface.c

#include< stdio.h>
#include< math.h>

void main(void)
{
	double a,b,x,dx,y,xmax,s,ds,dydx,r1,r2,r3,r4,p;

//  定数設定
	for(;;)
	{
		printf("定数 a の入力; a= ? ");
		scanf("%lf",&a);
		if(a>0)
		{
			break;
		}
		else
		{
printf("a を正の値にして再入力してください。\n");
		}
	}
	for(;;)
	{
		printf("定数 b の入力; b= ? ");
		scanf("%lf",&b);
		if(b>=0&&b<=a)
		{
			break;
		}
		else
		{
printf("b を 0 =< b <= a の値で再入力してください。\n");
		}
	}

	printf("\n");

	p=3.1415927;
	xmax=a;// x の最大値
	dx=xmax/1000000;

//  計算実行
	if(b==0)
	{
		s=p*a*a;
	}
	else
	{
		for(x=dx;x<=xmax;x=x+dx)
		{
y=sqrt(((a-b-2*x)+sqrt(4*b*x+(a-b)*(a-b)))*x/2);
// b=O のときの円の方程式 y=sqrt(a*a/2/2-(x-a/2)*(x-a/2)) を含む.

r1=sqrt((a-b-2.0*x+sqrt(4.0*b*x+(a-b)*(a-b)))/x);
r2=b/sqrt(4.0*b*x+(a-b)*(a-b))-1;
r3=sqrt(a-b-2.0*x+sqrt(4.0*b*x+(a-b)*(a-b)));
r4=r2/r3;

dydx=(1/(2.0*sqrt(2.0)))*r1+sqrt(x/2.0))*r2/r3;
ds=2*p*y*sqrt(1+dydx*dydx)*dx;
s=s+ds;
		}
	}

	printf("この卵の表面積は; S=%lf\n",s);
	printf("\n");

	printf("end\n");
}// the end of the program