ode45 dynamics rocket around earth equation of motion

27 views (last 30 days)
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.
  11 Comments

Sign in to comment.

Accepted Answer

James Tursa
James Tursa on 3 Mar 2020
Edited: James Tursa on 3 Mar 2020
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
  1 Comment
Jackson Hager
Jackson Hager on 3 Mar 2020
Thank you so much!! The answer was on my paper the whole time... you are a life saver.

Sign in to comment.

More Answers (1)

Jackson Hager
Jackson Hager on 3 Mar 2020
Eom1: r’’ + g(Re/r)^2 - r(theta’)^2=0 Re is the radius of the earth in km: 6371 km Eom2: r(theta’’)+ 2(r’)(theta’)=0
Mass canceled when solving for the eom’s. It was also not given in the problem statement.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!