MATLAB: Plotting Help (rocket motion around a planet

MATLABode45plottingrocket motion

Hello! Currently I am trying to model the flight path of a rocket around the planet at various given velocities, however whenever I go to plot it just leaves me with a straight line, I know for a fact with the given conditions I have it should be a circle so I'm just wondering where I'm going wrong, and as always I appreciate the help and thank you for your time.

timespan = 0:0.0008:10;
Re = 6371;
IC = [5*Re 0 0 (sqrt((g*Re^2)/5*Re)/5*Re)];
options = odeset('reltol',1e-10,'abstol',1e-10);
g = 9.81;
[t,y] = ode45(@(t,y) planetmotion(t,y,g),timespan,IC,options)
r =y(:,1);
theta = y(:,3)
x = r.*cos(theta); y= r.*sin(theta);
                
plot(x,y);
axis equal ; xlabel('x,m'); ylabel('y,m');
function ydot = planetmotion(t,y,g)
Re = 6371;
r=y(1); rdot=y(2); theta=y(3); thetadot=y(4);
rddot = r*thetadot^2-g*(Re/r)^2;
thetaddot = -2*rdot*thetadot/r;
ydot = [rdot ; rddot ; thetadot ; thetaddot];
end

From how I've tweaked it so far I think my problem is within my initial condtions particularly the last component, but a fresh set of eyes is always helpful.

Best Answer

Figured out problem,it was the last input for the initial condition like I thought.

Related Solutions

MATLAB: Ode45 dynamics rocket around earth equation of motion

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

MATLAB: Is ode45 changing the value of the state variable between calls to the dynamics function

ode45 implements the Dormand-Prince method which is an embedded 4th / 5th order Runge-Kutta method. For every "step" in the solution, ode45 needs to execute your function six times; i.e. there are six "stages" per every solution "step". Each stage will involve a change in the time and/or the state values.

So, at the end of one step, what you have is the result of a combibnation of all six stages, not just the last call to your function.

Related Question