// ［このプログラムの目的］卵の体積と表面積の 短径/長径 依存性のグラフを描くこと，2011年7月21日（木）

// 曲線の方程式：(x*x+y*y)**2=a*(x**3)+(a-b)*x*(y**2), 以下、a=1　として, b の値を 0=< b<=a に変えてグラフを計算する。

// file name: volume_&_surface.c

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

void main(void)
{
	double v,s,ds,a,b,q;
	double l1,l2;
	double pi;
	int i,imax;
	double bmax,db;
	double x,y,ypast;
	double xmax,dx,dydx;
	double l12[5001],vv[5001],ss[5001];

	FILE *fp;

//  定数設定
	pi=3.1415927;
    a=1.0;
	l1=a;// 長径
	 
//  他のパラメータ設定
	bmax=a;// b の最大値
	db=bmax/2000;// b のプロット間隔

	xmax=a;// x の最大値
	dx=xmax/2000;// 卵型の横軸座標 x のプロット間隔

//  計算実行
	i=0;
	for(b=0;b<=bmax;b=b+db)
	{
		i++;

//  体積計算
		if(b==0)
		{
			v=pi*a*a*a/6;
		}
		else
		{
			v=(pi/2)*((a+b)*(a+b)*(a+b)*a/(6*b)-a*a*a/6-a*a*b/2-((a+b)*(a+b)*(a+b)*(a+b)*(a+b)-(a-b)*(a-b)*(a-b)*(a-b)*(a-b))/(60*b*b));
		}

 //  表面積計算
		if(b==0)
		{
			s=pi*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)) を含んでいる．
				dydx=(1/(2.0*sqrt(2.0)))*sqrt((a-b-2.0*x+sqrt(4.0*b*x+(a-b)*(a-b)))/x)+(sqrt(x/2.0))*((b/sqrt(4.0*b*x+(a-b)*(a-b))-1)/sqrt(a-b-2.0*x+sqrt(4.0*b*x+(a-b)*(a-b))));
				ds=2*pi*y*sqrt(1+dydx*dydx)*dx;
				s=s+ds;
			}
		}

//  短径計算
		ypast=-1;
		for(x=0;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)) を含んでいる．
			if(y< ypast)
			{
				q=ypast;
				break;
			}
			ypast=y;
		}
		l2=2*q;// 短径

//  データのまとめ
		vv[i]=v;
		ss[i]=s;
		l12[i]=l2/l1;

		printf("i=%d,l12=%f,v=%f,s=%f\n",i,l12[i],v,s);

		s=0;
	}
	imax=i;

// 計算データのテキストファイルへの書き込み
	fp=fopen("volume_&_surface.txt","w");
	if(fp==NULL)
	{
		printf("FILE OPEN ERROR\n");
	}
	else
	{
		for(i=1;i<=imax;i++)
		{
			fprintf(fp,"%f,%f,%f\n",l12[i],vv[i],ss[i]);
		}
		fflush(fp);
		fclose(fp);
	}
	printf("end\n");
}// the end of the programm