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

//#include 

#define Pi 3.141592653

double fai(double r);
double F(double r);
void calc_force(int n_atom, double x[], double f_x[], double y[], double f_y[], double *pot_e);
void initial_cond(int n_atom, double x[], double y[], double v_x[], double v_y[], double v0);
void make_fcc(int n, double x[], double y[]);

void main(void){
	FILE *fpe, *fp, *fpxy;
	
	double x[1000], y[1000], v_x[1000], v_y[1000], f_x[1000], f_y[1000];
	double dt;
	int i, n, n_max;
	int n_atom;
	int z;
	double m, v0;     //mass
	double kin_e, pot_e;
	double constant;

        if((fpe = fopen("./energy.xls","w")) == NULL){
		printf("File not open\n");
		return;
	}
     
	if((fp = fopen("./sample.mgf","w")) == NULL){
		printf("File not open\n");
		return;
	}

	if((fpxy = fopen("./xy.xls","w")) == NULL){
		printf("File not open\n");
		return;
	}

	n_atom = 100;
	dt = 0.01;
	m = 1.0;
	v0 = 1.0;
	n_max = 1000;
	pot_e = 0;
	constant = 0;

	fprintf(fpe,"n\tpot_e + kin_e\tpot_e\tkin_e\t\n");
	fprintf(fp, "# Micro AVS Geom:2.00\n");
	fprintf(fp, "%d\n", n_max);

	initial_cond(n_atom, x, y, v_x, v_y, v0);

	calc_force(n_atom, x, f_x, y, f_y, &pot_e);

	for(n = 0; n <= n_max; n++){
		for(i = 1; i <= n_atom; i++){
			x[i] = x[i] +v_x[i]*dt + 0.5*dt*dt*f_x[i]/ m;
			y[i] = y[i] +v_y[i]*dt + 0.5*dt*dt*f_y[i]/ m;
		}

		for(i = 1; i <= n_atom; i++){
			v_x[i] = v_x[i] + dt*0.5* f_x[i] / m;
			v_y[i] = v_y[i] + dt*0.5* f_y[i] / m;
		}

		calc_force(n_atom, x, f_x, y, f_y, &pot_e);
        
		kin_e = 0;
		for(i = 1; i <= n_atom; i++){
			v_x[i] = v_x[i] + dt*0.5* f_x[i] / m;
			v_y[i] = v_y[i] + dt*0.5* f_y[i] / m;
			kin_e += pow(v_x[i], 2) + pow(v_y[i], 2);
		}
		
		kin_e *= 0.5 * m;

		constant = pot_e + kin_e;

		printf("%d : ", n);
		printf("%f", constant);
		printf("%f %f", pot_e, kin_e);
//      for(i = 1; i <= n_atom; i++) printf("%f ", x[i]);
//      for(i = 1; i <= n_atom; i++) printf("%f", v_x[i]);
        printf("\n");

              if( (n % 100) == 0){	
		fprintf(fpxy,"%d\n", n);
		for(i = 1; i <= n_atom; i++) fprintf(fpxy, "%f\t%f\n", x[i], y[i]);
        fprintf(fpxy,"\n");
              }

		fprintf(fp, "step%d\n", n + 1);
        fprintf(fp, "sphere\n");
        fprintf(fp, "sample\n");
        fprintf(fp, "color\n");
        fprintf(fp, "%d\n", n_atom);
        for(i = 1; i <= n_atom; i++) fprintf(fp,"%9f  %f  0.0  0.5   0.0  1.0  0.0\n", x[i], y[i]);
//      fprintf(fp, "\n");

		fprintf(fpe, "%d\t%f\t%f\t%f\t\n", n, constant ,pot_e, kin_e);

	}
	scanf("&z");
	fclose(fpe);fclose(fp);	fclose(fpxy);
}

double fai(double r){
	return - 2.0 * pow(r, - 6) + pow(r, - 12);
}

double F(double r){
	return 12.0 * (pow(r, - 13) - pow(r, - 7));
}

void calc_force(int n_atom, double x[], double f_x[], double y[], double f_y[], double *pot_e){
	int i, j;
	double r, fx, fy;

	*pot_e = 0;

	for(i = 1; i <= n_atom; i++){
			f_x[i] = 0;
			f_y[i] = 0;
	}

	for(i = 1; i <= n_atom - 1; i++){
		for(j = i + 1; j <= n_atom; j++){
			if(i == j)continue;
			else{
				r = sqrt(pow(x[i] - x[j], 2)+pow(y[i] - y[j], 2));

				if( r > 1.5) continue;
				
				*pot_e += fai(r);

				fx = (x[i] - x[j]) * F(r) / r;
				fy = (y[i] - y[j]) * F(r) / r;

				f_x[i] += fx;
				f_x[j] -= fx; 
				f_y[i] += fy;
				f_y[j] -= fy;
			}
		}
	}
	return;
}

void initial_cond(int n_atom, double x[], double y[], double v_x[], double v_y[], double v0){
    int modul, multi, inc, seed;
	double r, theta, v_av_x, v_av_y;
	int i;

	modul = 2147483647;
	multi = 16897;
	inc = 0;
	seed = 4711;

	make_fcc(n_atom, x, y);

	r = (double)seed;

	for(i = 1; i <= n_atom; i++){
		r=fmod( multi*r+inc, (double)modul);
	    theta = r * 2 * Pi;
	    v_x[i] = v0 * cos(theta);
	    v_y[i] = v0 * sin(theta);
	}

    v_av_x=0; v_av_y=0;
	for(i = 1; i <= n_atom; i++){
	    v_av_x += v_x[i];
	    v_av_y += v_y[i];
	}
	v_av_x /= n_atom;
	v_av_y /= n_atom;

	for(i = 1; i <= n_atom; i++){
	    v_x[i] -= v_av_x;
	    v_y[i] -= v_av_y;
	}


}
	
void make_fcc(int n_atom, double x[], double y[]){
	int i, j;
	int a = 1;
	double n;

	n = (int) sqrt((double)n_atom);

	for(i = 1; i <= n; i++){
		if(i / 2 == (i + 1) / 2){      //偶数
			for(j = 0; j < n; j++){
				x[a] = j;
				y[a] = 0.5*sqrt(3.0)*i;				
				a++;
                
			}
		}
		else{                            //奇数
			for(j = 0; j < n; j++){
				x[a] = j + 0.5;
				y[a] = 0.5*sqrt(3.0)*i;
				a++;
							
			}
		}
		
	}
}