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!
