Circular Restricted Three body problem in a No Inertial system

Question

Kevin Valencia 2022-12-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1882162-circular-restricted-three-body-problem-in-a-no-inertial-system

回答： Vinayak 2024-1-16

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

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Vinayak 2024-1-16

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1882162-circular-restricted-three-body-problem-in-a-no-inertial-system#answer_1390316

在 MATLAB Online 中打开

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 