// ［このプログラムの目的］カージオイドを用いて花の形を作る_i，2012年7月22日（日）



// file name: flower_i.c



#include< stdio.h>

#include< math.h>



void main(void)

{

	double a,pi;// pi は円周率

	double r,t;// 元のカージオイドの動径と位相角

	double f,df;// ハート曲線位相角、および、その変分

	double rr,ff;// 1個の花びらの動径と位相角

	double fmin,fmax;// f 最小値と最大値

	double x,y;// ハート曲線（花びらとして使う）の直交座標

	int n;// 花びらの数

	double b,c;// カージオイドからハート形曲線への変換係数

	double d;// カージオイドを縦方向に伸縮する伸縮係数

	double e,p;// 花びらの土台となる中心円の半径と位相角

	double dp;// p の変分

	double beta;// 1個の花びらが中心円の周囲を占有する角度範囲

	double dbeta;// beta の変分

	double k;// ハート曲線の角度を広げる

	double l;// 1個の花びらが土台となる芯円の周囲を占有する角度に対する補正係数(ただし、0 < l )

	int nf;// t および beta の分割数

	int m;// 1個の花びらが土台となる芯円の周囲を占有する角度 beta の変分 dbeta のカウント数

	int i,imax,j;

	double xx[20001],yy[20001];// メモリ容量の上限に注意

	double t1, t2;// t の仲介となる補助変数



	FILE *fp;



//  定数設定 [注意: 定数 b, c, d, e の数値の組み合わせによっては、解析的に計算エラーが生じる可能性が含まれている。しかし、詳細は不明。]

	pi=3.14159265;

	a=1;



	b=1;// b=1

	c=0.2;// c=-0.2

	d=1.3;// d=1.3

	e=0.2;// e=0.2

	k=0.75;// k=0.75

	l=1;// l=1



	printf("花びらの数を自然数で入力. \n n=? ");

	scanf("%d",&n);

	printf("n=%d\n",n);

	printf("\n");



	beta=2*pi/n;// 1個の花びらが土台となる芯円の周囲を占有する基本的な角度

	beta=l*beta;// beta　の補正された角度



//  他のパラメータ設定

	fmin=-pi/2;

	fmax=3*pi/2;



	nf=100;



	df=(fmax-fmin)/nf;// 変換前の位相角 t のプロット間隔

	dbeta=beta/nf;// beta の 変分



	dp=2*pi/n/100;



	i=0;



// 主計算

	for(j=1;j<=n;j++) // n個の花びらを掃引

	{

		m=0;



		for(f=fmin;f< fmax;f=f+df) //  1個の花びらを通過する位相角 f の掃引

		{

			i++;

			m++;



			if(f<=-pi/2)// 元のカージオイドからハート形への位相の変換(開始)

			{

				t1=0;

			}

			else

			{

				t1=b*sqrt(f+pi/2);

			}

			if(f>=3*pi/2)

			{

				t2=0;

			}

			else

			{

				t2=b*sqrt(3*pi/2-f);

			}

			

			t=t1-t2+(1-b*sqrt(2/pi))*f+b*sqrt(pi/2);// 元のカージオイドからハート形への位相の変換(終了)



			r=a*(1-sin(t));// ハート形の動径と位相の関係



			x=r*(1-c*sin(f)*fabs(cos(f)))*cos(f);// ハート形の形成と整形、および、直交座標への表示

			y=d*r*(1+c*sin(f))*sin(f)+2*a*d*(1-c);// ハート形の形成と整形、および、ハート形の底に原点を移動して直交座標に表示



			rr=sqrt(x*x+y*y);// 上式の動径の計算



			if(x==0)// 上式の位相の計算

			{

				ff=pi/2;

			}

			else

			{

				if(x>0)

				{

					ff=asin(y/rr);

				}

				else

				{

					ff=pi-asin(y/rr);

				}

			}

			

			ff=k*(ff-pi/2)+pi/2;// ハート曲線(花びら１枚)の角度を広げる



			ff=ff+2*pi*(j-1)/n;// j 番目の花びらの位相角



			xx[i]=rr*cos(ff)+e*cos(2*pi*(j-1)/n+pi/2-beta/2+m*dbeta);

			yy[i]=rr*sin(ff)+e*sin(2*pi*(j-1)/n+pi/2-beta/2+m*dbeta);



			printf("i=%d,x=%f,y=%f\n",i,xx[i],yy[i]);

		}



		if(l<=1)

		{

			for(p=2*pi*(j-1)/n+pi/2+beta/2;p< 2*pi*j/n+pi/2-beta/2+0;p=p+dp)

			{

				i++;



				xx[i]=e*cos(p);

				yy[i]=e*sin(p);



				printf("i=%d,x=%f,y=%f\n",i,xx[i],yy[i]);

			}

		}

		else

		{

			for(p=2*pi*(j-1)/n+pi/2+beta/2;p>2*pi*j/n+pi/2-beta/2-0;p=p-dp)

			{

				i++;



				xx[i]=e*cos(p);

				yy[i]=e*sin(p);



				printf("i=%d,x=%f,y=%f\n",i,xx[i],yy[i]);

			}

		}



	}



	i++;

	xx[i]=xx[1];// 更に、始めと終わりの隙間を埋めるための保存

	yy[i]=yy[1];





	for(p=pi/2-beta/2;p<2*pi*(n-1)/n+pi/2+beta/2-0;p=p+dp)// 中心円を描く

	{

		i++;

		xx[i]=e*cos(p);

		yy[i]=e*sin(p);

	}





	imax=i;



// 計算データ[n 個の花びらをもつ花の形の(x,y)座標]をテキストファイル（flower_i.txt）に書き込む

	fp=fopen("flower_i.txt","w");

	if(fp==NULL)

	{

		printf("FILE OPEN ERROR\n");

	}

	else

	{

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

		{

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

		}

		fflush(fp);

		fclose(fp);

	}

	printf("end\n");

}// the end of the program