#include <iostream> 
#include <cmath>
#include <cstdio>
#include <fstream>
#include <string>
#include <sstream>
#include <assert.h>
#include <GLUT/glut.h>
 
using namespace std;
 
//Global constants
const double PI = 4*atan(1.0);
double* xArray_ptr;
double* yArray_ptr;
GLuint circle;
string DIRECTORY = "/Users/rj/Desktop/";
 
//Math section
 
//Generate random trajectories
/***********************************************************************/
double linearCongruence(long int scale, unsigned long int& seed, unsigned long int offSet, unsigned long int modulo)
{
    double randomNum = 0;
    seed = (scale*seed+offSet)%modulo;
    randomNum = (double)seed/(double)(modulo-1);
    return(randomNum);
}
 
void genRandNum(double rand1[], double rand2[], int numRandomNums)
{
    //Constants for random number generator, choices from Numerical Recipes
    const unsigned long int modulo = pow(2.0, 32.0);
    const long int scale = 1664525;
    const unsigned long int offSet = 1013904223;
 
    //Random number generator seed and pointer
    unsigned long int seed1 = 1024;
    unsigned long int seed2 = 512;
 
    for(int i = 0; numRandomNums > i; i++)
    {
        rand1[i] = (linearCongruence(scale,seed1,offSet,modulo));
        rand2[i] = (linearCongruence(scale,seed2,offSet,modulo));
    }
}
 
void boxMuller(double gauss1[], double gauss2[], double rand1[], double rand2[])
{
    for(int i = 0; 10000 > i; i++)
    {
        gauss1[i] = sqrt(-2*log(rand1[i]))*cos(2.0*PI*rand2[i]);
        gauss2[i] = sqrt(-2*log(rand1[i]))*sin(2*PI*rand2[i]);
    }
}
/***********************************************************************/
 
//Calculate classical dynamics
/***********************************************************************/
void posClassical( double posArray[], double gauss1[], double gauss2[], const double FIELD, double currentTime)
{
	//mass is 1
 
	for(int i = 0; 10000 > i; i++)
	{
		posArray[i] = (gauss1[i] + gauss2[i] * currentTime - 0.5 * FIELD * pow(currentTime, 2));
	}
}
 
void momentumClassical(double momArray[], double gauss2[], const double FIELD, double currentTime)
{
 
	//mass is 1
 
	for(int i = 0; 10000 > i; i++)
	{
		momArray[i] = (gauss2[i] - FIELD * currentTime);
	}
}
/***********************************************************************/
 
//Calculate probability densities analitically
/***********************************************************************/
 
//Spatial probability density
double rho_X_T(double pos, double spatialWidth0, double momWidth0, double atomicTime, double potential)
{
    double spatialWidthT = sqrt(pow(spatialWidth0,2.0)+pow(momWidth0,2.0)*pow(atomicTime,2.0));
    double exponent = (-0.5)*pow(((pos+(potential/2.0)*pow(atomicTime,2.0))/spatialWidthT),2.0);
 
    return((1/(sqrt(2.0)*PI))*(1/spatialWidthT)*exp(exponent));
}
/***********************************************************************/
 
//Statistics
/***********************************************************************/
 
 
double arrayAverage(double array[], int arraySize)
{
	double sum = 0;
 
	for(int i = 0; arraySize > i; i++)
	{
		sum = sum + array[i];
	}
	//cout << sum << endl;
    return(sum/(double)arraySize);
}
 
double calculateAvgPosition(double posArray[], double arraySize)
{
	return(arrayAverage(posArray, arraySize));
}
 
double calculateAvgMomentum(double momArray[], double arraySize)
{
	return(arrayAverage(momArray, arraySize));	
}
 
