The Code-Bin
Links
Home
Add your code!
All Listings
About
Latest Entry
Featured Scripts
Author's Website
Latest Entries
FFMPEG Thumbnail Scr...
PHP, 0.8KB
Jul. 29, 10:24pm
John
Z80 Assembler, 190 bytes
Feb. 17, 3:36am
John
Z80 Assembler, 176 bytes
Sep. 13, 2:19am
John
Z80 Assembler, 77 bytes
Sep. 13, 2:18am
John
Z80 Assembler, 209 bytes
Sep. 13, 2:17am
untitled C++ Code
Posted by: RJ | September 11, 2009 @ 8:06pm
C++ Code
[
Download
]
#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; }
Syntax Highlighting
[
Open in new window
]
Author Comments
Generates data and plot.
Rating
4.59 / 8
63 Votes
http://codebin.yi.org/373
page generated in 0.00 seconds