Small Satellites

Home » Posts tagged 'TLE'

Tag Archives: TLE

J2 Propagator

In this example we will implement the J2 Perturbation (first-order) propagator which accounts for secular variations in the orbit elements due to Earth oblateness. We do not consider atmospheric drag or solar or lunar gravitational forces.It uses two-line element (TLE) sets as an input. Additionally as an example we wil use the implemented propagator to plot orbits of RADUGA satellites.

% TLE Data Source. http://www.celestrak.com/  January 20, 2012
clc;
clear all;
Earth3DPlot(1);

mu    = 398600;          % Earth’s gravitational parameter [km^3/s^2]
R_earth = 6378;          % Earth radius [km]
J2 = 0.0010836;
we = 360*(1 + 1/365.25)/(3600*24);      % Earth's rotation [deg/s]
fname = 'Raduga.txt';    % TLE file name
% Open the TLE file and read TLE elements
fid = fopen(fname, 'rb');
while ~feof(fid)
% Reading TLE File
L1 = fscanf(fid,'%23c%*s',1);
L2 = fscanf(fid,'%d%6d%*c%5d%*3c%*2f%f%f%5d%*c%*d%5d%*c%*d%d%5d',[1,9]);
L3 = fscanf(fid,'%d%6d%f%f%f%f%f%f%f%f',[1,9]);
epoch       = L2(1,4)*24*3600;        % Epoch Date and Julian Date Fraction
Db          = L2(1,5);                % Ballistic Coefficient
i           = L3(1,3);                % Inclination [deg]
RAAN        = L3(1,4);        % Right Ascension of the Ascending Node [deg]
e           = L3(1,5)/1e7;            % Eccentricity
omega       = L3(1,6);                % Argument of periapsis [deg]
M           = L3(1,7);                % Mean anomaly [deg]
n           = L3(1,8);                % Mean motion [Revs/day]
% Orbital parametres
a = (mu/(n*2*pi/(24*3600))^2)^(1/3);         % Semi-major axis [km]
T = 2*pi*sqrt(a^3/mu)/60;                    % Period in [min]
rp = a*(1-e);
h = (mu*rp*(1 + e))^0.5;                     % Angular momentum
E = keplerEq(M,e,eps);
theta =  acos((cos(E) -e)/(1 - e*cos(E)))*180/pi;    % [deg] True anomaly

% J2 Propagator
dRAAN = -(1.5*mu^0.5*J2*R_earth^2/((1-e^2)*a^3.5))*cosd(i)*180/pi;
domega = dRAAN*(2.5*sind(i)^2 - 2)/cosd(i);
eps = 1E-9;
dt  = 5;        % time step [sec]
ti  = 0;
j = 1;
while(ti <= T);
    E = 2*atan(tand(theta/2)*((1-e)/(1+e))^0.5);
    M = E  - e*sin(E);
    t0 = M/(2*pi)*T;
    t = t0 + dt;
    M = 2*pi*t/T;
    E = keplerEq(M,e,eps);
    theta = 2*atan(tan(E/2)*((1+e)/(1-e))^0.5)*180/pi;
    RAAN  = RAAN  +  dRAAN*dt ;
    omega = omega + domega*dt;
    [R V] = Orbital2State( h, i, RAAN, e,omega,theta);
    % Considering Earth's rotation
    fi_earth = we*ti;
    Rot = [cosd(fi_earth), sind(fi_earth),0;...
        -sind(fi_earth),cosd(fi_earth),0;0,0,1];
    R = Rot*R;
    ti = ti+dt;
    Rv(j,:) = R';
    j = j+1;
end
scatter3(Rv(:,1),Rv(:,2),Rv(:,3),'.');
end
fclose(fid);
hold off;

% Ground Track
% Open the TLE file and read TLE elements
EarthTopographicMap( 2,820,420);
fid = fopen(fname, 'rb');
while ~feof(fid)
[ h, i, RAAN, e,omega,theta ] = TLE2OE( fid );
[Rv, alfa,delta ] = J2PropagR( h, i, RAAN, e,omega,theta) ;
scatter(alfa,delta,'.');
end
fclose(fid);
figure(2);
legend('RADUGA 32', 'RADUGA-1 4', 'RADUGA-1 5', 'RADUGA-1 7', 'RADUGA-1M 1');
 text(270,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
	'HorizontalAlignment','left','FontSize',14 );
