Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

I Binary system simulation in C

  1. Oct 30, 2016 #1
    I'm writing a C program to simulate the orbits of a binary system with a 4th order Runge-Kutta method; basically I'm considering the 2 bodies rotating around a fictitious 3rd body(which would be their center of mass) and I turned down the reciprocal gravitational pulls between the bodies to 0(since the effect of those pulls is what keeps them rotating around their CoM, right?)
    in the beginning, to avoid over-complicating things, I computed the center of mass of the system in the initial state and considered it to stay fixed throughout the simulation: with "good" initial conditions I can get fairly reasonable orbits, but clearly there are cases(the majority of them actually) for which this approximation does not hold and funny things happen, like the outer body flying away while the inner one keeps orbiting about the initial CoM.
    Now I tried to get to the next level and at every iteration I would compute again the position of the CoM in the x-y plane as
    $$x_c = \frac{m_ax_a + m_bx_b}{m_a+m_b}$$
    $$y_c = \frac{m_ay_a + m_by_b}{m_a+m_b}$$
    but if I do it this way when I integrate with RK4 i get some weird results (NaN), so if anyone could give a look at my code I'd appreciate it, but also I wanted to know if there is some flaw in my reasoning that may be leading to this.

    I tried printing out pretty much everything but I can't really identify where the error occurs, even though I think it is most probably in the rk_temp function, since it is where the main calculations are done;
    I noticed that right after the first iteration, the mass of the CoM somehow turns to 0 but still it does not explain why i get NaN; I guess the NaNs could be due to r3 (also in rk_temp) becoming too small and it is approximated to 0 but have no idea how to avoid that;

    Here it is, stripped down to the very essential:

    Code (C):
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>

    #define G 1.981e-29 //UNITS [kg]*[UA]^3*[year]^-2

    //POSITION AND VELOCITY*****************************************************************
    struct state{
      double r;
      double v;

    struct planet{
      struct state x;
      struct state y;
      double m;

    struct planet_sys{
      struct planet P;
      struct planet Q;

    //STRUCT NEEDED FOR EASE OF USE IN rk_temp FUNCTION*******************************
    struct temp{
      double r1, r2, r3, r4;
      double v1, v2, v3, v4;

    struct planet *init(double x0, double y0, double vx0, double vy0, double m);
    struct planet *centerofMass(struct planet_sys S);
    struct planet_sys *rk4(struct planet P, struct planet Q, struct planet C, double dt);
    struct temp *rk_temp(struct state xP, struct state xC, double r, double dt, double Mc);
    double dist_axis(struct state xP, struct state xQ);
    double dist(struct planet A, struct planet B);

    int main(int argc, char *argv[]){
      double x0p, y0p, vx0p, vy0p, x0q, y0q, vx0q, vy0q;
      double mp, mq;
      double t, T, dt;
      int N, i;
      struct planet P, Q, C;
      struct planet_sys PQ;

      //INPUT FOR BODY 1********************************************************************
      printf("P: Inserire m[kg], x0, y0[UA], vx0, vy0[UA/year]: ");
      scanf("%lf %lf %lf %lf %lf", &mp, &x0p, &y0p, &vx0p, &vy0p);
      P = *init(x0p, y0p, vx0p, vy0p, mp);

      //INPUT FOR BODY 2********************************************************************
      printf("P: Inserire m[kg], x0, y0[UA], vx0, vy0[UA/year]: ");
      scanf("%lf %lf %lf %lf %lf", &mq, &x0q, &y0q, &vx0q, &vy0q);
      Q = *init(x0q, y0q, vx0q, vy0q, mq);

      PQ.P = P;
      PQ.Q = Q;
      C = *centerofMass(PQ);

      printf("Inserire T[years], dt[years]: ");
      scanf("%lf %lf", &T, &dt);
      N = (int)(T/dt);

        //APPLY RUNGE-KUTTA AND RE-COMPUTE CENTER OF MASS***************************
      for(i = 1; i < N; i++){
        t = i*dt;
        PQ = *rk4(PQ.P, PQ.Q, C, dt);
        C = *centerofMass(PQ);

    //INITIALIZE PLANET TO THE GIVEN INITIAL DATA***************************************
    struct planet *init(double x0, double y0, double vx0, double vy0, double m){
      static struct planet P;
      P.x.r = x0;
      P.y.r = y0;
      P.x.v = vx0;
      P.y.v = vy0;
      P.m = m;
      return &P;

    //INITIALIZE THE COM AS A 3RD PLANET****************************************************
    struct planet *centerofMass(struct planet_sys S){
      static struct planet C;
      C.m = S.P.m + S.Q.m;
      C.x.r = (S.P.x.r*S.P.m + S.Q.x.r*S.Q.m)/C.m;
      C.y.r = (S.P.y.r*S.P.m + S.Q.y.r*S.Q.m)/C.m;
      C.x.v = (S.P.x.v*S.P.m + S.Q.x.v*S.Q.m)/C.m;
      C.y.v = (S.P.y.v*S.P.m + S.Q.y.v*S.Q.m)/C.m;
      return &C;
    //RUNGE KUTTA METHOD*****************************************************
    struct planet_sys *rk4(struct planet P, struct planet Q, struct planet C, double dt){
      static struct planet_sys PQ_NEW;
      static struct planet Pnew, Qnew;
      static struct temp xP, yP, xQ, yQ;

      xP = *rk_temp(P.x, C.x, dist(P, C), dt, C.m);
      yP = *rk_temp(P.y, C.y, dist(P, C), dt, C.m);
      xQ = *rk_temp(Q.x, C.x, dist(Q, C), dt, C.m);
      yQ = *rk_temp(Q.y, C.y, dist(Q, C), dt, C.m);

      Pnew.x.r = P.x.r + (xP.r1 + 2*xP.r2 + 2*xP.r3 + xP.r4)/6.;
      Pnew.x.v = P.x.v + (xP.v1 + 2*xP.v2 + 2*xP.v3 + xP.v4)/6.;
      Pnew.y.r = P.y.r + (yP.r1 + 2*yP.r2 + 2*yP.r3 + yP.r4)/6.;
      Pnew.y.v = P.y.v + (yP.v1 + 2*yP.v2 + 2*yP.v3 + yP.v4)/6.;
      Qnew.x.r = Q.x.r + (xQ.r1 + 2*xQ.r2 + 2*xQ.r3 + xQ.r4)/6.;
      Qnew.x.v = Q.x.v + (xQ.v1 + 2*xQ.v2 + 2*xQ.v3 + xQ.v4)/6.;
      Qnew.y.r = Q.y.r + (yQ.r1 + 2*yQ.r2 + 2*yQ.r3 + yQ.r4)/6.;
      Qnew.y.v = Q.y.v + (yQ.v1 + 2*yQ.v2 + 2*yQ.v3 + yQ.v4)/6.;
      PQ_NEW.P = Pnew;
      PQ_NEW.Q = Qnew;
      return &PQ_NEW;

    //THE ERROR LIKELY OCCURS HERE*****************************************************************
    struct temp *rk_temp(struct state xP, struct state xC, double r, double dt, double Mc){
      static struct temp X;
      double H = G*Mc;
      double delta_x = dist_axis(xP, xC);
      double r3 = fabs(r*r*r);
      double k = H/r3;
      X.r1 = xP.v*dt;
      X.v1 = (-k*delta_x)*dt;
      X.r2 = (xP.v + 0.5*X.v1)*dt;
      X.v2 = (-k*(delta_x + 0.5*X.r1))*dt;
      X.r3 = (xP.v + 0.5*X.v2)*dt;
      X.v3 = (-k*(delta_x + 0.5*X.r2))*dt;
      X.r4 = (xP.v + X.v3)*dt;
      X.v4 = (-k*(delta_x + X.r3))*dt;
      return &X;

    //COMPUTE DISTANCE ON A SINGLE AXIS********************************************************
    double dist_axis(struct state xP, struct state xQ){
      return (xP.r - xQ.r);

    //COMPUTE DISTANCE IN THE X-Y PLANE******************************************************
    double dist(struct planet A, struct planet B){
      double x = dist_axis(A.x, B.x);
      double y = dist_axis(A.y, B.y);
      return sqrt(x*x + y*y);
  2. jcsd
  3. Oct 30, 2016 #2

    jim mcnamara

    User Avatar

    Staff: Mentor

    Just a comment - NaN can be the result of a computation that exceeds the limits of precision for a floating point datatype. You mention this, I think. Sometimes you can work around the problem by selecting a datatype with greater precision. long double and, depending on your platform, so-called quad precision floating point (128 bytes) might be a choice.

    FORTRAN also has native support for extended precision real (FP) datatypes, if that helps at all.

    macros: isinf(), isnan(), isnormal(), and some others. Check math.h for what your system has. You need to check for these conditions as appropriate.
    floating point is "interesting" in the sense of the ancient Chinese curse -- 'May you live in interesting times'

    Error checking is mandatory in C.
  4. Oct 30, 2016 #3
    I found the error: inside rk4 I created a new struct for each planet and assigned the new values for position and velocities but did not assign anything to the mass, so when I updated the CoM it would get random stuff and made everything go mad.
    Thanks for the help anyway, I'll be more careful in the future.
  5. Oct 30, 2016 #4


    User Avatar
    2017 Award

    Staff: Mentor

    Why don't you reduce the system fully to a 1-body problem? With a fixed CoM, if you know the position of one object you also know the position of the other one - it is sufficient to simulate the motion of one object. If the CoM is moving you can add this constant motion afterwards.
  6. Oct 30, 2016 #5
    Well honestly I thought about it only when I was just terminating, and since I see it more like a programming exercise than anything, I didn't really want to change it all, but I'll surely try and make some efficiency comparison in the next days.
    I still have to do a bit of finetuning, but I also mean to expand it adding more bodies so I'm fine with this approach right now :smile:
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted