"Index Exceeds Matrix Dimensions" Error Within ODE45 Events Function

I am writing a code to produce a simulation of a Hohmann transfer from a parking orbit to a geosynchronous orbit. I am using ODE45 to integrate the equations of motion, and I am trying to use the events feature to stop the integration when the simulation reaches a location corresponding to Y = 0. Here is my call to ODE45.

tspan2 = [0,tfinalphase2];
options2 = odeset('Events', @events2, 'RelTol', 1e-08, 'MaxStep', 100.0);
[t_out2, state_out2] = ode45(@two_body2, tspan2, y2, options2);

Here is the Events function

function [value, isterminal, direction] = events2(~, state_out2)
value = state_out2(end,2);
isterminal = 1;
direction = 0;
end

state_out2 is a state vector consisting of the x,y,z components of location and the x,y,z components of velocity. The second column is the y component of the location - I want the integration to stop when this value becomes 0 (spacecraft passes through X-Z plane).

I am getting the error "Index Exceeds Matrix Dimensions" from the line

value = state_out2(end,2);

What am I doing wrong?

Here is my full code

clear all
close all
clc
Astronautics_Project_Calculations
format long;
global r1
global v1
global T1
global r3
%global va2
global dv1
global vb2
%global vb3
%global dv2
%global a
%global en
%global e
%global T2
%global TT
%global T3
%global MtoKM 
MtoKM = 1.60934;                                        %conversion miles to kilometers
global mu                                            %(km/s)
rx0 = 0;                                                %initial x position
ry0 = r1*cos(0.488692);                                 %initial y position
rz0 = r1*sin(0.488692);                                 %initial z position
rdotx0 = -v1;                             %initial x velocity
rdoty0 = 0;                                             %initial y velocity
rdotz0 = 0;                             %initial z velocity
y = [rx0 ry0 rz0 rdotx0 rdoty0 rdotz0];
%FIRST ORBIT
options1 = odeset('RelTol', 1e-08,'MaxStep',100.0);
tfinalphase1 = 1.25*T1;
tspan1 = [0, tfinalphase1];                             %set time span for phase 1 integration 
[t_out1, state_out1] = ode45(@two_body1, tspan1, y, options1);
%SECOND ORBIT
rx1 = state_out1(end,1) ;
ry1 = state_out1(end,2);
rz1 = state_out1(end,3);
rdotx1 = state_out1(end,4);
rdoty1 = state_out1(end,5);
rdotz1 = state_out1(end,6);
y2 = [rx1 ry1 rz1 rdotx1 rdoty1-(dv1*cos(0.488692)) rdotz1-(dv1*sin(0.488692))];
tfinalphase2 = TT;
tspan2 = [0,tfinalphase2];
options2 = odeset('Events', @events2, 'RelTol', 1e-08, 'MaxStep', 100.0);
[t_out2, state_out2, te2, ye2, ie2] = ode45(@two_body2, tspan2, y2, options2);
%THIRD ORBIT
rx2 = state_out2(end,1);
ry2 = state_out2(end,2);
rz2 = state_out2(end,3);
rdotx2 = state_out2(end,4); 
rdoty2 = state_out2(end,5);
rdotz2 = state_out2(end,6);
y3 = [rx2 ry2 rz2 rdotx2 rdoty2+(dv2*cos(0.488692)) rdotz2+(dv2*sin(0.488692))];
tfinalphase3 = T3;
tspan3 = [0, tfinalphase3+1000];
options3 = odeset('Events', @events3, 'RelTol', 1e-08, 'MaxStep', 100.0);
[t_out3, state_out3, te3, ye3, ie3] = ode45(@two_body3, tspan3, y3, options3);
[k,j,h] = sphere(20);
figure
surf(6378*k,6378*j,6378*h,'edgecolor','k','facecolor','none','linestyle',':');
hold on;
axis equal;
title('Complete Orbit Sequence - Earth Relative Position')
xlabel('X Axis (km)')
ylabel('Y Axis (km)')
zlabel('Z Axis (km)')
%axis([-40000 40000 -40000 40000 -20000 20000])
scatter3(state_out1(:,1),state_out1(:,2),state_out1(:,3),'.','m');
scatter3(state_out2(:,1),state_out2(:,2),state_out2(:,3),'.','c');
scatter3(state_out3(:,1),state_out3(:,2),state_out3(:,3),'.','g');
scatter3(state_out1(1,1),state_out1(1,2),state_out1(1,3), 300,'red','x','Linewidth',3);
scatter3(state_out2(1,1),state_out2(1,2),state_out2(1,3), 300,[1 0.6 0],'x','LineWidth',3);
scatter3(state_out2(end,1),state_out2(end,2),state_out2(end,3), 300,'k','x','LineWidth',2);
legend('Earth','Initial Orbit - 400 Miles', 'Hohmann Transfer Orbit','Geosynchronous Orbit - 22236 Miles','Injection to Initial Orbit','Periapsis of Hohmann Ellipse','Apoapsis of Hohmann Ellipse')
figure
subplot(3,1,1)
hold on
title('Spacecraft Inertial Velocity')
xlabel('Time (s)')
ylabel('Spacecraft X Axis Velocity (km/s)')
scatter(t_out1(:,1),state_out1(:,4),'.','m','LineWidth',1);
scatter((t_out2(:,1)+t_out1(end,1)),state_out2(:,4),'.','c','LineWidth',1);
scatter((t_out3(:,1)+t_out2(end,1)+t_out1(end,1)),state_out3(:,4),'.','g','LineWidth',1);
scatter(t_out1(end,1),state_out1(end,4), 100,'d','k','LineWidth', 1);
scatter((t_out2(end,1)+t_out1(end,1)),state_out2(end,4),100,'x','k','LineWidth', 1);
legend('1st Orbit','2nd Orbit','3rd Orbit','Delta V 1', 'Delta V 2','location','bestoutside')
subplot(3,1,2)
hold on
xlabel('Time (s)')
ylabel('Spacecraft Y Axis Velocity (km/s)')
scatter(t_out1(:,1),state_out1(:,5),'.','m','LineWidth',1);
scatter((t_out2(:,1)+t_out1(end,1)),state_out2(:,5),'.','c','LineWidth',1);
scatter((t_out3(:,1)+t_out2(end,1)+t_out1(end,1)),state_out3(:,5),'.','g','LineWidth',1);
scatter(t_out1(end,1),state_out1(end,5),100,'d','k','LineWidth', 1);
scatter((t_out2(end,1)+t_out1(end,1)),state_out2(end,5),100,'x','k','LineWidth', 1);
legend('1st Orbit','2nd Orbit','3rd Orbit','Delta V 1', 'Delta V 2','location','bestoutside')
subplot(3,1,3)
hold on
xlabel('Time (s)')
ylabel('Spacecraft Z Axis Velocity (km/s)')
axis([0 120000 -10 10]);
scatter(t_out1(:,1),state_out1(:,6),'.','m','LineWidth',1);
scatter((t_out2(:,1)+t_out1(end,1)),state_out2(:,6),'.','c','LineWidth',1);
scatter((t_out3(:,1)+t_out2(end,1)+t_out1(end,1)),state_out3(:,6),'.','g','LineWidth',1);
scatter(t_out1(end,1),state_out1(end,6), 100,'d','k','LineWidth', 1);
scatter((t_out2(end,1)+t_out1(end,1)),state_out2(end,6),100,'x','k','LineWidth', 1);
legend('1st Orbit','2nd Orbit','3rd Orbit','Delta V 1','Delta V 2','location','bestoutside')