void writePhaseSpacePointToFile(double avgPosition, double avgMomentum, double currentTime)
{
	string filename = DIRECTORY + "phaseSpaceData.txt";
	ofstream fout;
	fout.open(filename.c_str(), ios_base::app);
 
	fout << avgPosition << ";" << avgMomentum << ";" << currentTime << endl;
 
	fout.close();
}
 
 
//If I switch to vectors I can get rid of all this arraySize crap...
void calculatAvgPhaseSpace(double totalTime, double numTimeSteps, 
						   double posArray[], double momArray[], 
						   int arraySize, const double FIELD, double gauss1[], double gauss2[])
{
 
	const double deltaT = totalTime/(double)numTimeSteps;
 
	double currentTime = 0;
	double avgPosition = 0; //These two seem to be a waste of memory
	double avgMomentum = 0;
 
	for(int i = 0; numTimeSteps > i; i++)
	{
		currentTime = deltaT * i;
 
		//load the position and momentum values for each particle
		posClassical(posArray, gauss1, gauss2, FIELD, currentTime);
		momentumClassical(momArray, gauss2, FIELD, currentTime);
 
		//calculate the averages
		avgPosition = calculateAvgPosition(posArray, arraySize);
		avgMomentum = calculateAvgMomentum(momArray, arraySize);
 
		//Write data to file
		writePhaseSpacePointToFile(avgPosition, avgMomentum, currentTime);
	}
}
 
 
/***********************************************************************/
 
//Bookkeeping
/***********************************************************************/
void writeArrayToFile(double xArray[], double yArray[], ofstream* filePtr)
{
	//Arrays must be same size, should add an assert here at some point.
	for(int i = 0; 10000 > i; i++)
	{
		*filePtr << xArray[i] << ";" << yArray[i] << endl;
	}
}
 
void buildSequentialFileNameArray(string nameArray[], string directory, string prefix, string extenstion, const int numFiles)
{
	stringstream ss;
	string suffix;
	for(int i = 0; numFiles > i; i++)
	{
		ss << i;
		suffix = "_" + ss.str() + "." + extenstion;
		nameArray[i] = directory + prefix + suffix;
	}
}
 
void buildFileStreamArray(ofstream* streamArray[], string nameArray[], int numStreams)
{
	for(int i = 0; numStreams > i; i++)
	{
		streamArray[i] = new ofstream(nameArray[i].c_str()) ;
	}
}
 
/***********************************************************************/
 
//Graphics part
/***********************************************************************/
void createCircle()
{
    float theta;
    int i;
    int div = 50;
 
 
    circle = glGenLists(1);
 
    glNewList(circle,GL_COMPILE);
    glBegin(GL_TRIANGLE_FAN);
    glVertex2f(0,0);
    for (i=0;i<div;i++)
    {
        theta = (float)i*2*3.14159/div;
        glVertex2f(cos(theta),sin(theta));
    }
    glVertex2f(1,0);
    glEnd();
    glEndList();
}
 
void drawAxis(int xSize, int ySize)
{
	glBegin(GL_LINES);
	//x-axis
	glColor3f(0,0,0);
	glVertex3f(-xSize,0,0);
	glVertex3f(xSize,0,0);
 
	//y-axis
	glColor3f(0,0,0);
	glVertex3f(0,-ySize,0);
	glVertex3f(0,ySize,0);
	glEnd();
}
 
void drawBaorder(int xSize, int ySize)
{
	glBegin(GL_LINES);
 
	//Left boarder
	glColor3f(0,0,0);
	glVertex3f(-xSize,-ySize,0);
	glVertex3f(-xSize,ySize,0);
 
	//Botom boarder
	glColor3f(0,0,0);
	glVertex3f(-xSize,-ySize,0);
	glVertex3f(xSize,-ySize,0);
 
	//Right boarder
	glColor3f(0,0,0);
	glVertex3f(xSize,-ySize,0);
	glVertex3f(xSize,ySize,0);
 
	//top boarder
	glColor3f(0,0,0);
	glVertex3f(-xSize,ySize,0);
	glVertex3f(xSize,ySize,0);
	glEnd();
}
 
void display(void)
{
    glClear(GL_COLOR_BUFFER_BIT);
    glMatrixMode(GL_MODELVIEW);
    glLoadIdentity();
 
    glColor3f(1,0,0); //pure red
    glPointSize(2.0);
/*
    glBegin(GL_POINTS);
    for (int i=0; i < 10000; i++)
    {
        glVertex2f(xArray_ptr[i],yArray_ptr[i]);
    }
    glEnd();
*/    
	drawAxis(35,35);
	drawBaorder(35,35);
 
 
    glColor3f(1,0,0);
    for (int i=0; i < 10000; i++)
    {
        glPushMatrix();
        glTranslatef(xArray_ptr[i],yArray_ptr[i],0);
        glScalef(.05,.05,1);
        glCallList(circle);
        glPopMatrix();
    }
 
    glutSwapBuffers();
}
 
