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

#include "stdafx.h"
#include <stdlib.h>
#include <../glut-3.7.6-bin/glut.h>
#include <math.h>
#include <conio.h>
#include <omp.h>

#define G 6.67428E-11
#define TIMESTEP 250
#define INTERVAL 20
#define DIMENSION 400
#define	ITERATION 1000
#define SCALING_FACTOR 10000

typedef struct ParticleProperty {

	unsigned char r;
	unsigned char g;
	unsigned char b;
	unsigned char a;
}Property;

int SIZE;
int loop = 0;

double * x;
double * y;

Property * property;

/*
void timer(int value) {

	int i = loop * SIZE;
	int bound = (loop +1) * SIZE;
	for (i ; i < bound; i++) {

		//particlesToDisplay[i*2] = particle[i].x * SCALING_FACTOR;
		//particlesToDisplay[i*2 +1] = particle[i].y * SCALING_FACTOR;
	}

	loop++;
	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]);
	glColorPointer(4, GL_UNSIGNED_BYTE, sizeof(Particle), &particlesToDisplay[0]);
	glPointSize(3.0);
	glDrawArrays(GL_POINTS, 0, 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(800, 600);
	glutCreateWindow("N-Body Simulation");

	glutDisplayFunc(display);
	glutTimerFunc(INTERVAL, timer, 1);
	*/
	int i;
	int j;
	int k;

	double ax;
	double ay;
	double rx;
	double ry;
	double r_square;

	double *mass;
	double *ux;
	double *uy;

	FILE *infile;
	FILE *outfile;

	fopen_s(&infile, "input.txt", "r");

	fscanf_s(infile, "%d", &SIZE);

	property = (Property*)malloc(SIZE * sizeof(Property));
	mass = (double*)malloc(SIZE * sizeof(double));
	ux = (double*)malloc(SIZE * sizeof(double));
	uy = (double*)malloc(SIZE * sizeof(double));

	x = (double*)malloc(SIZE * ITERATION * sizeof(double));
	y = (double*)malloc(SIZE * ITERATION * sizeof(double));

	// Set up the particle properties
	for (i = 0; i<SIZE; i++) {

		fscanf_s(infile, "%lf %lf %lf", &x[i], &y[i], &mass[i]);

		property[i].r = rand() % 255;
		property[i].g = rand() % 255;
		property[i].b = rand() % 255;
		property[i].a = rand() % 100 + 125;

		ux[i] = 0;
		uy[i] = 0;
	}

	// Simulation Data

#pragma omp parallel for private()

	for (i = 1; i < ITERATION; i++) {

		for (j = SIZE; j--;) {

			ax = 0;
			ay = 0;

			for (k = SIZE; k--;) {

				if (j == k) continue;

				rx = x[SIZE*i + k] - x[SIZE*i + j];
				ry = y[SIZE*i + k] - y[SIZE*i + j];
				r_square = rx*rx + ry*ry;

				ax += G * mass[j] * rx / r_square;
				ay += G * mass[j] * ry / r_square;
			}

			x[SIZE*i + j] += (ux[j] + ax * TIMESTEP / 2) * TIMESTEP;
			y[SIZE*i + j] += (uy[j] + ay * TIMESTEP / 2) * TIMESTEP;

			ux[j] += ax * TIMESTEP;
			uy[j] += ay * TIMESTEP;
		}
	}

	//glutMainLoop();

	printf_s("finished");
	_getch();

	free(x);
	free(y);
	free(ux);
	free(uy);
	free(mass);

	return 0;
}