global earthRADIUS gravity
gravity = 0.00981;
earthRADIUS = 6371;
a = 0;
b = 800e3;
c= 80e3;
d = 5;
t_Vec = linspace(a,b,c);
r_0 = (d)*(earthRADIUS);
Vr0 = a;
r_d = Vr0;
theta = a;
% Case 1
vTHETA = sqrt(gravity*earthRADIUS^2/r_0);
thetaDOT = vTHETA/r_0;
p1 = [r_0,r_d,theta,thetaDOT];
% Case 2
vTHETA = sqrt(1.8*gravity*earthRADIUS^2/r_0);
thetaDOT = vTHETA/r_0;
p2 = [r_0,r_d,theta,thetaDOT];
% Case 3
vTHETA = sqrt(2*gravity*earthRADIUS^2/r_0);
thetaDOT = vTHETA/r_0;
p3 = [r_0,r_d,theta,thetaDOT];
% Case 4
vTHETA = sqrt(3*gravity*earthRADIUS^2/r_0);
thetaDOT = vTHETA/r_0;
p4 = [r_0,r_d,theta,thetaDOT];
% Solve ODE
[t1,S1] = ode45(@rODE , t_Vec, p1);
[t2,S2] = ode45(@rODE , t_Vec, p2);
[t3,S3] = ode45(@rODE , t_Vec, p3);
[t4,S4] = ode45(@rODE , t_Vec, p4);
% Extracting ri from each Si
r1 = S1(:,1);
r2 = S2(:,1);
r3 = S3(:,1);
r4 = S4(:,1);
% Similarly extracting theta_i from each Si,
our assumption theta = S(3)
the_1 = S1(:,3);
the_2 = S2(:,3);
the_3 = S3(:,3);
the_4 = S4(:,3);
% x & y vectors from r vectors & theta
vectors
x1 = r1.*cos(the_1);
y1 = r1.*sin(the_1);
x2 = r2.*cos(the_2);
y2 = r2.*sin(the_2);
x3 = r3.*cos(the_3);
y3 = r3.*sin(the_3);
x4 = r4.*cos(the_4);
y4 = r4.*sin(the_4);
rdot1 = S1(:,2);
rdot2 = S2(:,2);
rdot3 = S3(:,2);
rdot4 = S4(:,2);
Vr1 = rdot1;
Vr2 = rdot2;
Vr3 = rdot3;
Vr4 = rdot4;
% For Vtheta_i's we require thetadot_i's
tD_1 = S1(:,4);
tD_2 = S2(:,4);
tD_3 = S3(:,4);
tD_4 = S4(:,4);
vt_1 = r1.*(tD_1);
vt_2 = r2.*(tD_2);
vt_3 = r3.*(tD_3);
vt_4 = r4.*(tD_4);
%--- Plot For Each Case ---
time_set = {t1,t2,t3,t4};
x_set = {x1,x2,x3,x4};
y_set = {y1,y2,y3,y4};
Vr_set = {Vr1,Vr2,Vr3,Vr4};
Vtheta_set = {vt_1,vt_2,vt_3,vt_4};
titles = {'[Case 1]','[Case 2]', '[Case
3]','[Case 4]'};
for i = 1:4
figure(i)
% Draw Trajectory
subplot(2,2,[1,2])
plot(x_set{i},y_set{i}, ); '-b'
xlabel('x (km)');
ylabel('y (km)');
title('Trajectory');
axis equal; grid on;
% Draw vr
subplot(2,2,3)
plot(time_set{i},Vr_set{i}, ); '-m'
title('V_r vs time');
xlabel('t (s)');
ylabel('V_r [km/s]');
grid on
% Draw V theta
subplot(2,2,4)
plot(time_set{i},Vtheta_set{i}, ); '-k'
title('V_\theta vs. time');
xlabel('t (s)');
ylabel('V_\theta [rad/s]');
grid on
% Set title for figure
sgtitle(titles{i});
end
function sD = rODE(t,k)
global Re g
r = k(1);
rdot = k(2);
theta = k(3);
thetaD = k(4);
rDD = (r*thetaD^2)-g*(Re/r)^2 ;
thetaDD = -(2*rdot*thetaD)/r;
sD = [rdot;rDD;thetaD;thetaDD];
end