void reshape(int width, int height)
{
    glMatrixMode(GL_PROJECTION);
	glClearColor(1,1,1,1);
    glLoadIdentity();
    gluOrtho2D(-40,40,-40,40); //you need to get these
    glViewport(0,0,width,height);
 
    createCircle();
}
 
void idle(void)
{
    //Need to generate x and y
 
    glutPostRedisplay();
}
 
/***********************************************************************/
 
//main program
/***********************************************************************/
int main(int argc, char *argv[])
{
    //OpenGL stuff
    glutInit(&argc, argv);
 
    glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_DEPTH);
    glutInitWindowSize(1024, 768);
 
    glutCreateWindow("Classical Statistical Dynamics");
 
    glutDisplayFunc(display);
    glutReshapeFunc(reshape);
    //glutIdleFunc(idle);
 
	//Problem specific constants
    const int numRandomNums = 10000;
	const double EFIELD = 1.0;
    const double spatialWidth0 = 1.0;
    const double momWidth0 = 0.5;
	const double numTimeSteps = 100;
	double atomicTime[] = {0,2,4};
	double atomicTime_t0 = 0;
	double atomicTime_t1 = 2;
	double atomicTime_t2 = 4;
	const int numFiles = 3;
	string fileNameArray[numFiles];
	ofstream* streamArray[numFiles];
	string filename_t0 = DIRECTORY + "data_t0.txt";
	string filename_t1 = DIRECTORY + "data_t1.txt";
	string filename_t2 = DIRECTORY + "data_t2.txt";
 
	buildSequentialFileNameArray(fileNameArray, DIRECTORY, "data_t", "txt", numFiles);
	buildFileStreamArray(streamArray, fileNameArray, numFiles);
 
	//Stream variables for classical trajectories
	ofstream fout0(filename_t0.c_str());
	ofstream fout1(filename_t1.c_str());
	ofstream fout2(filename_t2.c_str());
 
    //temp arrays for storing uniformly distributed random numbers
    double rand1[numRandomNums];
    double rand2[numRandomNums];
 
    //temp arrays for storing gaussian distributed random numbers
    double gauss1[numRandomNums]; 
    double gauss2[numRandomNums];
 
	//Arrays to hold classical position and momentum
	double posArray[numRandomNums];
	double momArray[numRandomNums];
 
    //Array pointers
    xArray_ptr = &posArray[0];
    yArray_ptr = &momArray[0];
 
    //Generate 10000 uniformly distributed random numbers
    genRandNum(rand1, rand2, numRandomNums);
 
    //Convert to Gaussian distributed random numbers
    boxMuller(gauss1, gauss2, rand1, rand2);
 
	//PROBLEM 1.2F
    /***********************************************************************/	
	calculatAvgPhaseSpace(atomicTime_t2, numTimeSteps, posArray, momArray, numRandomNums, EFIELD, gauss1, gauss2);
 
    /***********************************************************************/
 
	for(int i = 0; numFiles > i; i++)
	{
		posClassical(posArray, gauss1, gauss2, EFIELD, atomicTime[i]);
		momentumClassical(momArray, gauss2, EFIELD, atomicTime[i]);
 
		writeArrayToFile(posArray, momArray, streamArray[i]);
	}
/*	
    //PROBLEM 1.2D T0
    /***********************************************************************/
 
	//Calculate first itteration of classical position and momentum for t = 0
	posClassical(posArray, gauss1, gauss2, EFIELD, atomicTime_t0);
	momentumClassical(momArray, gauss2, EFIELD, atomicTime_t0);
 
	writeArrayToFile(posArray, momArray, &fout0);
    /***********************************************************************/
 
    //Problem 1.2D T1
    /***********************************************************************/
	//Calculate first itteration of classical position and momentum for t = 2
	posClassical(posArray, gauss1, gauss2, EFIELD, atomicTime_t1);
	momentumClassical(momArray, gauss2, EFIELD, atomicTime_t1);
 
	writeArrayToFile(posArray, momArray, &fout1);
    /***********************************************************************/
 
	//Problem 1.2D T2
    /***********************************************************************/
	//Calculate first itteration of classical position and momentum for t = 4
	posClassical(posArray, gauss1, gauss2, EFIELD, atomicTime_t2);
	momentumClassical(momArray, gauss2, EFIELD, atomicTime_t2);
 
	writeArrayToFile(posArray, momArray, &fout2);
    /***********************************************************************/
 
    //Need to make sure arrays are assigned and changing correctly before this function is called.
    glutMainLoop();
 
    return 0;
}