#include <iostream>
#include <cmath>
#include <cstdio>
#include <fstream>
#include <string>
#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;
}
}
/***********************************************************************/
//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(.025,.025,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_t0 = 0;
double atomicTime_t1 = 2;
double atomicTime_t2 = 4;
string filename_t0 = DIRECTORY + "data_t0.txt";
string filename_t1 = DIRECTORY + "data_t1.txt";
string filename_t2 = DIRECTORY + "data_t2.txt";
//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);
/***********************************************************************/
//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;
}