Plot the error in Jacobi constant along the orbit with respect to the Jocabi constant of the initial point on the orbit.
5 次查看(过去 30 天)
显示 更早的评论
I try to plot the following code:
---------------------------------------------------------------------------------------------------------------------------------
%% Clear the MATLAB Workspace.
clear
clc
close all
format compact
format long
%% Computation of Jacobi Constant for a Periodic Orbit.
% Constants, initial conditions, and integration setup.
s0_sc = [ 1.1745;
0;
0;
0;
-0.111503288786392;
0];
MU_Earth_Moon = 0.0121505856; % Mass parameter of the Earth-Moon system.
n = 1; % Mean-motion of the Moon's orbit around the Earth.
tolerance = 1e-13;
options = odeset('RelTol', tolerance, 'AbsTol', tolerance);
tspan = [ 0 3.393216168933766 ];
% Integrate the inital conditions along the full orbital period.
[t, s_sc] = ode45(@(t, x) CR3BP(t, x, MU_Earth_Moon, n), tspan, ...
s0_sc, options);
% Compute the Jacobi Constant along the orbit.
JC = zeros(size(s_sc, 1), 1);
for ii = 1 : size(s_sc, 1)
JC(ii, 1) = Jacobi_Constant(s_sc(ii, :), MU_Earth_Moon);
end
% Evaluate the Jacobi Constant at the L2 coordinates, s_L2, and initial
% conditions, s0_sc.
s_L2 = [ 1.55682165407869;
0;
0;
0;
0;
0];
JC_s0sc = Jacobi_Constant(s0_sc, MU_Earth_Moon);
JC_sL2 = Jacobi_Constant(s_L2, MU_Earth_Moon);
% Compute the error in Jacobi Constant along the orbit wrt to initial
% point.
size(JC_s0sc)
size(JC_sL2)
ii
JC_error = JC_s0sc(ii, :) - JC_sL2(ii, :);
% Plot the results.
% Non-dimensional distance from Moon to the L2 point.
gamma2 = 0.167832751007869;
figure(1)
scatter3(1 - MU_Earth_Moon, 0, 0, 'k*', 'LineWidth', 2.5, ...
'DisplayName', 'Moon'), hold on,
scatter3(1 - MU_Earth_Moon + gamma2, 0, 0, 'r*', 'LineWidth', 1.5, ...
'DisplayName', 'L_2'),
plot3(s_sc(:,1), s_sc(:,2), s_sc(:,3), 'r', 'LineWidth', 1,...
'DisplayName', 'L_2 Lyapunov Orbit'), grid on,
xlabel('x [NDU]'), ylabel('y [NDU]'), zlabel('z [NDU]'),
set(gca,'fontsize',10), set(gca,'fontweight','bold'),
set(gca,'gridalpha',1), fig1.Color = 'white';
legend('location','northeast')
figure(2)
plot(t, JC_error, 'r', 'LineWidth', 1), grid on,
xlabel('Time [NTU]'), ylabel('Jacobi Constant Error [NDU]'),
set(gca,'fontsize',10), set(gca,'fontweight','bold'),
set(gca,'gridalpha',1), fig2.Color = 'white';
---------------------------------------------------------------------------------------------
and I get the following error:
Index in position 1 exceeds array bounds. Index must not
exceed 1.
Error in Jacobi_Constant_Main (line 60)
JC_error = JC_s0sc(ii, :) - JC_sL2(ii, :);
I have looked at this code forever and I cannot figure what the problem is.
Other functions:
%% ODE45 function for integrating CR3BP EoM.
function [dS] = CR3BP(t, s, MU_Earth_Moon, n)
% Mass parameter of the Earth-Moon system, non-dimensional units.
mu = 0.0121505856;
% Initialize state at input epoch t.
dS = zeros(6,1);
% Retrieve components from current state.
x = s(1,1);
y = s(2,1);
z = s(3,1);
vx = s(4,1);
vy = s(5,1);
vz = s(6,1);
% Compute the distance from the smaller primary.
r_1 = sqrt((x + mu - 1)^2 + y^2 + z^2);
% Compute the distance from the larger primary.
r_2 = sqrt((x + mu)^2 + y^2 + z^2);
% Assign velocities along x, y and z axes of the CR3BP rotating-frame.
dS(1,1) = vx;
dS(2,1) = vy;
dS(3,1) = vz;
% Assign accelerations along x, y and z axes of the CR3BP rotating-frame,
% using equations of motion from the CR3BP model.
dS(4,1) = x + 2*vy - (1 - mu)*(x + mu)/(r_2^3) - ...
mu*(x - 1 + mu)/(r_1^3);
dS(5,1) = y - 2*vx - ((1 - mu)*y)/(r_2^3) - (mu*y)/(r_1^3);
dS(6,1) = -((1 - mu)*z)/(r_2^3) - (mu*z)/(r_1^3);
end
function val = Jacobi_Constant(s, mu)
r_1 = sqrt((s(1) + mu)^2 + s(2)^2 + s(3)^2);
r_2 = sqrt((s(1) - 1 + mu)^2 + s(2)^2 + s(3)^2);
Omega = 1;
U = (1 - mu)/r_1 + mu/r_2 + (1/2)*(Omega^2)*(s(1)^2 + s(2)^2);
Vsq = (s(4)^2 + s(5)^2 + s(6)^2);
val = 2*U - Vsq;
end
1 个评论
回答(1 个)
Zuber
2023-4-10
Hi,
I understand that you are facing an error while plotting the error in Jacobi constant. As per the error message, the index at position ‘1’ of variable ‘JC_s0sc’ and ‘JC_sL2’ in the following line of code is exceeding the array bounds.
JC_error = JC_s0sc(ii, :) - JC_sL2(ii, :);
This means that the array is being accessed at an index which is outside the range of its specified dimension.
Note that the value of variable ‘ii’ is 1095 while the size of ‘JC_s0sc’ and ‘JC_sL2’ is 1x2. Hence, when we try to access both the arrays at 1095
position, MATLAB throws an error since both the arrays have only one row.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1350914/image.png)
Replacing the above line of code with the following will fix the error:
JC_error = JC_s0sc(1, :) - JC_sL2(1, :);
I hope this resolves the error.
0 个评论
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!