Plot the orbit of a satellite

569 visualizaciones (últimos 30 días)
Bora K
Bora K el 9 de En. de 2020
Respondida: Christian Bellinazzi el 13 de Nov. de 2022
Hello everyone,
I have this MATLAB function satellit(t,x,model) provides the system of differential equations for the orbit elements x = (a e i O w M) of a satellite in an Earth orbit.
I have to write a function to compute the orbit of a satellite. In this function the user should provide model initial position x(0) and flight time tf. Finally I have to plot the satellite orbit in cartesian coordinates. However I'm struggling to do it. Could anyone help me ?
function dx = satellit(t,x,model)
dx = zeros(6,1);
a = x(1); % semi-major axis [km]
e = x(2); % eccentricity
i = x(3); % inclination [rad]
O = x(4); % longitude of the ascending node [rad]
w = x(5); % argument of periapsis [rad]
M = x(6); % mean anomaly [rad]
GM = 398600.4418; % graviational parameter in km^3/s^2
R = 6378.1370; % radius at equator in km
J2 = 0.0010826267d0; % Earth's J2 (WGS-84)
J3 = -0.0000025327d0; % Earth's J3 (WGS-84)
J4 = -0.0000016196d0; % Earth's J4 (WGS-84)
ac = a*a*a;
n = sqrt(GM/ac);
Thanks in advance

Respuesta aceptada

Meg Noah
Meg Noah el 9 de En. de 2020
Editada: Meg Noah el 28 de Oct. de 2021
Edited 10/28/2021 - Thanks to a bug found by @Frank Epstein in the interpretation of Mdot - this is updated in TLE2OrbitalElements.m
Here you go. The attached scripts produce these plots for the International Space Station. You can supply your own TLE for the satellite of interest. This does not include atmospheric drag or solar impacts. Would you like that as well - the special perturbations applied over time? Or did you just need the basic Kepler COE -> Cartesian transformations?
  35 comentarios
Meg Noah
Meg Noah el 28 de Oct. de 2021
Editada: Meg Noah el 28 de Oct. de 2021
@Frank Epstein Thank you for pointing out the bug in the Mdot interpretation - I edited and uploaded the code. I appreciate you taking the time to alert me and the user community.
86400 is the number of seconds in a day. This converts orbits per day (given by line2(53:63)) to orbits per second. The 2*pi converts orbits to radians as one orbit has 2*pi radians. https://en.wikipedia.org/wiki/Two-line_element_set
Walter Roberson
Walter Roberson el 29 de Oct. de 2021
Meg, thank you for your continued contributions on these topics!

Iniciar sesión para comentar.

Más respuestas (3)

MANISH SINGH
MANISH SINGH el 19 de Feb. de 2022
function [P] = Orbital_period ( R )
% orbital period ( R ) calculates the orbital period of the satellite orbitting
% around the earth in seconds for radius of orbit R in meters, R is the
% total Radius of earth and height of satellite from earth
% To call this function type[ P ] = orbital period ( R )
Ge = 6.673e-11; %Earth Gravitational constant in N*m^2/kg^2;
Me = 5.98e24; %Earth Mass in Kg
P = sqrt((4*(pi^2)*R^3)/(Ge*Me));
end

