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 allclose allclcAstronautics_Project_Calculationsformat long;global r1global v1global T1global r3%global va2
global dv1global 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);figuresurf(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')figuresubplot(3,1,1)hold ontitle('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 onxlabel('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 onxlabel('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)];returnendfunction 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)];returnendfunction 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)];returnendfunction [value, isterminal, direction] = events2(~, state_out2)value = state_out2(end,2);isterminal = 1;direction = 0;endfunction [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 mumu = 398600; % km^3/s^2
alt1 = 400;global r1r1 = 6378 + alt1*1.60934; %Convert to km
global v1 v1 = (mu/r1)^0.5; %Calculate Orbital Velocity
global T1T1 = ((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 r3r3 = 6378 + alt3*1.60934; %Convert to kmh1 = ((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 va2va2 = h2/r1; %velocity at point A (periapsis) of orbit 2
global dv1dv1 = va2-va1; %delta V necessary to get from orbit 1 to orbit 2
global vb2vb2 = h2/r3; %velocity at point B (apoapsis) of orbit 2
global vb3vb3 = h3/r3; %velocty at point B of orbit 3 (circular)
global dv2dv2 = vb3-vb2; %delta V necessary to get from orbit 2 to orbit 3
global aa = (r1+r3)/2; %semi-major axis of orbit 2 (r1 = rp , r3 = ra)
global enen = -(mu/(2*a)); %orbital energy of orbit 2
global ee = (r3-r1)/(r3+r1); %eccentricity of orbit 2
global T2T2 = ((2*pi)/(mu^0.5))*(a^(3/2)); %period of orbit 2
global TTTT = 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 T3T3 = ((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.1DeltaM2 = 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.
Best Answer