Tuesday, May 18, 2010

Anthony Reporting...

After spending a day and-a-half learning C++ I managed to write my own program to solve the 2-body and N*(2-body) numerically. I actually had some programming experience with Python prior to this summer project so learning C++ was not too difficult.

#include <iostream>

#include <cmath>

#include <cstdlib>

#include <fstream>


using namespace std;


int main ()

{

ofstream dataFile("2body.data");

double r,ax,ay,E,a;

double x = 1; //initial conditions

double y = 0;

double vx = 0;

double vy = 1;


double dt=0.5;

double t = 0;

for (int step =0; t<60 ;step ++) //integrator

{

r = sqrt(x*x + y*y);

ax = -x/(r*r*r);

ay = -y/(r*r*r);

vx += ax * dt ;

vy += ay * dt ;

x += vx * dt ;

y += vy * dt ;

t += dt;

dataFile x << '\t' << y << '\n';//writes data

}


return 0;

}


So this is my code in solving the 2-body problem numerically in C++ (hope you all can read this...). Note that the orbit of the particle is purely circular (e=0) and in Cartesian. In my program you can see that I had to output a data file so that I can plot the x-y data points using gnuplot. Currently, I am only plotting the numerical results and not doing any fancy visualization *yet*. I still need to learn OpenGL! Perhaps Josh can give us all a lecture on OpenGL? ;-)

So using gnuplot, this is my result...



Looks good so far. Moreover, I assigned the particle to move approx. 3 loops and the time-step was 0.1 s (instead of 0.02 s that I was assigned to do). The reason I used 0.1 s because I want you all the see the slight deviation in the radius from my numerical results. The integrator I *think* I used in my program was the leapfrog method (or Euler) so you can see that this integrator is quite stable in this time-step. When I assigned the time-step to be 0.5 s...



This looks interesting because it appears that the particle is undergoing epicyclic motion (is it really undergoing epicyclic motion?). Despite of that you can see that there is a larger deviation in the radius with this time-step. I should also post the energies and angular momentum data as well to see how well they conserve.

The next logical step is to generate N>>1000 particles and have each of them undergo the same motion; hence, N*(2-body) problem. So having just learned the concept of arrays in C++ and the psuedo-random number generator I was able to apply it in this problem.

#include <iostream>

#include <cmath>

#include <cstdlib>

#include <fstream>


using namespace std;


int main ()

{

ofstream dataFile("nbody.data");

int m,n= 4;

double x[n][m],y[n][m],vx[n][m],vy[n][m],ax[n][m],ay[n][m],r[n][m],dt;

dt = 0.05;

srand(time(0));

for (int i=0; i//filling each variable array

{

x[0][i]=0.5+(rand()/(RAND_MAX+1.0))*2.5;

y[0][i]=0;

vx[0][i]=0;

vy[0][i]=sqrt(1/x[0][i]);

double t=0;

for (int j=1; t<40; j++)

{ r[j-1][i]=sqrt(x[j-1][i]*x[j-1][i]+y[j-1][i]*y[j-1][i]);

ax[j-1][i]=-x[j-1][i]/(r[j-1][i]*r[j-1][i]*r[j-1][i]);

ay[j-1][i]=-y[j-1][i]/(r[j-1][i]*r[j-1][i]*r[j-1][i]);

vx[j][i] = vx[j-1][i]+(ax[j-1][i])*(dt);

vy[j][i] = vy[j-1][i]+(ay[j-1][i])*(dt);

x[j][i] = x[j-1][i]+vx[j-1][i]*dt;

y[j][i] = y[j-1][i]+vy[j-1][i]*dt;

dataFile x[j][i] << '\t' << y[j][i] << '\n';//Need to find another way to write the data properly

t += dt;

}

}

dataFile.close();

return 0;


}


Looks messy and long, doesn't it? Remember, I'm still a novice in C++! So initially, I assigned to compute for 4 particles and here's the output...



Still looks good. Actually, there were some problems in the coding prior to this result where particles located farther than 1 unit would fly off its orbit. Here's a sample of what I'm referring to...



It does look kind of funny but interesting at the same time. Notice that each particle seem to have a different center in their orbit, which puzzled me for a while. I sort of consulted with Prof. Artymowicz and he suggested it might be the initial conditions so I went and check and indeed it was! My code for the N*(2-body) program was based on the 2-body program I wrote earlier and it I might have forgotten to change the initial conditions. Silly mistake!

So then I tried to compute for 100 particles to see how it looks but apparently there is some bug in XCode (or in my code). Oh yes, I should mention that I am using XCode as my IDE to compile my programs since I'm mainly using a macbook pro which has two nvidia cards that are suitable for CUDA (perhaps not suitable for dealing really massive computations). I believe this bug has to do with with the memory storage in the arrays so I would need to look at pointers for C++ or something? Maybe one of my colleagues could improve my program to make it more efficient. This is where I start to appreciate more with the idea of CUDA where it would allow us to do some fast, efficient computations for modeling in astrophysics.

So this is where I'm currently at with the programming and now I'm learning OpenGL on my own so that I can do some 'real' animated, visual work with my program. Hopefully it won't not take too long to learn.

Cheers,

Anthony

No comments:

Post a Comment