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

#define IADD   453806245
#define IMUL   314159269
#define MASK   2147483647
#define SCALE  0.4656612873e-9

#define DIM 2

typedef struct {
	double x, y;
} Vect2D;

typedef struct {
	Vect2D r, v, a;
} Mol;

void SingleStep(void);
void LeapfrogStep(int);
void ApplyBoundaryCond(void);
void ComputeForces(void);
void InitVels(void);
void InitCoords(void);
void VRand(Vect2D *);
double RandR(void);
void InitRand(int);

int randSeedP = 17;
const int Nx = 20;
const int Ny = 20;
const int N = Nx*Ny;

Mol *mol;
Vect2D region; //obaszr symulacji, prostokąt Lx x Ly; dalej jest Lx = Ly
double dt = 0.005, time_now, density = 0.82, temperature = 1;
double dr_cut = pow(2, 1.0/6.0), dr2_cut = pow(dr_cut, 2); // pomocnicze do obcinania zasięgu sił
double velMag = sqrt (2 * (1.0 - 1.0 / N) * temperature); // moduł prędkości początkowych
double uSum = 0, virSum = 0, virSumStep = 0, uSumStep = 0;
int more_cycles, step_count, step_equil, step_limit = 10000;
//int Nx = (int)pow(N, 0.5), Ny = Nx;

int main() {
	int step_count = 0, i, j, k, k_avg;
	double vx2_sum = 0, vy2_sum = 0;
	
	/* alokacja tablicy mol[N], to jest nasz układ;
	   komórka i-ta, Mol[i] zawiera r, v, a dla molekuły i-tej */
	mol = (Mol*) malloc(N*sizeof(Mol));
	region.x = Nx/pow(density, 0.5); // Lx
	region.y = Ny/pow(density, 0.5); // Ly
		
	for (i=0; i < N; i++) {
		mol[i].r.x = mol[i].r.y = 0.0;
		mol[i].v.x = mol[i].v.y = 0.0;
		mol[i].a.x = mol[i].a.y = 0.0;
	}
	
	InitCoords();
	InitVels();
		
	more_cycles = 1; k = 1; k_avg = 100;
	while (more_cycles) {
		step_count++;
		SingleStep();
		
		// Show a sample trajectory
		if (step_count % k_avg == 0) {
			printf("%6.4f\t%6.4f\t", mol[315].r.x, mol[315].r.y);
			printf("%6.4f\t%6.4f\n", mol[317].r.x, mol[317].r.y);
		}
		/* Thermodynamics properties: kinetic, potenital, total energy, pressure
		for (i=0; i < N; i++) { vx2_sum += pow(mol[i].v.x, 2); vy2_sum += pow(mol[i].v.y, 2); }
		uSum += uSumStep;
		virSum += virSumStep;
		if (step_count % k_avg == 0) {
			printf("step: %d, Ekin = %6.4f ", step_count, (0.5/(N*k_avg))*(vx2_sum + vy2_sum));
			printf("Epot = %6.4f, Etot = %6.4f ", uSum/(k_avg*N), uSum/(k_avg*N) + (0.5/(N*k_avg))*(vx2_sum + vy2_sum));
			printf("P = %6.4f\n", density*(vx2_sum + vy2_sum + virSum)/(DIM*N*k_avg));
			vx2_sum = vy2_sum = 0.0;
			uSum = 0.0;
			virSum = 0.0;
		}*/

		if (step_count >= step_limit) more_cycles = 0;
	}
	
	return 0;
}

void SingleStep() {
  time_now = step_count * dt;
  LeapfrogStep(1);
  ApplyBoundaryCond();
  ComputeForces();
  LeapfrogStep(2);
}

