title: generalizing the force in n-body simulation
body: <p>I have a problem generalizing the gravitational force when programming an n-body problem in c++. If i try to rewrite my program to be easily adaptable for an arbitrary number of bodies i ran into a problem when trying to generalize the force.</p>
  
  <p>I wrote the following code:</p>
  
  <pre><code>#include &lt;cstdlib&gt;
  #include &lt;iostream&gt;
  #include &lt;cmath&gt;
  #include &lt;fstream&gt;
  #define h  1000.0
  #define G  6.67384*pow(10.0,-11)
  
  using namespace std;
  
  
  class particle{
        public:
        double kx1,kx2,kx3,kx4, kv1, kv2, kv3, kv4;
        double ky1, ky2, ky3, ky4, kvy1, kvy2, kvy3, kvy4;
        double x,y,vx,vy,m;
  
  
        double dist(particle aap){
               double dx = x - aap.x;
               double dy = y - aap.y;
               return sqrt(pow(dx,2.0)+pow(dy,2.0));
               }
  
        double g(double x1, double y1,particle aap){
               return G*aap.m*(aap.x-x1)/pow(dist(aap),3.0);
               }
  
        double p(double x1, double y1, particle aap){
               return G*aap.m*(aap.y-y1)/pow(dist(aap),3.0);
         }
  
  
        void update(){            //zet het object 1 stap vooruit
             x = x + (1/6.0)*(kx1+2*kx2+2*kx3+kx4);
             vx = vx + (1/6.0)*(kv1+2*kv2+2*kv3+kv4);
             y = y + (1/6.0)*(ky1+2*ky2+2*ky3+ky4);
             vy = vy + (1/6.0)*(kvy1+2*kvy2+2*kvy3+kvy4);
             }
  
       void create(double x1, double y1, double vx1, double vy1, double m1){
                        x = x1;
                        y = y1;
                        vx = vx1;
                        vy = vy1;
                        m =m1;
                        }
  
       bool operator ==(particle &amp;other){
            if(x == other.x &amp;&amp; y == other.y &amp;&amp; vx == other.vx &amp;&amp; vy == other.vy){
                 return true;
                 }
                 }
  
  
        };
  
  
  
  particle zon, maan, aarde;
  
  void set(){
       zon.create(1, 1, -2, 1, 2*pow(10.0,30));
       aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
       maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx +  -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
       }
  
  
  
  double xforce(double x1, double y1, particle aap){     //kracht in de x-richting 
  
         particle bodies[] = {zon, aarde, maan};
  
         double fx;
         for (int i = 0; i &lt; 3; i++){
             if (bodies[i].x == aap.x &amp;&amp; bodies[i].y == aap.y &amp;&amp; bodies[i].vx == aap.vx &amp;&amp; bodies[i].vy == aap.vy ){;}
  
             else{
                  fx += aap.g(x1,y1,bodies[i]);
                  }
                  }
         return fx;
         }
  
  double yforce(double x1, double y1, particle aap){ //kracht in de y-richting 
  
         particle bodies[] = {zon, aarde, maan};
  
         double fy;
         for (int i = 0; i &lt;= 3; i++){
             if (bodies[i].x == aap.x &amp;&amp; bodies[i].y == aap.y &amp;&amp; bodies[i].vx == aap.vx &amp;&amp; bodies[i].vy == aap.vy) {;}
  
             else{
                  fy += aap.p(x1,y1,bodies[i]);
                  }
                  }
         return fy;
         }
  
  
  void corr(particle&amp; body){
       body.kx1 = h*body.vx;
       body.kv1 = h*xforce(body.x, body.y, body);
       body.ky1 = h*body.vy;
       body.kvy1 = h*yforce(body.x, body.y, body);
  
       body.kx2 = h*(body.vx + 0.5*body.kv1);
       body.kv2 = h*xforce(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, body);
       body.ky2 = h*(body.vy + 0.5*body.kvy1);
       body.kvy2 = h*yforce(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, body);
  
       body.kx3 = h*(body.vx+ 0.5*body.kv2);
       body.kv3 = h*xforce(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2, body);
       body.ky3 = h*(body.vy+ 0.5*body.kvy2);
       body.kvy3 = h*yforce(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2,body);
  
       body.kx4 = h*(body.vx+body.kv3);
       body.kv4 = h*xforce(body.x+ body.kx3, body.y + body.ky3, body);
       body.ky4 = h*(body.vy + body.kvy3);
       body.kvy4 = h*yforce(body.x + body.kx3, body.y + body.ky3, body);
       }
  
  
  
  
  
  
  
  void bereken(){
       set();
  
       zon.create(1, 1, -2, 1, 2*pow(10.0,30));
       aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
       maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx +  -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
       ofstream file;
       file.open("3body.txt");
       for(int i =0; i &lt;=30000; i++){
               corr(maan);
               corr(zon);
               corr(aarde);
               zon.update();
               aarde.update();
               maan.update();
               file &lt;&lt; i*h &lt;&lt;"  "&lt;&lt; zon.x &lt;&lt; "  "&lt;&lt; zon.y &lt;&lt; "    "&lt;&lt; zon.vx&lt;&lt; "  "&lt;&lt; zon.vy &lt;&lt;"  "&lt;&lt; aarde.x &lt;&lt; "    " &lt;&lt; aarde.y &lt;&lt;"    "&lt;&lt; aarde.vx &lt;&lt;"    " &lt;&lt; aarde.vy &lt;&lt;"   "&lt;&lt; maan.x&lt;&lt;"  "&lt;&lt;maan.y&lt;&lt;"\n";
               }
       file.close();
       }
  
  
  int main()
  {   
      bereken();
      system("pause");
      return 0;
  }
  </code></pre>
  
  <p>I believe the problem lies in the functions xforce() and yforce(). I don't know what exactly goes wrong here. I'll include a working version of my code for 3 bodies here so it will be easier to compare the two.</p>
  
  <pre><code>#include &lt;cstdlib&gt;
  #include &lt;iostream&gt;
  #include &lt;cmath&gt;
  #include &lt;fstream&gt;
  #define h  1000.0
  #define G  6.67384*pow(10.0,-11)
  
  using namespace std;
  
  
  class particle{
        public:
        double kx1,kx2,kx3,kx4, kv1, kv2, kv3, kv4;
        double ky1, ky2, ky3, ky4, kvy1, kvy2, kvy3, kvy4;
        double x,y,vx,vy,m;
  
  
        double dist(particle aap){
               double dx = x - aap.x;
               double dy = y - aap.y;
               return sqrt(pow(dx,2.0)+pow(dy,2.0));
               }
  
        double g(double x1, double y1,particle aap, particle bever){
               return G*aap.m*(aap.x-x1)/pow(dist(aap),3.0) + G*bever.m*(bever.x-x1)/pow(dist(bever),3.0);
               }
  
        double p(double x1, double y1, particle aap, particle bever){
               return G*aap.m*(aap.y-y1)/pow(dist(aap),3.0) + G*bever.m*(bever.y-y1)/pow(dist(bever),3.0);
         }
  
  
        void update(){            //zet het object 1 stap vooruit
             x = x + (1/6.0)*(kx1+2*kx2+2*kx3+kx4);
             vx = vx + (1/6.0)*(kv1+2*kv2+2*kv3+kv4);
             y = y + (1/6.0)*(ky1+2*ky2+2*ky3+ky4);
             vy = vy + (1/6.0)*(kvy1+2*kvy2+2*kvy3+kvy4);
             }
  
       void create(double x1, double y1, double vx1, double vy1, double m1){
                        x = x1;
                        y = y1;
                        vx = vx1;
                        vy = vy1;
                        m =m1;
                        }
  
       bool operator ==(particle &amp;other){
            if(x == other.x &amp;&amp; y == other.y &amp;&amp; vx == other.vx &amp;&amp; vy == other.vy){
                 return true;
                 }
                 }
  
  
        };
  
  
  
  particle zon, maan, aarde;
  
  void set(){
       zon.create(1, 1, -2, 1, 2*pow(10.0,30));
       aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
       maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx +  -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
       }
  
  
  
  double xforce(double x1, double y1, particle&amp; aap){     //kracht in de x-richting 
  
         particle bodies[] = {zon, aarde, maan};
  
         double fx;
         for (int i = 0; i &lt; 3; i++){
             if (bodies[i].x == aap.x &amp;&amp; bodies[i].y == aap.y &amp;&amp; bodies[i].vx == aap.vx &amp;&amp; bodies[i].vy == aap.vy ){;}
  
             else{
                  fx += aap.g(x1,y1,bodies[i],aap);
                  }
                  }
         return fx;
         }
  
  double yforce(double x1, double y1, particle&amp; aap){ //kracht in de y-richting 
  
         particle bodies[] = {zon, aarde, maan};
  
         double fy;
         for (int i = 0; i &lt;= 3; i++){
             if (bodies[i].x == aap.x &amp;&amp; bodies[i].y == aap.y &amp;&amp; bodies[i].vx == aap.vx &amp;&amp; bodies[i].vy == aap.vy) {;}
  
             else{
                  fy += aap.p(x1,y1,bodies[i],aap);
                  }
                  }
         return fy;
         }
  
  
  void corr(particle&amp; body, particle aap, particle bever){
       body.kx1 = h*body.vx;
       body.kv1 = h*body.g(body.x, body.y, bever ,aap);
       body.ky1 = h*body.vy;
       body.kvy1 = h*body.p(body.x, body.y, bever, aap);
  
       body.kx2 = h*(body.vx + 0.5*body.kv1);
       body.kv2 = h*body.g(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, bever,aap);
       body.ky2 = h*(body.vy + 0.5*body.kvy1);
       body.kvy2 = h*body.p(body.x + 0.5*body.kx1, body.y + 0.5*body.ky1, bever,aap);
  
       body.kx3 = h*(body.vx+ 0.5*body.kv2);
       body.kv3 = h*body.g(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2, bever,aap);
       body.ky3 = h*(body.vy+ 0.5*body.kvy2);
       body.kvy3 = h*body.p(body.x + 0.5*body.kx2, body.y + 0.5*body.ky2,bever,aap);
  
       body.kx4 = h*(body.vx+body.kv3);
       body.kv4 = h*body.g(body.x+ body.kx3, body.y + body.ky3, bever,aap);
       body.ky4 = h*(body.vy + body.kvy3);
       body.kvy4 = h*body.p(body.x + body.kx3, body.y + body.ky3, bever,aap);
       }
  
  
  
  
  
  
  
  void bereken(){
       set();
  
       zon.create(1, 1, -2, 1, 2*pow(10.0,30));
       aarde.create(1.5*pow(10.0,11), 0, 2, 29780, 6*pow(10.0,24));
       maan.create(aarde.x + 1, aarde .y + 3.844399*pow(10.0,8), aarde.vx +  -1022.0, aarde.vy + 1, 7.3347*pow(10.0,22));
       ofstream file;
       file.open("3body.txt");
       for(int i =0; i &lt;=30000; i++){
               corr(maan, aarde, zon);
               corr(zon, maan , aarde);
               corr(aarde, maan , zon);
               zon.update();
               aarde.update();
               maan.update();
               file &lt;&lt; i*h &lt;&lt;"  "&lt;&lt; zon.x &lt;&lt; "  "&lt;&lt; zon.y &lt;&lt; "    "&lt;&lt; zon.vx&lt;&lt; "  "&lt;&lt; zon.vy &lt;&lt;"  "&lt;&lt; aarde.x &lt;&lt; "    " &lt;&lt; aarde.y &lt;&lt;"    "&lt;&lt; aarde.vx &lt;&lt;"    " &lt;&lt; aarde.vy &lt;&lt;"   "&lt;&lt; maan.x&lt;&lt;"  "&lt;&lt;maan.y&lt;&lt;"\n";
               }
       file.close();
       }
  
  
  int main()
  {   
      bereken();
      system("pause");
      return 0;
  
  }
  </code></pre>
  
  <p>In this code the total force is calculated in the functions g() and p() instead of in xforce and yforce().</p>