Help plotting a circular orbit

63 views (last 30 days)
Hi so basically i would like to be able to produce a circular graph of sattelite position around earth. The code attached below starts with acceleration and use Eulers method to integrate once for velocity and again for displacement. The code works and produces results i am just unsure how to convert the outputs into a circular graph that updates displacement around earth as a result of the calculated velocity.
Code:
clear all
clc
G = 6.6743*10^-11; %Gravitational Constant, Units: m^3 kg^-1 s^-2
Mc = 5.972*10^24; %Mass of central body (earth), Units: kg
rE = 6371; %Radius of earth, Units: Km
height = 4000 %Altitude of sattelites orbit, Units: km
r = height + rE; %Radius of circular orbits, Units: km
mu = 3.986004418*10^5; %Earths Gravitational parameter, Units: km^3s^-2
x_initial = 0;
v_initial = 12000; %km/h
dt = 0.1;
t_vector = 0:dt:100000;
x(1) = x_initial;
v(1) = v_initial;
%%Math working out
%Fx = -G*x*(Mc/r^3);
%a = Fx/Mc
%a = (-G*x)/r^3
%a = -5.9833e-23*x m/s^2
%v(t) = 12000 + a*t
%y(t) = r + 12000t + 0.5a*t^2
x_initial = 0;
v_initial = 12000;
x(1)=x_initial;
v(1)=v_initial;
for i=1:length(t_vector)-1
a = -G/r^3;
v(i+1) = v(i) + a*dt;
xslope = v(i);
x(i+1) = x(i)+xslope*dt;
end
figure(1)
plot(t_vector,x)
figure(2)
plot(t_vector,v)
  1 Comment
Walter Roberson
Walter Roberson on 4 Jun 2022
You need to track x and y separately.

Sign in to comment.

Accepted Answer

Lateef Adewale Kareem
Lateef Adewale Kareem on 4 Jun 2022
G = 6.6743*10^-11; %Gravitational Constant, Units: m^3 kg^-1 s^-2
Mc = 5.972*10^24; %Mass of central body (earth), Units: kg
rE = 6371; %Radius of earth, Units: Km
height = 4000; %Altitude of sattelites orbit, Units: km
r = height + rE; %Radius of circular orbits, Units: km
x_initial = 0;
y_initial = r*1e3;
dt = 1;
T = sqrt(4*pi^2*(r*1e3)^3/(G*Mc)); % period in sec
t_vector = 0:dt:T;
dsdt = @(s)[s(3), s(4), -G*Mc*s(1)/norm(s([1,2]))^3, -G*Mc*s(2)/norm(s([1,2]))^3]; %satelite orbitaal dynamics
x = x_initial; y = y_initial; v = sqrt(G*Mc/(r*1e3));
xv = v*y_initial/(r*1e3); yv = -v*x_initial/(r*1e3); % initial velocities
state = [x_initial, y_initial, xv, yv];
for i=2:length(t_vector)
state = state + rk4(dsdt, state, dt);
x = [x, state(1)];
y = [y, state(2)];
xv = [xv, state(3)];
yv = [yv, state(4)];
end
figure
subplot(2,2,1)
plot(t_vector,x)
title('x position of satellite')
subplot(2,2,2)
plot(t_vector,y)
title('y position of satellite')
subplot(2,2,3)
plot(t_vector,xv)
title('x velocity of satellite')
subplot(2,2,4)
plot(t_vector,yv)
title('y velocity of satellite')
figure
plot(x, y)
title('orbit of satellite')
function dy = rk4(dydt, y, dt) %runge kutta integrator
k1 = dydt(y); k2 = dydt(y + dt*k1/2);
k3 = dydt(y + dt*k2/2); k4 = dydt(y + dt*k3);
dy = dt*(k1+2*k2+2*k3+k4)/6;
end
  1 Comment
Adam Lacey
Adam Lacey on 4 Jun 2022
This is amazing thank you so much

Sign in to comment.

More Answers (0)

Categories

Find more on Reference Applications in Help Center and File Exchange

Products


Release

R2021b

Community Treasure Hunt

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

Start Hunting!