// Calculation Program for calculating the coordinates data expressing Relationship between the Ratio of the minor to the major axes and the values of Volume and Surface area of egg ahaped body, 2011/07/21

// Equation of the egg shaped curve: (x*x+y*y)**2=a*(x**3)+(a-b)*x*(y**2)

// 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;

// Setting of the constants
	pi=3.1415927;
        a=1.0;
	l1=a;// the major axis
	 
// Setting of the other parameters
	bmax=a;// maximum value of b
	db=bmax/2000;// plotting interval of b

	xmax=a;// maximum value of x
	dx=xmax/2000;// plotting interval of x

// execution of calculation
	i=0;
	for(b=0;b<=bmax;b=b+db)
	{
		i++;

// calculation of volume of egg body
		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));
		}

 // calculation of surface area of egg body
		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);// including the equatopn of a circle with b=O
				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;
			}
		}

// calculation of the minor axis
		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);// including the equatopn of a circle with b=O
			if(y< ypast)
			{
				q=ypast;
				break;
			}
			ypast=y;
		}
		l2=2*q;// the minor axis

// summary of data
		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;

// writing the calculated coordinates data of the curve into a textfile
	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