void LeapfrogStep(int part) {
	int i;
	
	if (part == 1)
		for (i=0; i < N; i++) {
			mol[i].v.x += 0.5*dt*mol[i].a.x;
			mol[i].v.y += 0.5*dt*mol[i].a.y;
			mol[i].r.x += dt*mol[i].v.x;
			mol[i].r.y += dt*mol[i].v.y;
		}
	else if (part == 2)
		for (i=0; i < N; i++) {
			mol[i].v.x += 0.5*dt*mol[i].a.x;
			mol[i].v.y += 0.5*dt*mol[i].a.y;
		}
}

void ApplyBoundaryCond() {
	int i;
	for (i=0; i < N; i++) {
		if (mol[i].r.x >= 0.5*region.x) mol[i].r.x -= region.x;
		else if (mol[i].r.x <= -0.5*region.x) mol[i].r.x += region.x;
		
		if (mol[i].r.y >= 0.5*region.y) mol[i].r.y -= region.y;
		else if (mol[i].r.y <= -0.5*region.y) mol[i].r.y += region.y;
	}
}

void ComputeForces() {
	double drx, dry, dr, dr2, dr6, fc_val;
	int i, j;
	
	uSumStep = 0.0; virSumStep = 0.0;
	for (i=0; i < N; i++) { mol[i].a.x = 0; mol[i].a.y = 0; }
	
	for (i=0; i < N; i++) {
		for (j=0; j < N; j++)
			if (i != j) {
				drx = mol[i].r.x - mol[j].r.x;
				dry = mol[i].r.y - mol[j].r.y;
				
				if (drx >= 0.5*region.x) drx -= region.x;
				else if (drx < -0.5*region.x) drx += region.x;
				if (dry >= 0.5*region.y) dry -= region.y;
				else if (dry < -0.5*region.y) dry += region.y;
				
				dr2 = drx*drx + dry*dry;
				if (dr2 < dr2_cut) {
					dr6 = pow(dr2, 3);
					fc_val = 48.0*(1.0/dr6)*(1.0/dr6 - 0.5)*(1.0/dr2);
					mol[i].a.x += fc_val*drx;
					mol[i].a.y += fc_val*dry;
							
					if (i < j) {
						uSumStep += 4*(1.0/dr6)*(1.0/dr6 - 1) + 1;
						virSumStep += fc_val*dr2;
					}
				}
			}
	}
}

void InitCoords() {
	double Lx, Ly, dx, dy;
	int i, j;
	
	Lx = region.x; Ly = region.y;
	dx = Lx/Nx; dy = Ly/Ny;
	
	for (i=0; i < Nx; i++)
		for (j=0; j < Ny; j++) {
			mol[Ny*i + j].r.x = dx*(i + 0.5) - 0.5*Lx;
			mol[Ny*i + j].r.y = dy*(j + 0.5) - 0.5*Ly;
		}
	
}

void InitVels() {
	int i;
	Vect2D vSum;

	for (i=0; i < N; i++) { vSum.x = vSum.y = 0; }
	for (i=0; i < N; i++) {
    	VRand (&mol[i].v);
    	mol[i].v.x *= velMag;
    	mol[i].v.y *= velMag;
    	vSum.x += mol[i].v.x;
    	vSum.y += mol[i].v.y;
	}
	
	// korekta zapewniająca, że środek masy układu w chwili t = 0 ma prędkość zerową
	for (i=0; i < N; i++) {
		mol[i].v.x += (-1.0/N) * vSum.x;
  		mol[i].v.y += (-1.0/N) * vSum.y;
	}
}

void VRand(Vect2D *p) {
	double s;

	s = 2. * M_PI * RandR();
	p->x = cos(s);
	p->y = sin(s);
}

double RandR() {
	randSeedP = (randSeedP * IMUL + IADD) & MASK;
	return (randSeedP * SCALE);
}

void InitRand(int randSeedI) {
	struct timeval tv;

	if (randSeedI != 0) randSeedP = randSeedI;
	else {
		gettimeofday (&tv, 0);
		randSeedP = tv.tv_usec;
	}
}