- #1
mastrofoffi
- 51
- 12
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:
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:
#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;
};
//MADE UP OF 2 "STATE" STRUCTS SO I CAN WORK SEPARATELY ON X AND Y COMPONENTS
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);
//INPUT T(TOTAL INTEGRATION TIME) AND dt(INTEGRATION TIME STEP)****************
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);
}
exit(EXIT_SUCCESS);
}//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;
}
//FUNCTION NEEDED TO COMPUTE THE INCREMENT VARIABLES USED IN RK4*************
//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);
}