Plot the orbit of a satellite

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

Answer 1

[Meg Noah]

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?

[Bora K]

Thank you for the very good answer. I don't need to define the perturbations. However TLE files are a bit confusing because I want to see the behaviour of the satellite orbit for varying excentricty, semi major axis or inclination.

For example I have an initial position which is x(0) = (38000 km e 50 0 0 0) and excentricity is varying from 0.0001 to 0.6. Do you know how to examine this case?

Thanks in advance

[Meg Noah]

I posted the project here:

https://www.mathworks.com/matlabcentral/fileexchange/73873-satellite-orbit-coordinate-transformations

For just your own Orbital Elements you can use the functions provided at the above link (download it and rate it please) and this script. The script can easily be edited with your Orbital Elements.

%% *OrbitalElements2OrbitStateVector*
%% *Purpose*
%  To show how to convert the orbital elements to the inertial p-q orbital
%  plane.  The p-axis is through the center of the orbit to perigee.  The q
%  axis is through the focus (center of Earth) and normal to the p axis.
%% *History*
%  When       Who    What
%  ---------- ------ --------------------------------------------------
%  2020/01/09 mnoah  original code

%% *Standard Gravitational Parameter*
% The standard gravitational parameter $$\mu$$ of a celestial body is the
% product of the gravitational constant G and the mass M of the body. 
mu = 3.98618e14; % [m3/s2] Earth's geocentric gravitational constant

%% *Express your own Orbital Elements*
e = 0; % eccentricity
x = [38000 e 50 0 0 0]; 
OE.a_km = x(1);  % semi-major axis [km]
OE.e = x(2);  % eccentricity
OE.i_deg = rad2deg(x(3));  % inclination [deg]
OE.omega_deg = rad2deg(x(4));  % longitude of the ascending node [deg]
OE.Omega_deg = rad2deg(x(5));  % argument of periapsis [deg]
OE.M_deg = rad2deg(x(6));  % mean anomaly [deg]

fprintf(1,['Kepler Elements for satelliteID %d epoch %s:\n' ...
    '\ta [m] = %f\n\te = %f\n\ti [deg] = %f\n\tomega [deg] = %f\n' ...
    '\tOmega [deg] = %f\n\tM [deg] = %f\n'], floor(OE.satelliteID), ...
    datestr(OE.epoch),OE.a_km, OE.e, OE.i_deg, OE.omega_deg, ...
    OE.Omega_deg, OE.M_deg);
a_m = OE.a_km*1e3;
e = OE.e;
M_deg = OE.M_deg;

%% *Orbital Plane Coordinates*
% In the plane of the orbit:
%  
% $$p= a\left ( \cos E - e \right )$$
%  
% $$q=a\sqrt{1-e^{2}}\sin E$$
% 
% where
%  
% $$p$ - [m] coordinate along axis through center and perigee
%  
% $$q$ - [m] coordinate along axis passing through focus and perpendicular
% to the p-axis
%  
% $$E$ - [rad] eccentric anomaly
%  
% $$e$ - [unitless] eccentricity
%% *Mean Motion*
% The mean motion, n, is the angular speed required for a body to 
% complete one orbit = 2*pi/Period
n = sqrt(mu/a_m^3);  % [rad/s] mean motion
%% *Derivatives Yeild Velocity Components Along Orbit*
%
% $$\dot M = \dot E - e (\cos E) \dot E \\ = n$$
%  
% $$\dot E = \dot M / (1 - e \cos E) \\$$
%  
% $$\dot P = -a (\sin E) \dot E \qquad $$
%  
% $$\dot Q = a (\cos E) \dot E \sqrt{1 - e^2}$$
%
%% *Use Kepler Equation to Find Location Along Orbit*
% $$E-e\sin \left ( E \right )= n\left ( t-t_{p} \right )=M$$
%  
% where
%  
% $$t_{p}$ - [s] time of periapsis
%  
% $$n$ - [rad/s] mean motion (in the TLE it is orbits/day)
%  
% $$M$ - [rad] mean anomaly
%% *Use Newton Method to Solve the Kepler Equation*
% $$E_{i+1} = E_i - f(E_i) / f'(E_i) \\$$
%  
% $$f'(E) = 1 - e \cos E$$
%% *Output Orbit Plane State Vector*
%  p_m - [m] coordinate along axis through center and perigee
%  q_m - [m] coordinate passing through focus and perpendicular to p-axis
%  dpdt_m_per_s = [rad/s] p component velocity
%  dqdt_m_per_s = [rad/s] q component velocity
M_rad = deg2rad(M_deg);
E_rad = M_rad; 
dE = 99999;
eps = 1e-6; % [rad] control precision of Newton's method solution
while (abs(dE) > eps)
    dE = (E_rad - e * sin(E_rad) - M_rad)/(1 - e * cos(E_rad));
    E_rad = E_rad -  dE;
