MATLAB: Ode45 dynamics rocket around earth equation of motion

equations of motionode45

Hello,
I am trying to solve these 2 equations of motion for the radius (r) and theta, of the rocket with respect to time. Then I want to plot their x vs y trajectory. For some reason, my theta is stopping at 1.7678. Any help would be greatly appreciated.

Best Answer

You have the wrong initial conditions. They should be:
r0
rdot0
theta0
thetadot0
But you have vtheta0 in that last spot instead of thetadot0. So the value and units are all crazy for that variable, hence the garbage plot. To fix this, I would suggest using variables to help you see what is correct, and comment all of the initial condition calculations with the units used. E.g., with the initial condition correction and the rdot correction and annotated plots:
function orbit_polar
clc;
close all
clear all
timestep = 0:10:80000;
%\
% Case (i) v0 = sqrt(g*Re^2/r0)
%/
Re = 6371; % km


g = .00981; % km/s^2
r0 = 5*Re; % km
v0 = sqrt(g*Re^2/r0); % km/s
rdot0 = 0; % km
theta0 = 0; % rad
thetadot0 = v0 / r0; % rad/s
initial = [r0;rdot0;theta0;thetadot0]; % km, km/s, rad, rad/s
[t,y] = ode45(@ pew,timestep,initial);
r = y(:,1);
theta = y(:,3);
x = r.*1000.*cos(theta); % m

ypos = r.*1000.*sin(theta); % m
figure
plot(x,ypos);
axis square; % changed slightly, but you can use axis equal here as instructed
grid on
xlabel('X (m)');
ylabel('Y (m)');
title('v0 = sqrt(g*Re^2/r0)');
figure 
plot(t,y(:,1));
grid on
xlabel('Time (s)');
ylabel('r (km)');
title('Radius magnitude');
end
function ydot = pew(t,y)
r = y(1) ; rdot=y(2) ; theta = y(3) ; thetadot=y(4);
g = .00981;
Re=6371;
rdoubledot = (r*thetadot^2)-g*(Re/r)^2 ;
thetadoubledot = -(2*rdot*thetadot)/r; % fixed rdot
ydot = [rdot;rdoubledot;thetadot;thetadoubledot];
end
Related Question