title('RADUGA satellites ground track');

J2Propagator_02Raduga Orbit

Input TLE File: Source  www.celestrak.com

Raduga.txt
RADUGA 32               
1 23448U 94087A   13017.58278277  .00000000  00000-0  10000-3 0  5863
2 23448  12.7840  37.1920 0005679 124.1377 235.9487  1.00272544 83823
RADUGA-1 4              
1 25642U 99010A   13017.71286665  .00000130  00000-0  10000-3 0  6727
2 25642  11.7705  50.8546 0001080 106.9153 253.1213  1.00296963 50866
RADUGA-1 5              
1 26477U 00049A   13019.59753669 -.00000103  00000-0  10000-3 0  5404
2 26477   9.3910  51.0039 0002612 354.7657   5.2534  1.00203902 45393
RADUGA-1 7              
1 28194U 04010A   13018.70907761  .00000086  00000-0  10000-3 0  6657
2 28194   6.8757  64.2408 0001484  24.4802 335.5290  1.00365778 32290
RADUGA-1M 1             
1 32373U 07058A   13018.90229375 -.00000043  00000-0  00000+0 0  4653
2 32373   0.0263  70.5780 0002023 210.8419 231.8461  1.00272457 19621
Advertisements

Space Debris

The program plots track-able space debris distribution in the orbit from three  major accidents.

'Cosmos 2251/Iridium 33 ', 'BREEZE-M 2012-044C ', 'FENGYUN 1C'

TLE Data Source. http://www.celestrak.com/ December 08, 2012

clc;
clear all;
close all;
mu = 398600;             % Standard gravitational parameter for the earth
% Input TLE files
% Cosmos Debris TLE1,  Iridium Debris TLE2, BREEZE-M Debris TLE3'
% FENGYUN Debris TLE4
names = {'Cosmos 2251 ' 'Iridium 33 ' 'BREEZE-M 2012-044C ' 'FENGYUN 1C'};
tPieces = 0;
for ind = 1:4
filename = ['TLE' num2str(ind),'.txt'];
fid      = fopen(filename,'rb');
i = 1;
while ~feof(fid)
L1 = fscanf(fid,'%23c%*s',1);
L2 = fscanf(fid,'%d%6d%*c%5d%*3c%*2f%f%f%5d%*c%*d%5d%*c%*d%d%5d',[1,9]);
L3 = fscanf(fid,'%d%6d%f%f%f%f%f%f%f%f',[1,9]);

inc(i)= L3(1,3);                % Inclination [deg]
RAAN(i) = L3(1,4);              % Right Ascension of the Ascending Node [deg]
e(i)  = L3(1,5)/1e7;            % Eccentricity
w(i)  = L3(1,6);                % Argument of periapsis [deg]
M(i)  = L3(1,7);                % Mean anomaly [deg]
n     = L3(1,8);                % Mean motion [Revs per day]

% Orbital parametres
a(i) = (mu/(n*2*pi/(24*3600))^2)^(1/3);       % Semi-major axis [km]
T(i) = 2*pi*sqrt(a(i)^3/mu)/60;               % Period in [min]
Re = 6371;
h_p(i) = (1 - e(i))*a(i) - Re;                % Perigee Altitude [km]
h_a(i) = (1 + e(i))*a(i) - Re;                % Apogee Altitude [km]

% Orbital Elements
OE = [a(i) e(i) inc(i) RAAN(i) w(i) M(i) T(i) h_p(i) h_a(i)];
%fprintf(L1);
%fprintf('%4.2f  %7.7f   %4.4f  %4.4f   %4.4f  %4.4f  %4.2f  %4.2f  %4.2f \n', OE);
i = i+1;
end
nPieces = i - 1;            % Number of Debris Pieces
fclose(fid);
figure(ind);
scatter(a,inc,'b.');
grid on;
xlabel('Semi-major axis [km]');
ylabel('Inclination [deg]');
title_name = [names(ind),'Debris Cloud'];
legend(['Number of Pieces ', num2str(nPieces),'/Dec 08, 2012'])
title(title_name);

