Circular Restricted Three body problem in a No Inertial system

Hi!

I am trying to programming de CRTBP, but i am trying to simulate in the No intertial system but my code apparently gives me the solution of the Inertial system, i am trying to find the solution in the No inertial system, does anyone has an advice?

Apparently this is the solution for the Inertial system cause I obtained a spiral, i am using this motion equations for this:

Circular Restricted Three Body Problem

Here is my Matlab code:

clc;
clear all;
% Renombrar variables de las ecuaciones de movimientos
% x=x(1); x'=x(2); y=x(3); y'=x(4)
% x''=x*ome^2+2*ome*y'-G*m1*(x+mu2*r12)/(((x-r1)^2+y^2)^1.5)-G*m2*(x-mu2*r12)/(((x-r2)^2+y^2)^1.5)
% y''=y*ome^2-2*ome*x'-G*m1*y/(((x-r1)^2+y^2)^1.5)-G*m2*y/(((x-r2)^2+y^2)^1.5)
% x=x(1) , x1'=x(2) . y=x(3) , x3'=x(4)
% r es la posición de la masa "m" r(x,y)
% r1 es la posición de la masa "m1" r(x1,0)
% r1-r=(x1-x,0) la diferencia de los vectores
% Condiciones iniciales
x0 = [35,0]; %AU
y0= [0,35]; %AU
%Constantes
G = 1;
m1 =5.974e24; 
m2 = 7.348e22; 
 r1= -21.3; %AU
 r2= 33.7; %AU
 r12=r2-r1;
mu1=((m1)/(m1+m2));
mu2=((m2)/(m1+m2));
%ome=sqrt((m1+m2)/(-r1+r2)^3); %Velocidad angular, la ley de Kepler
ome=1;
[t,x]=ode45(@crtbp,[0 100],[35;0;0;0.03]);

figure
plot3(t,x(:,1),x(:,3));

figure
hold on;
plot(t,x(:,1),'r');
plot(t,x(:,3),'b');

pause(0.01);

figure
hold on
plot(t,x)

pause(0.01);

function dxdt=crtbp(t,x)
G=1;
m1 =5.974e24; 
m2 = 7.348e22;
r1=-21.3;
r2=33.7; 
r12=r2-r1;
mu1=((m1)/(m1+m2));
mu2=((m2)/(m1+m2));
%  ome=sqrt((m1+m2)/(-r1+r2)^3);
 ome=1;
% dxdt=[x(2); x(1)*ome^2+2*ome*x(4)-G*m1*(x(1)+r1)/(((x(1)-r1)^2+x(3)^2)^1.5)-G*m2*(x(1)-r2)/(((x(1)-r2)^2+x(3)^2)^1.5); x(4); x(3)*ome^2-2*ome*x(2)-G*m1*1*x(3)/(((x(1)-r1)^2+x(3)^2)^1.5)-G*m2*x(3)/(((x(1)-r2)^2+x(3)^2)^1.5)];
dxdt=[x(2); x(1)*ome^2+2*ome*x(4)-G*m1*(x(1)+mu2*r12)/(((x(1)-r1)^2+x(3)^2)^1.5)-G*m2*(x(1)-mu1*r12)/(((x(1)-r2)^2+x(3)^2)^1.5); x(4); x(3)*ome^2-2*ome*x(2)-G*m1*1*x(3)/(((x(1)-r1)^2+x(3)^2)^1.5)-G*m2*x(3)/(((x(1)-r2)^2+x(3)^2)^1.5)];
end

Answer 1

[Vinayak]

Hi Kevin,

Upon analysing the code, and the equations you wish to recreate, it is my understanding that for a non-inertial frame you need to consider the centrifugal force and Coriolis force.

Coriolis Term:

For the x-equation: +2 * omega * x(4)

For the y-equation: -2 * omega * x(2)

Centrifugal Term:

For the x-equation: -omega^2 * x(1)

For the y-equation: -omega^2 * x(3)

To achieve the resulting cylindrical spirals for the non-inertial frame, please modify the function to calculate ‘crtbp’ as shown below,

function dxdt = crtbp(t, x) 
    % Define constants for the two primary bodies and the rotating frame 
    G = 1; 
    m1 = 5.974e24; % Mass of the first primary body 
    m2 = 7.348e22; % Mass of the second primary body 
    r1 = -21.3;    % Position of the first primary body on the x-axis 
    r2 = 33.7;     % Position of the second primary body on the x-axis 
    r12 = r2 - r1; % Distance between the two primary bodies 
    mu1 = m1 / (m1 + m2); % Gravitational parameter for the first primary body 
    mu2 = m2 / (m1 + m2); % Gravitational parameter for the second primary body 
    ome = 1; % Angular velocity of the rotating frame 
    
    % Initialize the derivative vector 
    dxdt = zeros(4, 1); 

    % Equations of motion in the rotating frame: 
    % dxdt(1) represents the derivative of the x position, which is the x velocity. 
    dxdt(1) = x(2); 
    % dxdt(2) represents the derivative of the x velocity, including centrifugal force, 
    % Coriolis force, and the gravitational attractions from the two primary bodies. 
    dxdt(2) = x(1) * ome^2 + 2 * ome * x(4) - G * m1 * (x(1) + mu2 * r12) / (((x(1) - r1)^2 + x(3)^2)^1.5) - G * m2 * (x(1) - mu1 * r12) / (((x(1) - r2)^2 + x(3)^2)^1.5); 
    % dxdt(3) represents the derivative of the y position, which is the y velocity. 
    dxdt(3) = x(4); 
    % dxdt(4) represents the derivative of the y velocity, including centrifugal force, 
    % Coriolis force, and the gravitational attractions from the two primary bodies. 
    dxdt(4) = x(3) * ome^2 - 2 * ome * x(2) - G * m1 * x(3) / (((x(1) - r1)^2 + x(3)^2)^1.5) - G * m2 * x(3) / (((x(1) - r2)^2 + x(3)^2)^1.5); 
end 

Hope this helps!