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);
Euler rotation example, Rotation matrix, Quaternion, Euler Axis and Principal Angle
A classical Euler rotation involves first a rotation about e3 axis, then one about the e1 axis and finally a rotation about the e3 axis.
% Compute the elements of rotation matrix R
clc;
psi = pi/4;
theta = pi/4;
fi = pi/6;
R3_psi = [cos(psi) sin(psi) 0;
-sin(psi) cos(psi) 0;
0 0 1]
R1_theta = [1 0 0;
0 cos(theta) sin(theta);
0 -sin(theta) cos(theta)]
R3_fi = [cos(fi) sin(fi) 0;
-sin(fi) cos(fi) 0;
0 0 1]
R3_psi =
0.7071 0.7071 0
-0.7071 0.7071 0
0 0 1.0000
R1_theta =
1.0000 0 0
0 0.7071 0.7071
0 -0.7071 0.7071
R3_fi =
0.8660 0.5000 0
-0.5000 0.8660 0
0 0 1.0000
Multiplying each rotation matrix
R = R3_fi*R1_theta*R3_psi
R =
0.3624 0.8624 0.3536
-0.7866 0.0795 0.6124
0.5000 -0.5000 0.7071
From Rotation Matrix Compute Quaternion
q4 = 0.5*(1 + R(1,1)+ R(2,2) + R(3,3))^0.5; q1 = (R(2,3) - R(3,2))/(4*q4); q2 = (R(3,1) - R(1,3))/(4*q4); q3 = (R(1,2) - R(2,1))/(4*q4); q = [q1 q2 q3 q4] norm_q = norm(q) % Checking that the norm of q = 1 qv = [q1 q2 q3];
q =
0.3794 0.0500 0.5624 0.7330
norm_q =
1.0000
Compute Euler axis and its components along the axis E1,E2,E3 unit vector. Compute the Euler principle angle.
Euler_angle = 2*acos(q4)*180/pi % [deg] Euler_axis = qv/sind(Euler_angle/2) % [E1,E2,E3] norm_E = norm(Euler_axis)
Euler_angle =
85.7293
Euler_axis =
0.5577 0.0734 0.8268
norm_E =
1.0000
CubeSats currently in orbit
The program extracts CubeSats orbital elements from TLE sets.
TLE Data Source http://www.celestrak.com/ December 08, 2012
clc; clear all; mu = 398600; % Standard gravitational parameter for the earth fname = 'cubesat.txt'; % TLE file name % Open the TLE file and read TLE elements fprintf('a[km] e inc[deg] RAAN[deg] w[deg] M[deg] T[min] h_p[km] h_a[km]\n\n ') fid = fopen(fname, '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]); epoch = L2(1,4)*24*3600; % Epoch Date and Julian Date Fraction Db = L2(1,5); % Ballistic Coefficient 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); % names(i,:) = [L1]; i = i+1; end fclose(fid); scatter(a,inc,'b*'); grid on; xlabel('Semi-major axis [km]'); ylabel('Inclination [deg]'); title('CubeSats currently in orbit (Number of CubeSats 44, December 08, 2012)');
a[km] e inc[deg] RAAN[deg] w[deg] M[deg] T[min] h_p[km] h_a[km] DTUSAT 7198.85 0.0010129 98.6942 348.7037 60.0947 300.1226 101.31 820.56 835.14 CUTE-1 (CO-55) 7199.68 0.0010503 98.7014 347.9287 82.8078 277.4288 101.33 821.12 836.24 QUAKESAT 7200.51 0.0009603 98.7140 347.8920 102.1514 258.0737 101.35 822.60 836.43 AAU CUBESAT 7198.85 0.0010130 98.6942 349.1780 59.1696 301.0464 101.31 820.56 835.15 CANX-1 7198.96 0.0010048 98.6912 348.6908 60.2726 299.9452 101.31 820.73 835.20 CUBESAT XI-IV (CO-57) 7200.57 0.0010385 98.7106 347.8636 86.6453 273.5908 101.35 822.09 837.04 UWE-1 7068.55 0.0016644 97.8893 199.3059 330.1747 29.8495 98.57 685.78 709.31 CUBESAT XI-V (CO-58) 7068.46 0.0017244 97.8889 201.0839 328.4579 31.5593 98.57 685.27 709.65 NCUBE-2 7067.65 0.0016749 97.8885 200.8900 329.0837 30.9341 98.55 684.81 708.48 CSTB1 7081.95 0.0082588 97.8642 349.7817 272.4294 86.7438 98.85 652.46 769.44 MAST 7089.54 0.0092488 97.8744 343.5208 290.6610 68.4682 99.01 652.97 784.11 LIBERTAD-1 7093.55 0.0100613 97.8862 340.1948 303.8714 55.3008 99.10 651.18 793.92 POLYSAT CP3 7092.78 0.0099844 97.8854 340.6063 302.8126 56.3494 99.08 650.97 792.60 CAPE1 7093.28 0.0100140 97.8855 340.2673 303.5659 55.5990 99.09 651.25 793.31 POLYSAT CP4 7083.44 0.0083192 97.8614 348.9130 272.9657 86.2018 98.88 653.51 771.37 NTS (CANX-6) 6998.86 0.0014784 97.7781 37.8172 211.4274 148.6085 97.12 617.51 638.21 CUTE-1.7+APD II (CO-65)6996.61 0.0013710 97.7788 39.1393 212.2206 147.8143 97.07 616.02 635.21 COMPASS-1 6992.51 0.0014294 97.7785 39.5566 211.0446 148.9960 96.99 611.51 631.50 AAUSAT-II 6991.27 0.0013929 97.7791 40.9611 207.9487 152.0989 96.96 610.54 630.01 DELFI-C3 (DO-64) 6984.92 0.0013654 97.7929 43.2390 204.4542 155.6015 96.83 604.38 623.45 CANX-2 6996.94 0.0014274 97.7753 38.5852 213.9057 146.1242 97.08 615.95 635.92 SEEDS II (CO-66) 6994.77 0.0014702 97.7794 39.1404 212.0131 148.0200 97.03 613.48 634.05 SWISSCUBE 7092.97 0.0006783 98.3468 85.0891 301.0710 58.9603 99.08 717.16 726.78 BEESAT 7091.53 0.0007596 98.3461 84.9669 316.3983 43.6207 99.05 715.14 725.91 UWE-2 7091.94 0.0006165 98.3412 85.1781 317.1499 42.9051 99.06 716.56 725.31 ITUPSAT 1 7093.66 0.0007926 98.3553 85.1249 299.8155 60.2062 99.10 717.03 728.28 TISAT 1 6999.16 0.0014077 98.0502 60.7150 332.8485 27.1989 97.12 618.31 638.01 DICE-F 7003.61 0.0244060 101.7129 102.4166 215.4701 143.0058 97.22 461.68 803.54 DICE-Y 7001.82 0.0242032 101.7129 103.1190 214.0671 144.4790 97.18 461.35 800.28 RAX-2 7002.19 0.0242055 101.7075 102.3303 215.0375 143.4693 97.19 461.70 800.68 AUBIESAT-1 7000.58 0.0243161 101.7078 101.8965 215.8001 142.6690 97.15 459.36 799.81 M-CUBED & EXP-1 PRIME7001.26 0.0243566 101.7076 103.1161 213.4747 145.0857 97.17 459.74 800.79 E-ST@R 7174.06 0.0683773 69.4844 274.8017 26.2580 337.1524 100.79 312.51 1293.60 2012-006D 7173.48 0.0681813 69.4744 274.6267 26.5331 336.9035 100.78 313.38 1291.57 MASAT 1 7174.00 0.0683197 69.4772 274.9223 26.5308 336.9128 100.79 312.88 1293.13 XATCOBEO 7158.68 0.0664444 69.4820 272.6345 24.9101 338.2511 100.46 312.03 1263.34 PW-SAT 7165.74 0.0672363 69.4793 273.7501 25.6794 337.6108 100.61 312.95 1276.54 2012-006H 7177.12 0.0689452 69.4814 275.4853 26.6354 336.8405 100.85 311.29 1300.94 2012-006J 7178.71 0.0691250 69.4831 275.7959 26.7936 336.7126 100.89 311.48 1303.93 RAIKO 6777.72 0.0017840 51.6592 339.7729 11.4627 348.6772 92.55 394.63 418.82 FITSAT-1 (NIWAKA) 6775.20 0.0015861 51.6506 340.6207 23.6131 336.5585 92.50 393.45 414.94 TECHEDSAT 6769.70 0.0015823 51.6483 340.0468 26.7198 333.4605 92.39 387.99 409.41 F-1 6770.08 0.0016123 51.6502 335.2502 29.5738 330.6164 92.40 388.17 410.00 WE-WISH 6758.57 0.0018934 51.6466 334.4419 26.8172 333.3788 92.16 374.77 400.36
Radial Release
Example 10.8, Orbital Mechanics for Engineering Students, 2nd Edition.
A satellite is to be completely despin using a two-mass yo-yo device with tangential release. Assume the spin axis of moment of inertia of the satellite is C = 200 kg · m^2 and the initial spin rate w0 = 5 [rad/s]. The total yo-yo mass is 4 kg, and the radius of the spacecraft is 1 meter.
clc; C = 200; % [kg*m^2] w0 = 5; % [rad/s] m = 4; % [kg] R = 1; % [m] K = 1 + C/(m*R^2); % Nondimensional factor K
The required cord length l
l = R*(K)^0.5 %[m] % Time t to despin t = (K)^0.5/w0 %[s] % Tension in the yo-yo cables t = 0:0.001:5; N = 2*K*C*w0^3*t./(R*(K + w0^2*t.^2).^2); plot(t,N); xlabel('time [s]'); ylabel('Tension N [N]'); grid on; [maxN,i] = max(N);
l =
7.1414
t =
1.4283
The maximum tension in the yo-yo cables [N];
maxN % which ocurs at t = i*0.001 % [s]
maxN =
454.7542
t =
0.8260
The speed of the masses at release;
v = R*w0*(K)^0.5 % [m/s]
v = 35.7071
The angle rotated by the satellite during the despin;
theta = (pi/2 - 1)*(K)^0.5 % [rad]
theta =
4.0763
Cord length required for complete despin
l = R*(K^0.5 - 1) % [m]
l =
6.1414
Published with MATLAB® 7.10
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 