Small Satellites

Home » 2012 » December

Monthly Archives: December 2012

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

 

Advertisements

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

CubeSats currently in orbit (Number of CubeSats 44, December 08, 2012)

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

Tension in the yo-yo cables
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

 

%d bloggers like this: