Plot the error in Jacobi constant along the orbit with respect to the Jocabi constant of the initial point on the orbit.

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)

ans = 1×2

1 1

size(JC_sL2)

ans = 1×2

1 1

ii

ii =         1905

JC_error = JC_s0sc(ii, :) - JC_sL2(ii, :);

Index in position 1 exceeds array bounds. Index must not exceed 1.

% 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

[Torsten]

For the problem, see above.

ii has value 1905.

JC_s0sc and JC_sL2 are scalars.

Answer 1

[Zuber]

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.

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.