Guest User

nbodycode

a guest
Apr 19th, 2015
219
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
C++ 2.12 KB | None | 0 0
  1. #include <iostream>
  2. #include <ctime>
  3. #include <math.h>
  4. using namespace std;
  5. #define N 100
  6. #define G 6.673e-11
  7. #define TIMESTAMP 1e11
  8. struct Particle{
  9.     double rx, ry;//position components
  10.     double vx, vy;//velocity components
  11.     double fx, fy;//force components
  12.     double mass;//mass of the particle
  13.  
  14. };
  15. Particle Update(Particle p, double timestamp)
  16. {
  17.     p.vx += timestamp*p.fx / p.mass;
  18.     p.vy += timestamp*p.fy / p.mass;
  19.     p.rx += timestamp*p.vx;
  20.     p.ry += timestamp*p.vy;
  21.     return p;
  22. }
  23. void PrintParticle(Particle p)
  24. {
  25.     printf("rx == %f ry == %f vx == %f vy == %f mass == %f\n", p.rx,p.ry,p.vx,p.vy,p.mass);
  26. }
  27. //Reset the forces on particle
  28. Particle ResetForce(Particle p)
  29. {
  30.     p.fx = 0.0;
  31.     p.fy = 0.0;
  32.     return p;
  33. }
  34. //Add force to particle a by particle b
  35. Particle AddForce(Particle a,Particle b)
  36. {
  37.     double EPS = 3E4;      // softening parameter (just to avoid infinities)
  38.     double dx = b.rx - a.rx;
  39.     double dy = b.ry - a.ry;
  40.     double dist = sqrt(dx*dx + dy*dy);
  41.     double F = (G * a.mass * b.mass) / (dist*dist + EPS*EPS);
  42.     a.fx += F * dx / dist;
  43.     a.fy += F * dy / dist;
  44.     return a;
  45.  
  46. }
  47.  
  48. int main()
  49. {
  50.     Particle particles[N];
  51.     srand(time(NULL));
  52.     //randomly generating N Particles
  53.     for (int i = 0; i < N; i++){
  54.         double rx = 1e18*exp(-1.8)*(.5 - rand());
  55.         particles[i].rx = rx;
  56.         double ry = 1e18*exp(-1.8)*(.5 - rand());
  57.         particles[i].ry = ry;
  58.         double vx = 1e18*exp(-1.8)*(.5 - rand());
  59.         particles[i].vx = vx;
  60.         double vy = 1e18*exp(-1.8)*(.5 - rand());
  61.         particles[i].vy = vy;
  62.         double mass = 1.98892e30*rand()*10 + 1e20;
  63.         particles[i].mass = mass;
  64.  
  65.     }
  66.    
  67.     int numberofiterations = 10;
  68.     int count = 0;
  69.     while (count < numberofiterations){
  70.         for (int i = 0; i < N; i++)
  71.         {
  72.             particles[i] = ResetForce(particles[i]);
  73.             for (int j = 0; j < N; j++)
  74.             {
  75.                 if (i != j)
  76.                 {
  77.                     particles[i] = AddForce(particles[i], particles[j]);
  78.                 }
  79.  
  80.             }
  81.         }
  82.         //loop again to update the time stamp here
  83.         for (int i = 0; i < N; i++)
  84.         {
  85.             particles[i] = Update(particles[i], TIMESTAMP);
  86.         }
  87.         for (int i = 0; i < N; i++)
  88.         {
  89.             PrintParticle(particles[i]);
  90.         }
  91.         count++;
  92.     }
  93.  
  94.    
  95.    
  96.  
  97.     return 0;
  98. }
Advertisement
Add Comment
Please, Sign In to add comment