disp('Time of Periapsis Achievement (s)')
disp(t_out1(end,1))
disp('Simulation Spacecraft Velocity (km/s)')
disp(((state_out2(1,4)^2)+(state_out2(1,5)^2)+(state_out2(1,6)^2))^0.5)
disp('Planned Spacecraft Velocity (km/s)')
disp(v1+dv1)
disp('Simulation Time of Apoapsis Achievement (s)')
disp(t_out1(end,1)+t_out2(end,1))
disp('Planned Time of Apoapsis Achievement (s)')
disp(1.25*T1+TT)
disp('True Radius of Apoapsis (km)')
disp(((state_out2(end,1)^2)+(state_out2(end,2)^2)+(state_out2(end,3)^2))^0.5)
disp('Planned Radius of Apoapsis (km)')
disp(r3)
disp('True Velocity at Apoapsis (km/s)')
disp(((state_out2(end,4)^2)+(state_out2(end,5)^2)+(state_out2(end,6)^2))^0.5)
disp('Planned Velocity at Apoapsis (km/s)')
disp(vb2)

function ydot1 = two_body1(~,y)
mu = 398600;
r = norm(y(1:3));
ydot1 = [y(4); y(5); y(6); -(mu*y(1))/(r^3); -(mu*y(2))/(r^3); -(mu*y(3))/(r^3)];
return
end

