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');
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
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
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);
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
Small Satellites 