// NBody_Simulation.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"
#include <iostream>
#include <glut.h>
#include <windows.h>
#include <stdlib.h>
#include <vector>
#include <random>
#include <cmath>
#include <vector>
#include <omp.h>
#define GC 6.67408*pow(10,-11)
#define TIMESTAMP 150
#define INTERVAL 20
#define DIMENSION 400
#define SIZE 10
#define	ITERATION 1000
#define SCALING_FACTOR 10000

using namespace std;

struct Particle {
	double x, y, mass;
	float ux, uy;
	unsigned char r, g, b, a;
};
vector<Particle> particles;
vector<Particle> particlesToDisplay;
vector<double> distTravelX;
vector<double> distTravelY;

double distanceX(Particle p1, Particle p2) {
	return abs(p2.x - p1.x);
}
double distanceY(Particle p1, Particle p2) {
	return abs(p2.y - p1.y);
}

int ROUND = 0;
void timer(int value){
	for (size_t i = 0; i < particlesToDisplay.size()&& (ROUND*particlesToDisplay.size())<(ITERATION*SIZE); i++){
		particlesToDisplay[i].x += distTravelX[i+(ROUND*particlesToDisplay.size())] * SCALING_FACTOR;
		particlesToDisplay[i].y += distTravelY[i + (ROUND*particlesToDisplay.size())] * SCALING_FACTOR;
	}
	++ROUND;
	glutPostRedisplay();
	glutTimerFunc(INTERVAL, timer, 1);
}

void display(void){
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(-DIMENSION, DIMENSION, -DIMENSION, DIMENSION, -1, 1);

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();

	glColor3ub(255, 255, 255);
	glEnableClientState(GL_VERTEX_ARRAY);
	glEnableClientState(GL_COLOR_ARRAY);
	glVertexPointer(2, GL_DOUBLE, sizeof(Particle), &particlesToDisplay[0].x);
	glColorPointer(4, GL_UNSIGNED_BYTE, sizeof(Particle), &particlesToDisplay[0].r);
	glPointSize(3.0);
	glDrawArrays(GL_POINTS, 0, particlesToDisplay.size());
	glDisableClientState(GL_VERTEX_ARRAY);
	glDisableClientState(GL_COLOR_ARRAY);
	glFlush();
	glutSwapBuffers();
}

int main(int argc, char *argv[]){
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_RGBA | GLUT_DEPTH | GLUT_DOUBLE);

	glutInitWindowSize(1920, 1010);
	glutCreateWindow("NBody");

	glutDisplayFunc(display);
	glutTimerFunc(INTERVAL, timer, 1);

	int i, j;
	double distance;
	for (i = 0;i<SIZE;i++) {
		Particle pt;
		pt.x = -100 + (rand() % 100);
		pt.y = -100 + (rand() % 100);
		pt.mass = (rand() % 255) + 50;
		pt.ux = 0;
		pt.uy = 0;
		pt.r = rand() % 255;
		pt.g = rand() % 255;
		pt.b = rand() % 255;
		pt.a = 255;
		particles.push_back(pt);
	}
	//particles[3].mass = 10000;
	vector<double> ForcesX;
	vector<double> ForcesY;
	particlesToDisplay = particles;
	for (int k = 0;k < ITERATION;k++) {
		ForcesX.clear();
		ForcesY.clear();

		double force = 0.0;
		for (i = 0;i < SIZE;i++) {
			force = 0.0;
			#pragma omp parallel for private(distance) reduction(+:force)
			for (j = 0;j < particles.size();j++) {
				if (i == j)
					continue;
				distance = distanceX(particles[i], particles[j]);
				force += (-1) * ((GC*(particles[i].mass*particles[j].mass)) / pow(distance, 2))*((particles[j].x - particles[i].x) / abs(particles[j].x - particles[i].x));
			}
			cout << force << endl;
			ForcesX.push_back(force);
		}
		cout << endl;
		for (i = 0;i < SIZE;i++) {
			force = 0.0;
			#pragma omp parallel for private(distance) reduction(+:force)
			for (j = 0;j < particles.size();j++) {
				if (i == j)
					continue;
				distance = distanceY(particles[i], particles[j]);
				force += (-1) * ((GC*(particles[i].mass*particles[j].mass)) / pow(distance, 2))*((particles[j].y - particles[i].y) / abs(particles[j].y - particles[i].y));
			}
			ForcesY.push_back(force);
		}
		for (i = 0;i < SIZE;i++) {
			double distX = (particles[i].ux*TIMESTAMP) + (0.5*(ForcesX[i] / particles[i].mass))*pow(TIMESTAMP, 2);
			distTravelX.push_back(distX);
			particles[i].x += distX;
			particles[i].ux = particles[i].ux + ((ForcesX[i] / particles[i].mass)*TIMESTAMP);
		}
		for (i = 0;i < SIZE;i++) {
			double distY = (particles[i].uy*TIMESTAMP) + (0.5*(ForcesY[i] / particles[i].mass))*pow(TIMESTAMP, 2);
			distTravelY.push_back(distY);
			particles[i].y += distY;
			particles[i].uy = particles[i].uy + ((ForcesY[i] / particles[i].mass)*TIMESTAMP);
		}
	}
	glutMainLoop();
	return 0;
}