// ［このプログラムの目的］浅井氏による「卵型曲線」と「おむすび曲線」を描くこと，2011年5月27日（金）



// ( a x2 + b x + c )2 / ( d x + e )2 + ( y / f )2 = 1

// 推奨式は  ( ( x - h ) ( g x + 1 ) )2 / ( h ( p x + 1) )2 + ( y / f )2 = 1, 卵形曲線は　x ≧0、おむすび曲線は　x < 0

// 定数の推奨値は f=1.1, g=1.0, h=1.35, p=1.3



// file name: egg_curve_by_Asai_III.c



#include< stdio.h>

#include< math.h>



void main(void)

{

	double x,y;

	double f,g,h,p;

	int i,imax,j,m;

	double x1,x2,x3,x4,dx;

	double xx[10001],yy[10001];



	FILE *fp;



//  定数設定

	f=1.1; //  推奨値 f=1.1

    g=1.0; //  推奨値 g=1.0

	h=1.35; // 推奨値 h=1.35

	p=1.3; //  推奨値 p=1.3



	printf("g=%f\n\n",g);



//  他のパラメータ設定



	x1=-(1.0+h*p-g*h)/g;

	x2=(h*(g+p)-1.0-sqrt((1.0-h*(g+p))*(1.0-h*(g+p))+8.0*g*h))/(2*g);

	x3=0;

	x4=(h*(g+p)-1.0+sqrt((1.0-h*(g+p))*(1.0-h*(g+p))+8.0*g*h))/(2*g);



	dx=(x2-x1)/4000;// x のプロット間隔



//  卵形曲線の計算実行

	i=0;

	for(x=x1+dx;x<=x2;x=x+dx)

		{

		i++;



		xx[i]=x;



		y=f*sqrt(1.0-((x-h)*(g*x+1.0)/(h*(p*x+1.0)))*((x-h)*(g*x+1.0)/(h*(p*x+1.0))));

		yy[i]=y;

		printf("i=%d,x=%f,y=%f\n",i,x,y);

		}

	imax=i;

	j=0;



	for(i=imax+1;i<=2*imax;i++)

		{

		j++;

		m=imax-j;



		xx[i]=xx[m];

		yy[i]=-yy[m];

		}



// 卵形曲線の x-y 座標デーをテキストファイル egg-shaped curve.txt に保存

	fp=fopen("egg-shaped curve.txt","w");

	if(fp==NULL)

		{

		printf("FILE OPEN ERROR\n");

		}

	else

		{

		for(i=1;i<=2*imax-1;i++)

			{

			fprintf(fp,"%f,%f\n",xx[i],yy[i]);

			}

		fflush(fp);

		fclose(fp);

		}

	printf("end\n");



//  x のプロット間隔 dx の再設定

	dx=(x4-x3)/4000;



//  付随する曲線の計算実行

	i=0;

	for(x=x3+dx;x<=x4;x=x+dx)

		{

		i++;



		xx[i]=x;



		y=f*sqrt(1.0-((x-h)*(g*x+1.0)/(h*(p*x+1.0)))*((x-h)*(g*x+1.0)/(h*(p*x+1.0))));

		yy[i]=y;

		printf("i=%d,x=%f,y=%f\n",i,x,y);

		}

	imax=i;

	j=0;



	for(i=imax+1;i<=2*imax;i++)

		{

		j++;

		m=imax-j;



		xx[i]=xx[m];

		yy[i]=-yy[m];

		}



// 付随する曲線の x-y 座標データをテキストファイル egg-shaped curve_2.txt に保存

	fp=fopen("egg-shaped curve_2.txt","w");

	if(fp==NULL)

		{

		printf("FILE OPEN ERROR\n");

		}

	else

		{

		for(i=1;i<=2*imax-1;i++)

			{

			fprintf(fp,"%f,%f\n",xx[i],yy[i]);

			}

		fflush(fp);

		fclose(fp);

		}

	printf("end\n");

}// the end of the program