end
p_m = a_m*(cos(E_rad) - e);
q_m = a_m*sqrt(1 - e^2)*sin(E_rad);

dMdt_rad_per_s = n;
dEdt_rad_per_s = dMdt_rad_per_s/(1 - e*cos(E_rad));
dpdt_m_per_s = -a_m*sin(E_rad)*dEdt_rad_per_s;
dqdt_m_per_s = a_m*cos(E_rad)*dEdt_rad_per_s*sqrt(1 - e^2);

%% *Plot the Plane of the Orbit*
ra_m = a_m*(1 + e);  % [m] apogee 
rp_m = a_m*(1 - e);  % [m] perigee
rEarth_m = 6378e3; % [m] Earth radius at equator
Evals = 0:0.01:360.0; % [deg] values of the eccentric anomaly around orbit 
pvals = a_m*(cosd(Evals)-e); % [m] orbit positions
qvals = a_m*sqrt(1 - e^2)*sind(Evals); % [m] orbit positions

figure('color','white');
fill(rEarth_m.*cosd(Evals),rEarth_m.*sind(Evals),[0.75 1.00 0.75]);
hold on;
plot(pvals,qvals,'.b');
plot(p_m,q_m,'.r','MarkerSize',25);
p1_m = p_m+dpdt_m_per_s*300;
q1_m = q_m+dqdt_m_per_s*300;
plot([p_m p1_m],[q_m q1_m],'r');
theta = -atan2d(dqdt_m_per_s,dpdt_m_per_s);
plot(p1_m+[0 -0.1*rEarth_m*cosd(theta+30)],q1_m+[0  0.1*rEarth_m*sind(theta+30)],'r');
plot(p1_m+[0 -0.1*rEarth_m*cosd(theta-30)],q1_m+[0  0.1*rEarth_m*sind(theta-30)],'r');

% axes
plot([-1.25*rp_m 1.25*rp_m],[0 0],'k');
plot([0 0],[-1.25*ra_m 1.25*ra_m],'k');
text(1.25*rp_m, 0.1*rEarth_m, 'p');
plot(1.25*rp_m+[0 -0.1*rEarth_m*cosd(30)],[0 0.1*rEarth_m*sind(30)],'k');
plot(1.25*rp_m+[0 -0.1*rEarth_m*cosd(30)],[0 -0.1*rEarth_m*sind(30)],'k');
text(0.1*rEarth_m, 1.25*ra_m, 'q');
plot([0 -0.1*rEarth_m*sind(30)],1.25*ra_m+[0 -0.1*rEarth_m*cosd(30)],'k');
plot([0 +0.1*rEarth_m*sind(30)],1.25*ra_m+[0 -0.1*rEarth_m*cosd(30)],'k');

axis equal
axis tight
axis off
box on
title({'Initial State Vector of Satellite in Orbital Plane'; ...
    ['TLE Epoch: ' datestr(OE.epoch)]});


%% *Rotate To ECI*
Rz_Omega = [ ...
    [cosd(OE.Omega_deg) sind(OE.Omega_deg) 0]; ...
    [-sind(OE.Omega_deg) cosd(OE.Omega_deg) 0]; ...
    [0 0 1]];
