ode45 dynamics rocket around earth equation of motion

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.

[Jackson Hager]

clc;
close all
clear all
timestep = 0:10:80000;
initial = [31855;0;0;3.53552];


[t,y] = ode45(@ pew,timestep,initial)

r = y(:,1);
theta = y(:,3);
x = r.*1000.*cos(theta);
ypos = r.*1000.*sin(theta);
figure
axis equal ; plot(x,ypos);

figure 
axis equal; plot(t,y(:,1))

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*r*thetadot)/r;
ydot = [rdot;rdoubledot;thetadot;thetadoubledot];
end

[Jackson Hager]

[darova]

Show your equations

Where is the mass?

[darova]

[James Tursa]

+1 darova, this should be an Answer

[darova]

Eveything looks correct. Do you have any data? Maybe a picture?

[Jackson Hager]

[Jackson Hager]

that is the problem statement

[Jackson Hager]

Darova, yes that should have been rdot, thank you. It is still not running properly though, I should be getting a circle for the x and y trajectory but I am getting a linear graph.

[darova]

Do you have original equations? How do they look like?

[Jackson Hager]

Answer 1

[James Tursa]

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

[Jackson Hager]

Thank you so much!! The answer was on my paper the whole time... you are a life saver.

Answer 2

[Jackson Hager]

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.