/* theta = k*phi + theta0 */

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

/* parameter */
#define Pi        4*atan(1.0)
//#define theta0    0.0//Pi/2.0
//#define k         1
#define x_max     5.0
#define y_max     5.0
#define bunkatsu  10.0
#define xkizami   x_max/bunkatsu
#define ykizami   y_max/bunkatsu
#define hosei     4.0

double delta_x(double r,double x,double y,int k, double theta0)
{
	double dx,phi;

	if(y < 0.0)
		phi = 2.0*Pi - acos(x/r);
	else
		phi = acos(x/r);

	if(r == 0.0)
		dx = 0.0;
	else
		dx = (cos(k*phi+theta0));

	return dx;
}

double delta_y(double r,double x,double y,int k, double theta0)
{
	double dy,phi;
  
	//phi = asin(y/r);
	if(y < 0.0)
		phi = 2.0*Pi - acos(x/r);
	else
		phi = acos(x/r);
	if(r == 0.0)
		dy = 0.0;
	else
		dy = (sin(k*phi+theta0));

	return dy;
}



int main()
{
	FILE *fp1;
	int i,j;
	double x,y,dx,dy,r;
	double theta0;
	int k;
	
	printf("k = ");
	scanf("%d",&k);
	printf("theta0 =");
	scanf("%lf",&theta0);
	theta0 = theta0*Pi;
    
	fp1 = fopen("uzu_data.dat","w");
	if(fp1 == NULL){
		printf("failure\n");
		return 0;
	}
	
	fprintf(fp1,"# x     y      x+dx     y+dy  \n");
  
	for(i=-bunkatsu;i<=bunkatsu;i++){
		x = i * xkizami;
		for(j=-bunkatsu;j<=bunkatsu;j++){
			y =  j * ykizami;
			r = sqrt(x*x + y*y);
			dx = delta_x(r,x,y,k,theta0);
			dy = delta_y(r,x,y,k,theta0);
			fprintf(fp1,"%5.3f %5.3f %8.4f %8.4f \n",x-dx/(2.0*hosei),y-dy/(2.0*hosei),dx/hosei,dy/hosei);      
		}
	}
  
	return 0;
}