Rx_i = [ ...
    [1 0 0]; ...
    [0 cosd(OE.i_deg) sind(OE.i_deg)]; ...
    [0 -sind(OE.i_deg) cosd(OE.i_deg)]];
Rz_omega = [ ...
    [cosd(OE.omega_deg) sind(OE.omega_deg) 0]; ...
    [-sind(OE.omega_deg) cosd(OE.omega_deg) 0]; ...
    [0 0 1]];

% time of epoch
[Year,Month,Day,H,M,S] = datevec(OE.epoch);
HourUTC = H + M/60.0 + S/3600.0;
jd = juliandate(Year,Month,Day,HourUTC,0,0);
jd0 = juliandate(Year,Month,Day,0,0,0);
% form time in Julian centuries from J2000
T = (jd - 2451545.0d0)./36525.0d0;
% D0 = (jd0 - 2451545.0d0);
% % [deg] GMST = Sidereal Time at Greenwich
% GMST = 6.697374558 + 0.06570982441908*D0 + 1.00273790935*HourUTC + 0.000026*T^2;
% % [deg] Sidereal Time
% ST_deg = 100.46061837 + 36000.770053608*T + 0.000387933*T^2 - T^3/38710000;
% GST_deg = ST_deg + 1.00273790935*HourUTC*15;
% %LST_deg = GST_deg - ObsLongitude_degW;

Theta_deg = 100.460618375 + 36000.770053608336*T + 0.0003879333*T^2 + 15*H + M/4 + mod(S/240,360);

Rz_hour = [ ...
    [cosd(Theta_deg) sind(Theta_deg) 0]; ...
    [-sind(Theta_deg) cosd(Theta_deg) 0]; ...
    [0 0 1]];

% position of satellite at epoch in the orbit pq plane
r_pq = [p_m q_m 0]';
% position of satellite at epoch in ECI coordinates

omega_deg = OE.omega_deg;
Omega_deg = OE.Omega_deg;
i_deg = OE.i_deg;

%{
px = cosd(omega_deg)*cosd(Omega_deg) - sind(omega_deg)*cosd(i_deg)*sind(Omega_deg);
py = cosd(omega_deg)*sind(Omega_deg) + sind(omega_deg)*cosd(i_deg)*cosd(Omega_deg);
pz = sind(omega_deg)*sind(i_deg);

qx = -sind(omega_deg)*cosd(Omega_deg) - cosd(omega_deg)*cosd(i_deg)*sind(Omega_deg);
qy = -sind(omega_deg)*sind(Omega_deg) + cosd(omega_deg)*cosd(i_deg)*cosd(Omega_deg);
qz = cosd(omega_deg)*sind(i_deg);

wx = sind(i_deg)*cosd(omega_deg);
wy = -sind(i_deg)*sind(Omega_deg);
wz = cosd(i_deg);

r_ECI = [(p_m*px+q_m*qx) (p_m*py+q_m*qy) (p_m*pz+q_m*qz)];
%}