Christian Bellinazzi
Christian Bellinazzi el 13 de Nov. de 2022
clc
clear all
%% Definizione dei parametri del problema
G = 66.7 * 10^(-12); % [m^3/(Kg * s^2)]
M = 5.972 * 10^(24); % [Kg]
mu = (G * M);
a = 7*10^6; % semiasse maggiore dell'ellisse [km]
tpK = linspace(0,5830,1000);
zenith_angle = 89*(pi/180); %[rad]
r_E = 6378; %Earth radius [km]
r_M=1738; % Moon radius [Km]
phi = pi/2-zenith_angle;
% z = input('Inserire la quota desiderata in [Km]: ');
% v_BO = input('Inserire velocità al burn-out, inferiore a 11.2 [Km/s]: ');
% e = input('Inserire il valore di eccentricità dell orbita [a]: ');
% tpK_mio = input('Inserire in quale istante di tempo si vuole visualizzare il satellite, [0:5800] [s]: ');
% inclinazione = input('Inserire angolo eclittica-equatoriale, compreso tra 0 e 180 gradi: ');
z = 1000;
v_BO=10;
e=0.1;
tpK_mio=2000;
inclinazione=5;
r_BO = r_E+z;
h = r_BO*v_BO*cos(phi);
%% Problema di Keplero generalizzato all'intervallo di tempo
E1 = tpK * sqrt(mu/(a^3));
n_iter = length(E1);
tol = 10^(-8);
E_1 = zeros(1,n_iter);
f1 =@(E) E - e*sin(E);
f2 =@(E) 1 - e*cos(E);
kk = zeros(1,n_iter);
nuK_1 = zeros(1,n_iter);
for P=1:n_iter
E_1(P) = E1(P);
for i=2:n_iter
a1 = E_1(i-1);
f1_NEW = f1(a1);
f2_NEW = f2(a1);
E_1(i) = a1 - (f1_NEW-E_1(P))/f2_NEW;
if E_1(i) - a1 < tol
kk(P) = E_1(i) * 180/pi;
b_1 = E_1(i);
if kk(P)>=0 && kk(P)<180
nuK_1(P) = acos(-(e - cos(b_1))/(1 - e * cos(b_1))) * 180/pi;
i=n_iter+1;
else
nuK_1(P) = -acos(-(e - cos(b_1))/(1 - e * cos(b_1))) * 180/pi;
i=n_iter+1;
end
end
end
end
%% Problema di Keplero per il punto considerato all'intervallo specifico
E1_mio = tpK_mio * sqrt(mu/(a^3));
n_iter = 3;
f1_mio =@(E) E - e*sin(E);
f2_mio =@(E) 1 - e*cos(E);
E_mio = zeros(1,n_iter);
E_mio(1) = E1_mio;
for i=2:n_iter
a1_mio = E_mio(i-1);
f1_mio_new = f1_mio(a1_mio);
f2_mio_new = f2_mio(a1_mio);
E_mio(i) = a1_mio - (f1_mio_new-E1_mio)/f2_mio_new;
if E_mio(i) - a1_mio < tol
disp('Il valore di Anomalia Eccentrica selezionato risulta essere [deg]:');
kk_mio = E_mio(i) * 180/pi
b_1_mio = E_mio(i);
if kk_mio>=0 && kk_mio<180
disp('Il valore di Anomalia Vera selezionato risulta essere [deg]:');
nuK_1_mio = acos(-(e - cos(b_1_mio))/(1 - e * cos(b_1_mio))) * 180/pi
else
disp('Il valore di Anomalia Vera selezionato risulta essere [deg]:');
nuK_1_mio = -acos(-(e - cos(b_1_mio))/(1 - e * cos(b_1_mio))) * 180/pi
end
end
end
%% Perifocal Reference System
T=2*pi*sqrt(a^3/mu); %Orbital Period
t=linspace(0,T,length(E_1)); %Time Discretization
nu=zeros(length(t),1); %True anomaly discretization (initialization)
r=zeros(length(t),1); %Position vector discretization
r_P=zeros(length(t),1); %Position vector components (initialization)
r_Q=zeros(length(t),1);
r_W=zeros(length(t),1);
r_PQW=zeros(length(t),3);
alpha=zeros(length(t),1);
%% Studio dell'inclinazione dell'orbita
for i=1:length(t)
r(i)=((h^2)/mu)/(1+e*cosd(nuK_1(i)))*10^9;
r_P(i)=r(i)*cosd(nuK_1(i));
r_Q(i)=r(i)*sind(nuK_1(i));
if inclinazione>=0 && inclinazione<=90
% primo quadrante deve essere positivo
if nuK_1(i)>=0 && nuK_1(i)<=90
alpha(i) = +inclinazione -(inclinazione/10)*(nuK_1(i))/9;
r_W(i) = r(i)*sind(alpha(i));
end
% secondo quadrante deve essere negativo
if nuK_1(i)>=+90 && nuK_1(i)<=180
alpha(i) = inclinazione -(inclinazione/10)*(nuK_1(i))/9;
r_W(i)=r(i)*sind(alpha(i));
end
% terzo quadrante deve essere negativo
if nuK_1(i)>=-180 && nuK_1(i)<=-90
alpha(i) = inclinazione -(inclinazione/10)*(-nuK_1(i))/9;
r_W(i)=r(i)*sind(alpha(i));
end
% quarto quadrante deve essere positivo
if nuK_1(i)>=-90 && nuK_1(i)<=0
alpha(i) = inclinazione -(inclinazione/10)*(-nuK_1(i))/9;
r_W(i)=r(i)*sind(alpha(i));
end
r_PQW(i,:)=[r_P(i);r_Q(i);r_W(i)];
end
end
% if inclinazione>+90 && inclinazione<+180
% % primo quadrante deve essere negativo
% if nuK_1(i)>0 && nuK_1(i)<90
% alpha = -inclinazione +2*(nuK_1(i))/9;
% r_W(i)=r(i)*sind(alpha);
% end
% % secondo quadrante deve essere positivo
% if nuK_1(i)>+90 && nuK_1(i)<180
% alpha = -inclinazione -2*(-nuK_1(i))/9;
% r_W(i)=r(i)*sind(alpha);
% end
% % terzo quadrante deve essere positivo
% if nuK_1(i)<-90 && nuK_1(i)>-180
% alpha = -inclinazione -2*(nuK_1(i))/9;
% r_W(i)=r(i)*sind(alpha);
% end
% % quarto quadrante deve essere negativo
% if nuK_1(i)>-90 && nuK_1(i)<0
% alpha = -inclinazione -2*(nuK_1(i))/9;
% r_W(i)=r(i)*sind(alpha);
% end
% end
% r_PQW(i,:)=[r_P(i);r_Q(i);r_W(i)];
%
r_mio = ((h^2)/mu)/(1+e*cosd(nuK_1_mio))*10^9;
r_P_mio = r_mio*cosd(nuK_1_mio);
r_Q_mio = r_mio*sind(nuK_1_mio);
if inclinazione>=0 && inclinazione<=90
% primo quadrante deve essere positivo
if nuK_1_mio>=0 && nuK_1_mio<=90
alpha = inclinazione -(inclinazione/10)*(nuK_1_mio)/9;
r_W_mio=r_mio*sind(alpha);
end
% secondo quadrante deve essere negativo
if nuK_1_mio>+90 && nuK_1_mio<=180
alpha = inclinazione -(inclinazione/10)*(nuK_1_mio)/9;
r_W_mio=r_mio*sind(alpha);
end
% terzo quadrante deve essere negativo
if nuK_1_mio>-180 && nuK_1_mio<=-90
alpha = inclinazione -(inclinazione/10)*(-nuK_1_mio)/9;
r_W_mio=r_mio*sind(alpha);
end
% quarto quadrante deve essere positivo
if nuK_1_mio>-90 && nuK_1_mio<0
alpha = inclinazione -(inclinazione/10)*(-nuK_1_mio)/9;
r_W_mio=r_mio*sind(alpha);
end
end
r_PQW_mio = [r_P_mio,r_Q_mio,r_W_mio];
%% Definizione della superficie della Terra
x_R = linspace(-r_E,r_E,4000);
x_M = linspace(-r_M,r_M,4000);
funct1 = sqrt(r_E^2-x_R.^2);
funct2 = -sqrt(r_E^2-x_R.^2);
funct3 = sqrt(r_M^2-x_M.^2);
funct4 = -sqrt(r_M^2-x_M.^2);
%% Rappresentazione nel piano
figure(1)
plot(r_P(:,1),r_Q(:,1),'linewidth',1.5)
hold on
plot(x_R,funct1,x_R,funct2,LineWidth=2,Color='c')
hold on
plot(38460+x_M,funct3,38460+x_M,funct4,LineWidth=2,Color='b')
hold on
plot(r_P_mio,r_Q_mio,'r*',LineWidth=2)
hold on
plot(0,0,'*',LineWidth=2)
hold on
plot(38460,0,'*',LineWidth=2)
legend('Plane Orbit','Earth surface','','Moon surface','','Design point','Earth center','Moon center')
xlabel('X-coordinate')
ylabel('Y-coordinate')
title('Orbit Perifocal Plane')
axis equal
grid minor
[x1,y1,z1]=sphere(100);
E=imread('earth.jpg');
x_E=x1*r_E;
y_E=y1*r_E;
z_E=z1*r_E;
figure(2)
hold on
plot3(r_PQW(:,1),r_PQW(:,2),r_PQW(:,3),'g',LineWidth=2)
hold on
plot3(r_P_mio,r_Q_mio,r_W_mio,'r*',LineWidth=2)
hold on
warp(x_E,y_E,z_E,E)
grid minor
xlabel('asse x')
ylabel('asse y')
zlabel('asse z')
I have some problems about the inclination of the orbit. Anyone to help me? thanks

Christian Bellinazzi
Christian Bellinazzi el 13 de Nov. de 2022
I would like also to use getframe for tthe moving marker on the orbit. Anyone could hel me?

Categorías

Más información sobre Satellite Mission Analysis en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by