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, 1e10);
% Calculation of the real anomaly using elliptical coordinates
es = ((sqrt(1e*e))*sin(E))./(1e*cos(E));
ec = (cos(E)e)./(1e*cos(E));
v = atan2(es,ec);
% Actual distance of the satellite from the center of the earth
r = a*(1e*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 ( 1e16 default)
%
% Output
% u: Resulting eccentric anomaly
function [u]=kepler(M,e,epsl);
if nargin<3
epsl = 1e12;
end
M = rem(M,2*pi);
% Initial approximation
u = M;
% Start condition (i.e. u0 != u)
u0 = u+1;
while (sum(abs(uu0))>epsl)
u0 = u;
u = u  (ue.*sin(u)M)./(1e.*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;

