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

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);

Published with MATLAB® 7.10