function ydot2 = two_body2(~,y2)
mu = 398600;
r = norm(y2(1:3));
ydot2 = [y2(4); y2(5); y2(6); -(mu*y2(1))/(r^3); -(mu*y2(2))/(r^3); -(mu*y2(3))/(r^3)];
return
end

function ydot3 = two_body3(~,y3)
mu = 398600;
r = norm(y3(1:3));
ydot3 = [y3(4); y3(5); y3(6); -(mu*y3(1))/(r^3); -(mu*y3(2))/(r^3); -(mu*y3(3))/(r^3)];
return
end

function [value, isterminal, direction] = events2(~, state_out2)
value = state_out2(end,2);
isterminal = 1;
direction = 0;
end

function [value, isterminal, direction] = events3(~, state_out3)
value = (state_out3(end,2));
isterminal = 1;
direction = 0;
end

Here is the script called at the beginning (Astronautics_Project_Calculations)

clc
%INITIAL ORBIT
global mu
mu = 398600;                                        % km^3/s^2
alt1 = 400;
global r1
r1 = 6378 + alt1*1.60934;                           %Convert to km
global v1   
v1 = (mu/r1)^0.5;                                   %Calculate Orbital Velocity
global T1
T1 = ((2*pi)/(mu^0.5))*(r1^(3/2));                  %Calculate Orbital Period
disp 'ORBIT #1'
disp ' '
disp 'Orbital Radius (km) = '
disp (r1)
disp 'Orbital Velocity (km/s) = '
disp (v1)
disp 'Orbital Period (s) = '
disp (T1)
%HOHMANN TRANSFER
alt3 = 22236; 
global r3
r3 = 6378 + alt3*1.60934;                           %Convert to km
h1 = ((2*mu)^0.5)*((r1^2)/(2*r1))^(0.5);            %knowing apoapsis radius and 
h2 = ((2*mu)^0.5)*((r1*r3)/(r1+r3))^(0.5);          %periapsis radius for all 3 orbits, 
h3 = ((2*mu)^0.5)*((r3^2)/(2*r3))^(0.5);            %calculate angular momentum
va1 = h1/r1;                                        %velocity at point A of orbit 1 (same as v1)
global va2
va2 = h2/r1;                                        %velocity at point A (periapsis) of orbit 2
global dv1
dv1 = va2-va1;                                      %delta V necessary to get from orbit 1 to orbit 2
global vb2
vb2 = h2/r3;                                        %velocity at point B (apoapsis) of orbit 2
global vb3
vb3 = h3/r3;                                        %velocty at point B of orbit 3 (circular)
global dv2
dv2 = vb3-vb2;                                      %delta V necessary to get from orbit 2 to orbit 3
global a
a = (r1+r3)/2;                                      %semi-major axis of orbit 2 (r1 = rp , r3 = ra)
global en
en = -(mu/(2*a));                                   %orbital energy of orbit 2
global e
e = (r3-r1)/(r3+r1);                                %eccentricity of orbit 2
global T2
T2 = ((2*pi)/(mu^0.5))*(a^(3/2));                   %period of orbit 2
global TT
TT = T2/2;                                          %transfer time, assuming burn times are negligible, would be half the period of orbit 2
disp 'ORBIT #2'
disp ' ' 
disp 'Periapsis Radius (km) = '
disp (r1)
disp 'Periapsis Velocity (km/s) = '
disp (va2)
disp 'Apoapsis Radius (km) = '
disp (r3)
disp 'Apoapsis Velocity (km/s) = '
disp (vb2)
disp 'Orbital Energy = '
disp (en)
disp 'Semi-Major Axis (km) = '
disp (a)
disp 'Eccentricity = '
disp (e)
disp 'Orbital Period (s) = '
disp (T2)
disp 'Transfer Time (s) = '
disp (TT)
disp 'Delta V1 (periapsis)(km/s) = '
disp (dv1)
disp 'Delta V2 (apoapsis)(km/s) = '
disp (dv2)
global T3
T3 = ((2*pi)/(mu^0.5))*(r3^(3/2));                  %period of orbit 3
%FINAL ORBIT 
disp 'ORBIT #3'
disp ' '
disp 'Orbital Radius (km) = '
disp (r3)
disp 'Orbital Velocity (km/s) = '
disp (vb3)
disp 'Orbital Period (s) = '
disp (T3)
%ORBITAL TRANSFER PROPELLANT BUDGET ANALYSIS
g0 = 9.81;                                          %(m/s^2)
WGSM = 1500;                                        %(kg) Spacecraft Mass
DTSM = 2500;                                        %(kg) Delta Transfer Stage (dry) Mass
TLPM = 6249;                                        %(kg) Total Loaded Propellant Mass
Isp = 450.5;                                        %(s) Specific Impulse
Fs = 100;                                           %(kN) Transfer Stage Average Thrust
M1 = WGSM+DTSM+TLPM;                                %(kg) Total Spacecraft Mass Before 1st Burn
dv1m = dv1*1000;                                    %convert delta v from km/s to m/s for use in eqn 6.1
DeltaM1 = M1*(1-(exp(-dv1m/(Isp*g0))));              %calculate delta m required for 1st burn
disp 'Propellant Mass Required for 1st Burn (kg) = '
disp (DeltaM1)
M2 = M1-DeltaM1;                                    %(kg) Mass of Spacecraft After 1st Burn
dv2m = dv2*1000;                                    %convert delta v from km/s to m/s for use in eqn 6.1
DeltaM2 = M2*(1-(exp(-dv2m/(Isp*g0))));              %calculate delta m required for 2nd burn
disp 'Propellant Mass Required for 2nd Burn (kg) = '
disp (DeltaM2)
Mremaining = TLPM-DeltaM1-DeltaM2;               %calculate remaining propellant mass
disp 'Propellant Mass Remaining to be Discarded (kg) = '
disp (Mremaining);
M3 = M2 + DeltaM2;
b = (Fs/dv1)*log(M1/M2);                            %calculate burn rate from delta v
disp 'Average Propellant Burn Rate (kg/s) = '
disp (b)

Thanks in advance for any help.

[Jan]

The detail "simulation of a Hohmann transfer from a parking orbit to a geosynchronous orbit" is confusing only. The chances to get a useful answer grow, if you concentrate of the required details. Omit e.g. the huge block of code for displaying the results, because this does not matter the problem.

[Christian Harvie]

Thanks, I will be sure to simplify future questions.

Answer 1

[Jan]

The event function is called as

[value,isterminal,direction] = myEventsFcn(t, y)

The 2nd input is the current state vector. Then y(end,2) causes the error, because y is a vector, not a matrix.

You can find this using the debugger easily. Type this in the command window:

dbstop if error

Now start the code again until it stops at the error. Now check the size of y:

size(y)

I guess you want y(6) instead of y(end,2).

[Christian Harvie]

Thank You! I had assumed it was a matrix.