%guess  = [-5107606.49 -1741563.23 -4118273.08]';
r_ECI = inv(Rz_Omega)*inv(Rx_i)*inv(Rz_omega)*r_pq;
% r_ECR = Rz_hour*r_ECI;
% r_LLA = ecef2lla(r_ECI');
r_LLA = eci2lla(r_ECI',datevec(datenum(Year, Month, Day, H, M, S)),'IAU-2000/2006');
% 
% disp(num2str(r_LLA(1)));
% disp(num2str(wrapTo180(atan2d(r_ECI(2),r_ECI(1))-Theta_deg)));
% disp(datestr(OE.epoch));


Evals = 0:1:360.0; % [deg] values of the eccentric anomaly around orbit 
Orbit_p = a_m*(cosd(Evals)-e); % [m] orbit positions
Orbit_q = a_m*sqrt(1 - e^2)*sind(Evals); % [m] orbit positions
deltaT_s = ((Evals-E_deg_epoch) - e*sind(Evals-E_deg_epoch))/n_deg_per_s; % [s] time since epoch along orbit
% r_ECR_orbit = Rz_hour*r_ECI_orbit';
% r_LLA_orbit = ecef2lla(r_ECR_orbit');
% only due one at a time... <sigh>
Orbit_ECI = zeros(numel(deltaT_s),3);
Orbit_LLA = zeros(numel(deltaT_s),3);
for ipt = 1:size(Orbit_ECI,1)
    r_pq = [Orbit_p(ipt) Orbit_q(ipt) 0]';
    Orbit_ECI(ipt,:) = [inv(Rz_Omega)*inv(Rx_i)*inv(Rz_omega)*r_pq]'; %[Rz_Omega*Rx_i*Rz_omega*r_pq]';
    lla = eci2lla(Orbit_ECI(ipt,:),datevec(datenum(Year, Month, Day, H, M, S+deltaT_s(ipt))),'IAU-2000/2006');
    Orbit_LLA(ipt,:) = lla;
end

%% *Plot Cartesian Coordinates*
figure('color','white');
plot3(Orbit_ECI(:,1),Orbit_ECI(:,2),Orbit_ECI(:,3))
xlabel('ECI x [m]');
ylabel('ECI y [m]');
zlabel('ECI z [m]');
title('Satellite Orbit in ECI Coordinates');
grid on

%% *Plot Map*
figure('color','white');
earth = imread('ear0xuu2.jpg');
lv= size(earth,1);
lh= size(earth,2);
lats =  (1:lv)*180/lv - 90;
lons =  (1:lh)*360/lh - 180;
image(lons, -lats, earth)
hold on;
set(gca,'ydir','normal');
grid on
plot(Orbit_LLA(:,2),Orbit_LLA(:,1),'.r');
plot(r_LLA(2),r_LLA(1),'p','MarkerFaceColor',[1 0 0.7],'MarkerEdgeColor','none','Markersize',14);
set(gca,'XTick',[-180:30:180]);
set(gca,'YTick',[-90:30:90]);
set(gca,'Xcolor',0.3*ones(1,3));
set(gca,'Ycolor',0.3*ones(1,3));
title(['Ground Track for Epoch ' datestr(OE.epoch)])
xlabel('Longitude (degrees)')
ylabel('Latitude (degrees)')
axis equal
axis tight
set(gca,'fontweight','bold');

For the values you posted, and with e=0 (circular orbit), it produces these plots:

The script will let you play around with altitude and eccentricity and such. If you need the ECEF (aka ECR) cartesian coordinates, they can easily be uncommented and plotted. Just follow the way it is plotted for ECI.

[Meg Noah]

Was this answer OK? Are you going to accept it, or do you want some additional or different input?

[Bora K]

Thanks a lot this is the answer.

[KiranKumar Makam]

Hi Meg

I have a small question

Suppose I have a TLE of two objects can I find the point of collision or shortest path during their orbit?and plot them. I have been trying to use your codes at

https://in.mathworks.com/matlabcentral/fileexchange/73873-satellite-orbit-coordinate-transformations

and also the code given by to plot however I am failing to plot it on a 3D map or get time in future when both bodies will pass through the closest to each other

Can you help it out. Inputs are random TLEs of derbis floating around earth.

[Florian Anderl]

Thank you so much for these excellent and useful scripts !

[ali hassan]

plz help me coz i am running the same code but it is giving errors

output:

>> getSatelliteTLE

Unrecognized function or variable 'ID'.

Error in getSatelliteTLE (line 17)

if (ID == 43530)

code:

function [TLE] = getSatelliteTLE(ID,inEpochDatenum)
%% *purpose*
% return the TLE for Satellite based on epoch
%% *inputs*
%  ID - Spacecraft ID
%       73027 = Skylab
%  inEpochDatenum - TLEs can be selected based on epoch 
%% *outputs*
%  TLE - the two line element set corresponding to the satellite at that
%        epoch
%% *history*
%  When       Who    What
%  ---------- ------ --------------------------------------------------
%  2019/07/17 mnoah  original code
%  2020/01/19 mnoah  placeholder - edit to put your own TLE parser
%                    depending on your TLE source
if (ID == 43530)
    TLE = { ...
        'ISS (ZARYA)'; 
        '1 43530U 18056B   20340.83568049  .00000187  00000-0  33778-4 0  9995'; 
        '2 43530  97.9845  56.0782 0001443  79.9736 280.1626 14.76934526129992'};
else
    error('ISS Placeholder only - modify code for getting TLE of other spacecraft');
end
end
%% *OrbitalElements2OrbitStateVector*
%% *Purpose*
%  To show how to convert the orbital elements to the inertial p-q orbital
%  plane.  The p-axis is through the center of the orbit to perigee.  The q
%  axis is through the focus (center of Earth) and normal to the p axis.
%% *History*
%  When       Who    What
%  ---------- ------ --------------------------------------------------
%  2019/07/11 mnoah  original code
%% *Standard Gravitational Parameter*
% The standard gravitational parameter $$\mu$$ of a celestial body is the
% product of the gravitational constant G and the mass M of the body. 
mu = 3.98618e14; % [m3/s2] Earth's geocentric gravitational constant
%% *Get the Satellite TLE*
ID = 43530;
[TLE] = getSatelliteTLE(ID);
%% *Convert to Orbital Elements*
[OE] = TLE2OrbitalElements(TLE);
fprintf(1,['Kepler Elements for satelliteID %d epoch %s:\n' ...
    '\ta [m] = %f\n\te = %f\n\ti [deg] = %f\n\tomega [deg] = %f\n' ...
    '\tOmega [deg] = %f\n\tM [deg] = %f\n'], floor(OE.satelliteID), ...
    datestr(OE.epoch),OE.a_km, OE.e, OE.i_deg, OE.omega_deg, ...
    OE.Omega_deg, OE.M_deg);
a_m = OE.a_km*1e3;
e = OE.e;
M_deg = OE.M_deg;
%% *Orbital Plane Coordinates*
% In the plane of the orbit:
%  
% $$p= a\left ( \cos E - e \right )$$
%  
% $$q=a\sqrt{1-e^{2}}\sin E$$
% 
% where
%  
% $$p$ - [m] coordinate along axis through center and perigee
%  
% $$q$ - [m] coordinate along axis passing through focus and perpendicular
% to the p-axis
%  
% $$E$ - [rad] eccentric anomaly
%  
% $$e$ - [unitless] eccentricity
%% *Mean Motion*
% The mean motion, n, is the angular speed required for a body to 
% complete one orbit = 2*pi/Period
n = sqrt(mu/a_m^3);  % [rad/s] mean motion
%% *Derivatives Yeild Velocity Components Along Orbit*
%
% $$\dot M = \dot E - e (\cos E) \dot E \\ = n$$
%  
% $$\dot E = \dot M / (1 - e \cos E) \\$$
%  
% $$\dot P = -a (\sin E) \dot E \qquad $$
%  
% $$\dot Q = a (\cos E) \dot E \sqrt{1 - e^2}$$
%
%% *Use Kepler Equation to Find Location Along Orbit*
% $$E-e\sin \left ( E \right )= n\left ( t-t_{p} \right )=M$$
%  
% where
%  
% $$t_{p}$ - [s] time of periapsis
%  
% $$n$ - [rad/s] mean motion (in the TLE it is orbits/day)
%  
% $$M$ - [rad] mean anomaly
%% *Use Newton Method to Solve the Kepler Equation*
% $$E_{i+1} = E_i - f(E_i) / f'(E_i) \\$$
%  
% $$f'(E) = 1 - e \cos E$$
%% *Output Orbit Plane State Vector*
%  p_m - [m] coordinate along axis through center and perigee
%  q_m - [m] coordinate passing through focus and perpendicular to p-axis
%  dpdt_m_per_s = [rad/s] p component velocity
%  dqdt_m_per_s = [rad/s] q component velocity
M_rad = deg2rad(M_deg);
E_rad = M_rad; 
dE = 99999;
eps = 1e-6; % [rad] control precision of Newton's method solution
while (abs(dE) > eps)
    dE = (E_rad - e * sin(E_rad) - M_rad)/(1 - e * cos(E_rad));
    E_rad = E_rad -  dE;
end
p_m = a_m*(cos(E_rad) - e);
q_m = a_m*sqrt(1 - e^2)*sin(E_rad);
dMdt_rad_per_s = n;
dEdt_rad_per_s = dMdt_rad_per_s/(1 - e*cos(E_rad));
dpdt_m_per_s = -a_m*sin(E_rad)*dEdt_rad_per_s;
dqdt_m_per_s = a_m*cos(E_rad)*dEdt_rad_per_s*sqrt(1 - e^2);
%% *Plot the Plane of the Orbit*
ra_m = a_m*(1 + e);  % [m] apogee 
rp_m = a_m*(1 - e);  % [m] perigee
rEarth_m = 6378e3; % [m] Earth radius at equator
Evals = 0:0.01:360.0; % [deg] values of the eccentric anomaly around orbit 
pvals = a_m*(cosd(Evals)-e); % [m] orbit positions
qvals = a_m*sqrt(1 - e^2)*sind(Evals); % [m] orbit positions
figure('color','white');
fill(rEarth_m.*cosd(Evals),rEarth_m.*sind(Evals),[0.75 1.00 0.75]);
hold on;
plot(pvals,qvals,'.b');
plot(p_m,q_m,'.r','MarkerSize',25);
p1_m = p_m+dpdt_m_per_s*300;
q1_m = q_m+dqdt_m_per_s*300;
plot([p_m p1_m],[q_m q1_m],'r');
theta = -atan2d(dqdt_m_per_s,dpdt_m_per_s);
plot(p1_m+[0 -0.1*rEarth_m*cosd(theta+30)],q1_m+[0  0.1*rEarth_m*sind(theta+30)],'r');
plot(p1_m+[0 -0.1*rEarth_m*cosd(theta-30)],q1_m+[0  0.1*rEarth_m*sind(theta-30)],'r');
% axes
plot([-1.25*rp_m 1.25*rp_m],[0 0],'k');
plot([0 0],[-1.25*ra_m 1.25*ra_m],'k');
text(1.25*rp_m, 0.1*rEarth_m, 'p');
plot(1.25*rp_m+[0 -0.1*rEarth_m*cosd(30)],[0 0.1*rEarth_m*sind(30)],'k');
plot(1.25*rp_m+[0 -0.1*rEarth_m*cosd(30)],[0 -0.1*rEarth_m*sind(30)],'k');
text(0.1*rEarth_m, 1.25*ra_m, 'q');
plot([0 -0.1*rEarth_m*sind(30)],1.25*ra_m+[0 -0.1*rEarth_m*cosd(30)],'k');
plot([0 +0.1*rEarth_m*sind(30)],1.25*ra_m+[0 -0.1*rEarth_m*cosd(30)],'k');
axis equal
axis tight
axis off
box on
title({'Initial State Vector of ISS in Orbital Plane'; ...
    ['TLE Epoch: ' datestr(OE.epoch)]});

[Meg Noah]

Hi - that code only has a very old TLE for ISS. It has a switch statement for you edit to add more, or you could replace it altogether to read TLE from a file. TLE become stale after about 30 days due to perturbations.

http://celestrak.com/NORAD/elements/

The code I shared is only for a satellite at epoch. I'll be sharing a code that does orbit decays - there is a poster paper at the December AGU about it. Watch for a new code soon! :-)

[ali hassan]

what if i want to Find out the orbit and location of a PRSS1 satellite revolving around the earth at 11:38 hrs?

what will be the code?

[Meg Noah]

You need to search around for the TLE. If you can't find them on celestrak, they might be somewhere else. For example:

https://www.n2yo.com/satellite/?s=43530

Got this one from that site, and you can just replace the call to getSatelliteTLE with:

TLE = { ...
        'PRSS1'; ... 
'1 43530U 18056B   20341.51316904 +.00000187 +00000-0 [phone]'; ... 
'2 43530 097.9845 056.7485 0001442 079.7157 280.4217 14.76934886130094'; ... 
};

[ali hassan]

btw thanks for your time. is TLE are the kepler elements for any satellite at any specific time or they wont change?

[Meg Noah]

The epoch is the time of the state vector that is described by the TLE. Roughly, they are pretty good for about 30 days or so, but that is when used with propagation models such as SGP4 (Simplified General Perturbations Method). This method interprets the values of the TLE as a set of mean orbital elements given by Brouwer-Lyddane which involve spherical expansion of Earth's gravitation and give factors for accounting for atmosphere drag that are not bad for low density atmosphere and small objects but not great for space stations. The SGP4 factors in gravity of other planets, non-spherical earth gravity, drag, and solar radiation as well as other perturbations. The values of the orbital elements slowly evolve over time. The epoch uncertainty is about 1 km, but for that becomes about 3 km per day from epoch.

Watch my filelist for an upcoming contribution that is an educational lab that will use the RKF4/5 method to propagate a space station (Skylab) with proper atmosphere drag, and some approximations for the perturbations of Right Ascension of the Ascending Node and argument of perigiee. For this lab, I interpret the state vector as Kepler Elements from the TLE, but technically the TLE is using the Brouwer-Lyddane orbital elements.

I'm still editing the poster paper for AGU. But the lab and solutions will be published on matlab.

[ali hassan]

great help from your side.thanks indeed.

whats your degree or program level?

you did this project at undergraduate level??

currently, you are a researcher or master or Phd??

BEST REGARDS

[ali hassan]

and have you also worked on RTL SSDRs and raspberry pis??

do u have any idea about TDOA??

[Meg Noah]

I have completed my graduate work at UMass Lowell (PhD, 2014, Physics). I teach astronomy at Nashua Community College. The above code was written just to answer a question for someone here at matlab central, but I'm making a lab for students that will do orbit propagation using an RKF4/5. Currently on matlab central I've shared a similar propagator for baseball motion showing the 3D deflections due to magnus forces and spin in atmosphere. Generally, I do about TDOA. I've programmed small components before but only played with raspberry pi once a long time ago. Robotics is fun!

[ali hassan]

great work and effort indeed!!really amazing

actually i am undergrad student working on TDOA based localization using RTL SDRs. i am supposed to only implement it in software but i am interested in hardware too. i am using RTL SDRs+raspberry pie as receivers and GPS module to synchronize those receivers, a masterto perform calculations after delays are received. does it sound logical??

i have made the software part and i am attaching it as well.

any suggestions or improvements sir?

BEST REGARDS

[ali hassan]

pakistan has launched a missile on india and israel ? which missile will follow which orbiti.e polar,equatorialetc? i.e polar, equatorial, etc? is it easy for pakistan to launch a missile on india or for india to launch a missile on pakistan? prove mathematically

[Rik]

Your comment was flagged as spam. Although I don't think it is spam, I don't think it is appropriate for this forum either.

[ali hassan]

thanks @Rik. its okay either way. but do agree with you.

BEST REGARDS. i think maybe country names made it spam

[ali hassan]

Modified Rosetta mission to Churyumov-Garasimenko has been designed for scientific study. At burnout of last stage of launch vehicle following information has been extracted from two line element of satellite: Inclination=5.5 , Mean Motion=16.272 rev/day, eccentricity=0. Initial mass of satellite at launch is 3000 kg (including payload of 200kg). Engine installed on the satellite is capable of providing thrust of 250 N and specific impulse of 300 s at mass flow rate of 8.16 kg/s. Satellite follows equatorial hyperbolic departure from earth gravitational field with excess hyperbolic velocity of 2 km/s and arrive at an altitude of 100 km with velocity of 1 km/s at Churyumov-Garasimenko. Final parking orbit of Rosetta satellite is at an altitude of 20 km from Churyumov-Garasimenko. Mission designer has recommended raising of orbit by using following series (1, 1.7276, 7.0379,….) Analyze the mission to establish efficacy, so that mission is accomplished without violation of following constraints.

(a) Engine can’t be used for delta-v greater than 1.54 km/s and for more than 115 seconds for single operation

(b) Total number of burns cannot exceed 6

(c) Total ΔV for all maneuvers in earth gravitational field must not be greater than 4.4 km/s

(d)Empty mass ratio of satellite must not fall below 0.13

Instructions and additional data:-

(a) Clearly draw the figures for all the maneuvers with orbit number (Integer) and points

(b) Radius of Comet Churyumov-Garasimenko is to be taken as 2.92 km

[ali hassan]

can anybody help me out?

[Rik]

https://www.mathworks.com/matlabcentral/answers/6200-tutorial-how-to-ask-a-question-on-answers-and-get-a-fast-answer

https://www.mathworks.com/matlabcentral/answers/8626-how-do-i-get-help-on-homework-questions-on-matlab-answers

[Anirudh Chhabra]

Hi. Can someone point me in the direction as to how to shift the orbit to a different center. For example, here it is around (0,0,0). But I want it in a different location, let's say (1,1,1). How can I do that?

[Douglas Ross]

Hi Meg, thank you for the code, very useful. could you briefly explain what n_deg_per_s is representing? thank you

[Meg Noah]

Yes, Sir. In this context, the parameter 'n' represents the Mean Motion.

In the Two-Line-Element Set (TLE) it is given in orbits per day. Mean Motion is the number of times the satellite orbits the Earth in exactly 24 hours (one solar day). The theoretical limits are between 0 and 17 orbits per day.

The value of Mean Motion in radians per second can be found from:

n - Mean Motion

a - semi-major axis

= 3.98618e14; % [m3/s2] Earth's geocentric gravitational constant

G - universal gravitational constant

Me - mass of the Earth

The Mean Motion is the velocity of the satellite as if it were moving in a circular orbit in uniform circular motion with the same period that it has in its non-circular orbit.

The value of Mean Motion can be expressed in degrees per second by converting radians to degrees.

[Meg Noah]

Here's a satellite orbit with perturbations and space weather impacts (drag).

https://www.mathworks.com/matlabcentral/fileexchange/84310-skylab

[Henrique Chaves]

@Meg Noah sorry for my English and to bother you. But do you know or have you tried something similiar, but for the moon? Precisely lunar orbit satellite's?

Thanks!

[Meg Noah]

Are you looking for a code that predicts the location of the Moon, or are you looking for a code that predicts orbits around the Moon for satellites like Lunar Reconnaissance Orbiter?

For satellites around the Moon, the perturbations due to Earth need to be accounted for. But it is doable with the same basic approach. Space weather impacts are different as well.

[Henrique Chaves]

@Meg Noah

It's for satellites orbiting the moon like the one you mentioned. Yes I know about those perturbations, but for now I'm just looking for something to simulate with kepler's law.

A simple thing, then I will add the perturbations. Don't need to be a real satellite, just give parameters. Example: altitude, shape and period then I have orbit similiar to yours on old comment and output the near point moon and away point.

Do you know something?

[Meg Noah]

OK - just starting with the kepler's approach replace the mass of the earth with the mass of the moon.

G*Me = 3.98618e14; % [m3/s2] Earth's geocentric gravitational constant

G - universal gravitational constant

Me - mass of the Earth

That value will get you the first approximation to the values. You also will need to change the radius of the Earth to the radius of the Moon for calculations of altitude above Moon's surface.

But I don't know how good it will be for many orbits. The perturbations are different for the Moon.

[Henrique Chaves]

Thanks for your answer. Already did what you say.

Do you know somw forum or people so I can talk about it too help me?

[Frank Epstein]

I think there is an error in the code TLE2OrbitalElements.m

Line 72: should be OE.Mdot_orbit_per_day2 = str2double(line1(34:43));

[Frank Epstein]

Can someone expain where the constant 86400 comes from in calculation of the semi-major axis? I guess seconds per day.

146 % [km] semi major axis

147 OE.a_km = ( mu/(str2double(line2(53:63))*2*pi/86400)^2 )^(1/3);

[Meg Noah]

@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]

Meg, thank you for your continued contributions on these topics!

Answer 2

[MANISH SINGH]

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

Answer 3

[Christian Bellinazzi]

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

Answer 4

[Christian Bellinazzi]

I would like also to use getframe for tthe moving marker on the orbit. Anyone could hel me?