figure(ind+5);
scatter(h_p,h_a,'b.');
grid on;
xlabel('Perigee Altitude [km]');
ylabel('Apogee Altitude [km]');
title_name = [names(ind),'Debris Cloud'];
legend(['Number of Pieces ', num2str(nPieces),'/Dec 08, 2012'])
title(title_name);
av{ind}   = a;
incv{ind} = inc;
tPieces = tPieces + nPieces;
end

Spacedebri_01 Spacedebri_02 Spacedebri_03 Spacedebri_04 Spacedebri_05 Spacedebri_06 Spacedebri_07 Spacedebri_08

Combined Debris Plot

figure(5);
hold on;
grid on;
x = av{:,1};
y = incv{:,1};
scatter(x',y','k*');
x = av{:,2};
y = incv{:,2};
scatter(x',y','b*');
x = av{:,3};
y = incv{:,3};
scatter(x',y','g.');
x = av{:,4};
y = incv{:,4};
scatter(x',y','r.');
xlabel('Semi-major axis [km]');
ylabel('Inclination [deg]');
title_name = ['Combined Debris Cloud (Total number of pieces ',...
    num2str(tPieces),'/Dec 08, 2012 )'];
legend('Cosmos 2251 ', 'Iridium 33 ', 'BREEZE-M 2012-044C ', 'FENGYUN 1C');
title(title_name);

Spacedebri_09

 

Two-line element set (TLE)

 

clear all;
mu = 398600; %  Standard gravitational parameter for the earth
% TLE file name
fname = 'tle.txt';
% Open the TLE file and read TLE elements
fid = fopen(fname, 'rb');
L1c = fscanf(fid,'%24c%',1);
L2c = fscanf(fid,'%71c%',1);
L3c = fscanf(fid,'%71c%',1);
fprintf(L1c);
fprintf(L2c);
fprintf([L3c,'\n']);
fclose(fid);
% Open the TLE file and read TLE elements
fid = fopen(fname, 'rb');
L1 = fscanf(fid,'%24c%*s',1);
L2 = fscanf(fid,'%d%6d%*c%5d%*3c%*2f%f%f%5d%*c%*d%5d%*c%*d%d%5d',[1,9]);
L3 = fscanf(fid,'%d%6d%f%f%f%f%f%f%f',[1,8]);
fclose(fid);

epoch = L2(1,4)*24*3600;        % Epoch Date and Julian Date Fraction
Db    = L2(1,5);                % Ballistic Coefficient
inc   = L3(1,3);                % Inclination [deg]
RAAN  = L3(1,4);                % Right Ascension of the Ascending Node [deg]
e     = L3(1,5)/1e7;            % Eccentricity
w     = L3(1,6);                % Argument of periapsis [deg]
M     = L3(1,7);                % Mean anomaly [deg]
n     = L3(1,8);                % Mean motion [Revs per day]

% Orbital elements

a = (mu/(n*2*pi/(24*3600))^2)^(1/3);     % Semi-major axis [km]

% Calculate the eccentric anomaly using Mean anomaly
err = 1e-10;            %Calculation Error
E0 = M; t =1;
itt = 0;
while(t)
       E =  M + e*sind(E0);
      if ( abs(E - E0) < err)
          t = 0;
      end
      E0 = E;
      itt = itt+1;
end

% Six orbital elements
OE = [a e inc RAAN w E];
fprintf('\n a [km]   e      inc [deg]  RAAN [deg]  w[deg]    E [deg] \n ')
fprintf('%4.2f  %4.4f   %4.4f       %4.4f     %4.4f    %4.4f', OE);
ISS (ZARYA)             
1 25544U 98067A   12341.40243075  .00007074  00000-0  12274-3 0  7812
2 25544  51.6485 346.0291 0016039  10.8890 135.6421 15.51944191804768

 a [km]   e      inc [deg]  RAAN [deg]  w[deg]    E [deg] 
 6789.18  0.0016   51.6485       346.0291     10.8890    135.6432

Published with MATLAB® 7.10

 

%d bloggers like this: