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

#define NMAX 50	//粒子の個数
#define TMAX 10000	//最大時間

#define LX 10	//シミュレーション領域の横x方向長さ
#define LY 10	//シミュレーション領域の縦y方向の長さ

#define h 0.001	//時間刻み幅

double x_old[NMAX],y_old[NMAX],x[NMAX],y[NMAX];		//粒子の座標
double vx[NMAX],vy[NMAX];	//粒子の速度
double fx[NMAX],fy[NMAX];	//粒子に加わる力

void force();				//粒子に加わる力を求める関数
int roundoff(double n);		//四捨五入する関数
void note(int time);		//粒子の座標をファイルに書き込む関数

int main() {

	int i;
	int time;
	double x_new,y_new;
	double k;

	FILE *energy;
	energy=fopen("energy.dat","w");		//運動エネルギーのデータ用

/*粒子の初期配置と初速度を決める*/
	for(i=0;i<NMAX;i++) {
		x_old[i]=LX*(0.5-rand()/(RAND_MAX+1.0));
		y_old[i]=LY*(0.5-rand()/(RAND_MAX+1.0));

		x[i]=x_old[i];	//初期配置状態における力を求めるためにこちらも初期化する
		y[i]=y_old[i];

		vx[i]=0.0;
		vy[i]=0.0;
	}

	note(0);
	force();	//初期配置された状態における粒子にかかる力を求める

/*最初の１ステップ目の計算*/
	for(i=0;i<NMAX;i++) {
		x[i]=x_old[i]+vx[i]*h+fx[i]*h*h/2.0;
		y[i]=y_old[i]+vy[i]*h+fy[i]*h*h/2.0;
	}

/*Verlet法を用いて、時間ステップごとに粒子の座標を求めていく*/
	for(time=1;time<=TMAX;time++) {

		k=0.0;		//運動エネルギーを0にリセット
		force();
 
		for(i=0;i<NMAX;i++) {
				x_new=2.0*x[i]-x_old[i]+fx[i]*h*h;
				y_new=2.0*y[i]-y_old[i]+fy[i]*h*h;

				vx[i]=vx[i]+fx[i]*h;				//x方向の速度計算
				vy[i]=vy[i]+fy[i]*h;				//y方向の速度計算
				k+=(vx[i]*vx[i]+vy[i]*vy[i])/2.0;	//運動エネルギーの計算

				x_old[i]=x[i];	//古い分子の位置を格納
				y_old[i]=y[i];
				x[i]=x_new;		//いま求めた分子の位置を格納
				y[i]=y_new;

				x[i]+=-roundoff(x[i]/LX)*LX ;	//x方向の周期境界条件の考慮
				y[i]+=-roundoff(y[i]/LY)*LY ;	//y方向の周期境界条件の考慮
		}

/*一定時間おきに粒子の座標を記録する*/
		if(time%1000==0) {
			printf("time=%d\n",time);	//プログラム実行中に、現在のtime値を表示
			note(time);					//粒子の座標を記録
			
		}
		fprintf(energy,"%d\t%f\n",time,k);	//エネルギー値を記録
	}

	fclose(energy);

	return 0;
}

void force()
{
	int i,j;
	double x_i,y_i,x_ij,y_ij,r_ij;
	double fx_i,fy_i,fx_ij,fy_ij,f_ij;

/*力の値を0に初期化*/
	for(i=0;i<NMAX;i++) {
		fx[i]=0.0;
		fy[i]=0.0;
	}

/*作用反作用の法則より、j>iの分子対の計算だけで済む*/
	for(i=0;i<NMAX-1;i++) {
		x_i=x[i];
		y_i=y[i];
		fx_i=fx[i];
		fy_i=fy[i];

		for(j=i+1;j<NMAX;j++) {
			x_ij=x_i-x[j];		//i番目とj番目の粒子の間の距離(x方向)
			y_ij=y_i-y[j];		//i番目とj番目の粒子の間の距離(y方向)

			x_ij=x_ij-roundoff(x_ij/LX)*LX;	//周期境界条件を考慮
			y_ij=y_ij-roundoff(y_ij/LY)*LY;	//周期境界条件を考慮

			r_ij=sqrt(x_ij*x_ij+y_ij*y_ij);	//i番目とj番目の粒子の間の距離

			f_ij=1.0/(r_ij*r_ij*r_ij);	//距離の２乗に反比例する力を求める

			fx_ij=f_ij*x_ij;	//力のx成分の計算
			fy_ij=f_ij*y_ij;	//力のy成分の計算
			fx_i=fx_i+fx_ij;	//x方向の力を加算
			fy_i=fy_i+fy_ij;	//y方向の力を加算
			fx[j]=fx[j]-fx_ij;	//作用反作用の法則を考慮
			fy[j]=fy[j]-fy_ij;	//作用反作用の法則を考慮
			
		}
		fx[i]=fx_i;
		fy[i]=fy_i;
	}
}

/*四捨五入する関数*/
int roundoff(double n)
{
	int a;

	if(n>=0) {
		a=(int)(n+0.5);
	}
	else {
		a=(int)(n-0.5);
	}

	return a;
}

/*粒子の座標をファイルに書き込む関数*/
void note(int time)
{
	int i;
	char str[256];

	sprintf(str, "time%d.dat",time);	//ファイル名に、time値を含めるために必要な処理

/*粒子の座標を記録するファイルの作成*/
	FILE *output;
	output=fopen(str,"w");

	for(i=0;i<NMAX;i++) {
		fprintf(output,"%f\t%f\n",x[i],y[i]);
	}

	fclose(output);
}