// ［このプログラムの目的］５種の拡大系卵形曲線を描くこと，2008年8月13日（水）

// 曲線の方程式：x*x*x*x+2*d*x*x*y*y+f*y*y*y*y=a*x*x*x+c*x*y*y+e*y*y, 以下、a=4, d=1, f=1, c=a-b=0.8 の場合、なお、b=0 のときは円となる。
// ただし、d は y 方向の伸縮率の二乗、f も y 方向の伸縮率の四乗だが d から独立させて新定数とした。
// また、c は元々 a-b だが、y 方向の伸縮率の二乗を含めて新定数とした。

// 定数 c が可変

// file name: extended_egg_curve_2a.c

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

void main(void)
{
	double x,y,a,c,d,e,f,dc,cmin,cmax,u1,u2,u3;
	int i,imax,n,nmax,j,m;
	double xmax,dx;
	double xx[2005],yy[10][1005];// 配列を格納するメモリ容量の上限に注意

	FILE *fp;

//  定数設定
	a=4;
	d=0.25;
	f=0.5;
	e=0;

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

// 他のパラメータ設定
	xmax=a;// x の最大値
	dx=xmax/400;// x のプロット間隔

	cmin=-2;
	cmax=2;
	dc=(cmax-cmin)/4;

	i=0;

// 計算実行
	for(x=0;x<=xmax;x=x+dx)
		{
		i++;
		n=0;
		xx[i]=x;

		for(c=cmin;c<=cmax;c=c+dc)
			{
			n++;
u1=4*(d*d-f)*x*x*x*x+4*(a*f-c*d)*x*x*x;
u2=(c*c-4*d*e)*x*x+2*c*e*x+e*e;
u3=sqrt(u1+u2);
y=sqrt(((c*x+e-2*d*x*x)+sqrt(u3))/(2*f));

// y=sqrt(((c*x+e-2*d*x*x)+sqrt(4*(d*d-f)*x*x*x*x
// +4*(a*f-c*d)*x*x*x+(c*c-4*d*e)*x*x+2*c*e*x+e*e))/(2*f));
// d=f=1, c=a のときの円の方程式 
// y=sqrt(a*a/2/2-(x-a/2)*(x-a/2)) を含んでいる．
			yy[n][i]=y;

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

	imax=i;
	nmax=n;
	j=0;

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

		xx[i]=xx[m];

		for(n=1;n<=nmax;++n)
			{
			yy[n][i]=-yy[n][m];
			}
		}

// 計算データのテキストファイルへの書き込み
	fp=fopen("extended_egg_curve.txt","w");
	if(fp==NULL)
		{
		printf("FILE OPEN ERROR\n");
		}
	else
		{
		for(i=1;i<=2*imax;i++)
			{
fprintf(fp,"%f,%f,%f,%f,%f,%f\n",xx[i],yy[1][i],yy[2][i],yy[3][i],yy[4][i],yy[5][i]);
			}
		fflush(fp);
		fclose(fp);
		}
	printf("end\n");
}// the end of the program