sat_orbit

Contents

• SPG4 orbit propagation in Earth Centered Intertial (ECI) reference
• Convert orbit into ECEF
• Convert orbit into Lat, Lon, H (LLH)
• Check accuracy by compareìing propagated orbit with the reference one

function retv = sat_orbit()

% Orbit Propagator and check for accuracy
% It takes one TLE file and propagate the orbit
% referencing in ECI, ECEF and LLH references.
% It compares the result with the very high accuracy orbit,
% available for the ENVISAT satellite in the example.
% Detailes:
% 1) Read orbit from TLE file, see http://celestrak.com/ or
%   http://www.space-track.org
% 2) Propagate the orbit, by means of NORAD's SGP4
%   in
%   http://www.centerforspace.com/downloads/files/sgp4/AIAA-2006-6753Code.zip
% 3) Compare the propagated orbit with the reference one and plot the error
%
% Author: This interface is by Andrea Monti Guarnieri
%   [e-mail]
%   Orbit propagator is from centerforspace.

SPG4 orbit propagation in Earth Centered Intertial (ECI) reference

tlef='envisat.txt'; % tle file
reforb='doref.txt'; % file with doris reference orbits (high precision)
tsat ='ENVISAT';    % select satellite
dSV=1*60;           % sate Vecotrs interval, s (should be multiple of 1 min)
NCYCLES=10;         % number of orbits
model='SGP4';

% Read header
SAT = mytle(tlef,tsat); % get TLEF useful informations
NORB=floor(SAT.p/dSV);
Nsv = NORB*NCYCLES;  % number of state vectors

% Norad orbit propagation
[sat_rec rV] = norad4(tlef, tsat, floor(Nsv/2), dSV/60);

% Plot orbit in ECI reference
figure(1)
ii=1:Nsv;
plot3(rV.r(1,ii),rV.r(2,ii),rV.r(3,ii),'o-r')
grid;
title('Orbit ECI (inertial) (km)')

Convert orbit into ECEF

[R_ECEF, V_ECEF] = eci2ecef(rV.time, rV.r*1e3, rV.v*1e3); % and put into m, m/s
figure(2)
plot3(R_ECEF(1:end,1)/1e3,R_ECEF(1:end,2)/1e3,R_ECEF(1:end,3)/1e3,'-')
title('Orbit ECEF (km)')
xlabel('X');ylabel('Y');zlabel('Z');
grid

Convert orbit into Lat, Lon, H (LLH)

[lat_rad TH hs]=xyz2llh(R_ECEF(1:end,1), R_ECEF(1:end,2), R_ECEF(1:end,3));
lon= TH*180/pi;
lat= lat_rad*180/pi;

figure(5)
plot(lon,lat,'+');
title([tsat ' Orbit '])
xlabel('Lon')
ylabel('Lat')
grid
axis('equal')

Check accuracy by compareìing propagated orbit with the reference one

Load reference orbit

fp=fopen(reforb,'r');
if(fp == -1)
    error('Cannot open file');
end
jj=1;
while ~feof(fp)
    line=fgetl(fp);
    E.t(jj) = datenum(line(1:26));
    A=sscanf(line,'%*s %*s %g %g %g %g %g %g %g %g %g');
    E.r(jj,1:3) = A(3:5);
    E.v(jj,1:3) = A(6:8);
    jj=jj+1;
end
fclose(fp);

% find initial time
difftime=(E.t-rV.time(1));
indtp = find(difftime>=0);
indt0 = indtp(1);

% Resample on the same axis and check error
% limit the state vectors to Npd starting from nshift
% (sv spacing is 1 minute)
nshift=500;
Npd=15;
ii=nshift+indt0:nshift+indt0+Npd-1;
r_dor = E.r(ii,:);
t_dor = E.t(ii);

% use a larger set of vector for envisat, not to extrapolate
Npe=Npd+4;
ii=1+nshift:Npe+nshift;
r_spg = R_ECEF(ii,:);
t_spg = rV.time(ii)';

% Do orbit interpolation to compute error
% normalize and scale time axis
t0 = t_dor((Npd-1)/2);
ts = t_dor(end)-t0;
td = (t_dor-t0)/ts;
te = (t_spg-t0)/ts;

No=7; % oder for polinomial interpolation
% do polynomial fitting
px=zeros(3,No+1);
for jj=1:3;
    [px(jj,:)] = polyfit(te,r_spg(:,jj),No);
end

% interpolate on reference orbit time
r1_spg = zeros(Npd,3);
for jj=1:3;
    r1_spg(:,jj) = polyval(px(jj,:),td);
end

% compare the two sets
% compute the error for each sample
errd=zeros(1,Npd);
for jj=1:Npd
    errd(jj)=norm(r1_spg(jj,:)-r_dor(jj,:),2);
end
figure(7)
plot(t_dor-t0,errd,'r');
xlabel('Time (frac days)')
title('Error in the propagated orbit')
errm = min(errd)        % error in m
retv=0;

errm =

  367.0441

end

ans =

     0

Published with MATLAB® 7.7