The ECI Satellite Position Calculation in Matlab

 

Given the satellites orbital parameters, the computation of its orbit is quite easy. The Matlab code that follows shows how to perform this calculation.

Take into account that results are in the ECI reference system. To make the computation in the ECEF reference system, some more information are needed (Om = Om0 - OmEdot*(t - toe), where Om0 is the right ascension of the orbit in ECI, OmEdot is the Earth angular velocity [rad/s], t is the actual time and toe is the time passed from the last time that the ECEF reference system was identical to the ECI reference system.

 

The ouput of the function, called with the following parameters:

P = SatPos(20e6, 0.5, 0.3, 0.4, 0, 0, [0:20:28200]);

 

is depicted in the following picture

 

t = [0:20:28200] is the time vector for which the orbit is computed, i.e., every 20 sec, from 0 to 28200 seconds.

The result of the function, P, is the vector of ECI positions corresponding to the t input vector

 

% SATPOS
% 
% Function that computes the position of a satellite in the orbit
% This function doesn't take into account the orbit perturbations.
% 
% Input parameters:
% a: Semi major axis (meters)
% e: Eccentricity
% i: Inclination angle of the orbit [rad]
% Om: Right ascension [rad]
% om: Argument of perigee [rad]
% M0: Mean anomaly at time t = 0 [rad]
% t: time vector for which the orbit should be computed
% 
% Output parameters:
% P = [x,y,z]: ECI position of the satellite
% 
% (C) 2005, Leonardo Daga
function [P] = SatPos(a, e, i, Om, om, M0, t)

% Constants definition
muT=3.986005e14 ; % Earth Universal Gravitational Constant [m^3/s^2]

% Mean orbit angular speed
N = sqrt(muT/(a*a*a));

% Mean anomaly
M = M0 + N*t; 

% Eccentric anomaly evaluation using Kepler equation
E = kepler(M, e, 1e-10); 

% Calculation of the real anomaly using elliptical coordinates
es = ((sqrt(1-e*e))*sin(E))./(1-e*cos(E));
ec = (cos(E)-e)./(1-e*cos(E));
v = atan2(es,ec); 

% Actual distance of the satellite from the center of the earth
r = a*(1-e*cos(E)); 

% Rotation of the perigee
u = v + om;

% X and Y in the orbital plane
xp = r .* cos(u); 
yp = r .* sin(u); 

% Rotation of the orbital plane due to the inclination and right ascension
x = xp*cos(Om) - yp*cos(i)*sin(Om); 
y = xp*sin(Om) + yp*cos(i)*cos(Om); 
z = yp*sin(i); 

P = [x', y', z'];

% Disable this function if you don't want the drawing of the orbit
figure(1),plot3(P(:,1),P(:,2),P(:,3)),hold on, disEarth, hold off


% KEPLER Iterative calculation of the eccentric anomaly from
% the mean anomaly. The solved equation is the kepler equation for the
% elliptic orbits: u - e*sin(u) = M.
% 
% Input
% M: Mean Anomaly
% e: Eccentricity
% epsl: desired precision ( 1e-16 default)
%
% Output
% u: Resulting eccentric anomaly
function [u]=kepler(M,e,epsl);

if nargin<3
epsl = 1e-12;
end

M = rem(M,2*pi);

% Initial approximation
u = M;

% Start condition (i.e. u0 != u)
u0 = u+1;
while (sum(abs(u-u0))>epsl)
u0 = u;
u = u - (u-e.*sin(u)-M)./(1-e.*cos(u));
end
end 

% DISEARTH draws the earth globe
function disEarth()

rad = 6.38e6;

[xs,ys,zs]=sphere(15);

xs = xs*rad;
ys = ys*rad;
zs = zs*rad;

obj = surfl(xs,ys,zs)

grid on;
 

Leonardo Daga's Warehouseâ, http://leonardodaga.insyde.it
Send any Comments to: leonardo.daga@gmail.com