Small Satellites

Home » Posts tagged 'Orbital Mechanics'

Tag Archives: Orbital Mechanics

AGI STK 10 MATLAB INTERFACE: Satellite Ground Track

close all; clear all; clc

 AGI STK, is  as a software package from Analytical Graphics, Inc.(AGI) that
 allows to perform complex analyses of ground, sea, air, and space
 missions. More information can be found in AGI website.  https://www.agi.com
 To create this code we used educational code samples from % https://www.agi.com/resources/

% Establish the connection AGI STK 10
try
    % Grab an existing instance of STK 10
    uiapp = actxGetRunningServer('STK10.application');
catch
    % STK is not running, launch new instance
    % Launch a new instance of STK10 and grab it
    uiapp = actxserver('STK10.application');
end
%get the root from the personality
%it has two... get the second, its the newer STK Object Model Interface as
%documented in the STK Help
root = uiapp.Personality2;

% set visible to true (show STK GUI)
uiapp.visible = 1;
% From the STK Object Root you can command every aspect of the STK GUI
% close current scenario or open new one
try
    root.CloseScenario();
    root.NewScenario('SatelliteGroundTrack');
catch
    root.NewScenario('SatelliteGroundTrack');
end
% Set units to utcg before setting scenario time period and animation period
root.UnitPreferences.Item('DateFormat').SetCurrentUnit('UTCG');
% %set units to epoch seconds because this works the easiest in matlab
% root.UnitPreferences.Item('DateFormat').SetCurrentUnit('EPSEC');

% Set scenario time period and animation period
root.CurrentScenario.SetTimePeriod('25 May 2013 12:00:00.000', '26 May 2013 12:00:00.000');
root.CurrentScenario.Epoch = '25 May 2013 12:00:00.000';

% Create satellite
satObj = root.CurrentScenario.Children.New('eSatellite', 'SmallSats1');

% Propagate satellite
satObj.Propagator.InitialState.Representation.AssignClassical(...
    'eCoordinateSystemJ2000', 6750, 0.1, 53.4, 0, 0, 0);
% CoordinateSystem, Semimajor Axis, Eccentricity, Inclination,
% Arg. of Perigee, RAAN, Mean Anomaly

satObj.Propagator.StartTime = '25 May 2013 12:00:00.000';
satObj.Propagator.StopTime  = '25 May 2013 15:00:00.000';
satObj.Propagator.Propagate;

% Get Latitude, Longitude  for the satellite over the course of the mission.
LLAState = satObj.DataProviders.Item('LLA State').Group.Item('Fixed');
Elems = {'Time';'Lat';'Lon'};
satStartTime = satObj.Propagator.EphemerisInterval.FindStartTime;
satStopTime = satObj.Propagator.EphemerisInterval.FindStopTime;
Results = LLAState.ExecElements(satStartTime, satStopTime, 10, Elems);
time = cell2mat(Results.DataSets.GetDataSetByName('Time').GetValues);
Lat  = cell2mat(Results.DataSets.GetDataSetByName('Lat').GetValues);
Long = cell2mat(Results.DataSets.GetDataSetByName('Lon').GetValues);

Plot

figure(1);
hold on;
axis([0 360 -90 90]);
load('topo.mat','topo','topomap1');
contour(0:359,-89:90,topo,[0 0],'b')
axis equal
box on
set(gca,'XLim',[-180 180],'YLim',[-90 90], ...
    'XTick',[-180 -120 -60 0 60 120 180], ...
    'Ytick',[-90 -60 -30 0 30 60 90]);
image([-180 180],[-90 90],topo,'CDataMapping', 'scaled');
colormap(topomap1);
ylabel('Latitude [deg]');
xlabel('Longitude [deg]');
title('Satellite Ground Track');
scatter(Long,Lat,'.r');

Simplesat_01Ground track generated by AGI STK10  

Advertisements

Relative Motion of Satellites, Numerical Simulation

This example shows how to use the state vectors of spacecraft A and B to find the position, velocity and acceleration of B relative to A in the LVLH frame attached to A. We numerically solve fundamental equation of relative two-body motion to obtain the trajectory of spacecraft B relative to spacecraft A and the distance between two satellites. For more details about the theory and algorithm look at Chapter 7 of H. D. Curtis, Orbital Mechanics for Engineering Students,Second Edition, Elsevier 2010

clc;
clear all;
% Initial State vectors of Satellite A and B
RA0 = [-266.77,  3865.8, 5426.2];     % [km]
VA0 = [-6.4836, -3.6198, 2.4156];     % [km/s]
RB0 = [-5890.7, -2979.8, 1792.2];     % [km]
VB0 = [0.93583, -5.2403, -5.5009];    % [km/s]
mu  = 398600;            % Earth’s gravitational parameter [km^3/s^2]
t   = 0;                 % initial time
dt  = 15;                % Simulation time step [s]
dT  = 4*24*3600;         % Simulation interval  [s]
% Using fourth-order Runge–Kutta method to solve fundamental equation
% of relative two-body motion
F_r = @(R) -mu/(norm(R)^3)*R;
VA = VA0; RA  = RA0;
VB = VB0; RB  = RB0;
ind = 1;
while (t <= dT)
    % Relative position
    hA = cross(RA, VA);     % Angular momentum of A
    % Unit vectors i, j,k of the co-moving frame
    i = RA/norm(RA);  k = hA/norm(hA); j = cross(k,i);
    % Transformation matrix Qxx:
    QXx   = [i; j; k];
    Om    = hA/norm(RA)^2;
    Om_dt = -2*VA*RA'/norm(RA)^2.*Om;
    % Accelerations of A and B,inertial frame
    aA = -mu*RA/norm(RA)^3;
    aB = -mu*RB/norm(RB)^3;
    % Relative position,inertial frame
    Rr = RB - RA;
    % Relative position,LVLH frame attached to A
    R_BA(ind,:) = QXx*Rr';

    % A Satellite
    k_1 = dt*F_r(RA);  k_2 = dt*F_r(RA+0.5*k_1);
    k_3 = dt*F_r(RA+0.5*k_2);  k_4 = dt*F_r(RA+k_3);
    VA  = VA + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    RA  = RA + VA*dt;
    % B Satellite
    k_1 = dt*F_r(RB);  k_2 = dt*F_r(RB+0.5*k_1);
    k_3 = dt*F_r(RB+0.5*k_2);     k_4 = dt*F_r(RB+k_3);
    VB = VB + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    RB  = RB + VB*dt;

    R_A(ind,:)  = RA;
    R_B(ind,:)  = RB;
    time(ind)   = t;
    t   = t+dt;
    ind = ind+1;
end
r_BA = (R_BA(:,1).^2+R_BA(:,2).^2+R_BA(:,3).^2).^0.5;
close all;
figure(1);
hold on;
plot3(R_A(:,1),R_A(:,2),R_A(:,3),'r');
plot3(R_B(:,1),R_B(:,2),R_B(:,3),'y');
title('Satellites orbits around earth');
legend('Satellite A','Satellites B');
xlabel('km');ylabel('km');
% Plotting Earth
load('topo.mat','topo','topomap1');
colormap(topomap1);
% Create the surface.
radius_earth=6378;
[x,y,z] = sphere(50);
x =radius_earth*x;
y =radius_earth*y;
z =radius_earth*z;
props.AmbientStrength = 0.1;
props.DiffuseStrength = 1;
props.SpecularColorReflectance = .5;
props.SpecularExponent = 20;
props.SpecularStrength = 1;
props.FaceColor= 'texture';
props.EdgeColor = 'none';
props.FaceLighting = 'phong';
props.Cdata = topo;
surface(x,y,z,props);
hold off;

figure(2);
plot3(R_BA(:,1),R_BA(:,2),R_BA(:,3),'k');
title('The trajectory of spacecraft B relative to spacecraft A');
xlabel('km');ylabel('km');zlabel('km');

figure(3);
plot(time/3600,r_BA);
title('Distance between two satellites');
xlabel('hour');ylabel('km')

min_r = min(r_BA);
max_r = max(r_BA);
fprintf('Max distance between two satellites %6.4f km \n',max_r);
fprintf('Min distance between two satellites %6.4f km \n',min_r);
Max distance between two satellites 13850.3054 km 
Min distance between two satellites 262.0271 km

Relative_motion_sim_01 Relative_motion_sim_02 Relative_motion_sim_03

Relative Motion of Satellites

clc; clear all; close all;

% This example shows how to use the state vectors of spacecraft A and B
% to find the position, velocity and acceleration of B relative
% to A in the LVLH frame attached to A. For more details about the theory and
% algorithm look at Chapter 7 of  H. D. Curtis, Orbital Mechanics for
% Engineering Students,Second Edition, Elsevier 2010

Input

State vectors of Satellite A and B

RA = [-266.77,  3865.8, 5426.2];     % [km]
VA = [-6.4836, -3.6198, 2.4156];     % [km/s]
RB = [-5890.7, -2979.8, 1792.2];     % [km]
VB = [0.93583, -5.2403, -5.5009];    % [km/s]

% Earth gravitational parameter
mu = 398600;                        % [km^3/s^2]

Algorithm

hA = cross(RA, VA);                 % Angular momentum of A
% Unit vectors i, j,k of the co-moving frame
i = RA/norm(RA);  k = hA/norm(hA); j = cross(k,i);

% Transformation matrix Qxx:
QXx   = [i; j; k];
Om    = hA/norm(RA)^2;
Om_dt = -2*VA*RA'/norm(RA)^2.*Om;

% Accelerations of A and B,inertial frame
aA = -mu*RA/norm(RA)^3;
aB = -mu*RB/norm(RB)^3;

% Relative position,inertial frame
Rr = RB - RA;
% Relative position,LVLH frame attached to A
R_BA = QXx*Rr';

% Relative velocity,inertial frame
Vr = VB - VA - cross(Om,Rr);
% Relative velocity,LVLH frame attached to A
V_BA = QXx*Vr';

% Relative acceleration, inertial frame
ar = aB - aA - cross(Om_dt,Rr) - cross(Om,cross(Om,Rr))- 2*cross(Om,Vr);
% Relative acceleration,LVLH frame attached to A
a_BA = QXx*ar';

fprintf('Position of B relative to A in LVLH frame attached to A \n');
fprintf('R_BA = [%4.2f %4.2f %4.2f] km \n\n', R_BA);

fprintf('Velocity of B relative to A in LVLH frame attached to A \n');
fprintf('V_BA = [%6.4f %6.4f %6.4f] km/s \n\n', V_BA);

fprintf('Acceleration of B relative to A in LVLH frame attached to A \n');
fprintf('a_BA = [%8.8f %8.8f %8.8f] km/s^2 \n', a_BA);
Position of B relative to A in LVLH frame attached to A 
R_BA = [-6701.22 6828.28 -406.24] km 

Velocity of B relative to A in LVLH frame attached to A 
V_BA = [0.3168 0.1120 1.2470] km/s 

Acceleration of B relative to A in LVLH frame attached to A 
a_BA = [-0.00022213 -0.00018083 0.00050590] km/s^2

 

Orbital Inclination Change

Transfer from a LEO 350 km circular orbit with 53.4 deg inclination to a Geostationary Equatorial Orbit(GEO)

clear all; clc;
close all;

Input

R_LEO = 6378 + 350;      % km
R_GEO = 42164;           % km
mu    = 398600;          % km^3/s^2
incl  = 53.4;            % deg

Method 1: Hohman transfer from LEO to GEO and after inclination change

Rp  = R_LEO;
Ra  = R_GEO;
e   = (Ra - Rp)/(Ra + Rp);  % transfer orbit eccentricity
a   = (Ra + Rp)/2;          % transfer orbit semimajor axis
V_LEO   = (mu/R_LEO)^0.5;
Vp      = (2*mu/R_LEO - mu/a)^0.5;
Va      = (2*mu/R_GEO - mu/a)^0.5;
V_GEO   = (mu/R_GEO)^0.5;
dV_LEO  = abs(Vp - V_LEO);
dV_GEO  = abs(V_GEO - Va);
dV_Hoff = dV_GEO + dV_LEO;

% Inclination change at GEO
dV_incl = 2*V_GEO*sind(incl/2);
dV_total1 = dV_incl + dV_Hoff;

fprintf('Method 1: Hohman transfer from LEO to GEO and after inclination change\n');
fprintf('dV_Hoff  = %6.2f [km/s] \n',dV_Hoff);
fprintf('dV_incl  = %6.2f [km/s] \n',dV_incl);
fprintf('dV_total = %6.2f [km/s] \n',dV_total1);
Method 1: Hohman transfer from LEO to GEO and after inclination change
dV_Hoff  =   3.87 [km/s] 
dV_incl  =   2.76 [km/s] 
dV_total =   6.64 [km/s]

Method 2: Inclination change in LEO and after Hohman transfer to GEO

dV_incl   = 2*V_LEO*sind(incl/2);
dV_total2 = dV_incl + dV_Hoff;

fprintf('\nMethod 2: Inclination change in LEO and after Hohman transfer to GEO\n');
fprintf('dV_incl  = %6.2f [km/s] \n',dV_incl);
fprintf('dV_Hoff  = %6.2f [km/s] \n',dV_Hoff);
fprintf('dV_total = %6.2f [km/s] \n\n',dV_total2);
fprintf('Method 2/Method 1 = %6.2f \n',dV_total2/dV_total1);
Method 2: Inclination change in LEO and after Hohman transfer to GEO
dV_incl  =   6.92 [km/s] 
dV_Hoff  =   3.87 [km/s] 
dV_total =  10.79 [km/s] 

Method 2/Method 1 =   1.63

Kepler Equation Solver Without Transcendental Function Evaluations

clc;clear all;

tic;
j = 1;
for M = 0:0.005:pi
    i = 1;
    for e = 0:0.005:1
        c3 = 5/2 + 560*e;
        a = 15/c3*(1-e);
        b = -M/c3;
        y = (b^2/4+ a^3/27)^0.5;
        x = (-b/2 + y)^(1/3) - (b/2 + y)^(1/3);
        c15 = 3003/14336 + 16384*e;
        c13 = 3465/13312 - 61440*e;
        c11 = 945/2816 + 92160*e;
        c9 = 175/384 -70400*e;
        c7 = 75/112 + 28800*e;
        c5 = 9/8-6048*e;
        x2 = x^2;x3 = x2*x;x4 = x3*x;
        x5 = x4*x;x6 = x5*x;x7 = x6*x;
        x8 = x7*x;x9 = x8*x;
        x10 =x9*x; x11 = x10*x; x12 = x11*x;
        x13 = x12*x; x14 = x13*x; x15 = x14*x;

        f = c15*x15 + c13*x13 + c11*x11 + c9*x9 + c7*x7 + c5*x5 + c3*x3 +...
            15*(1 - e)*x -M;
        f1 = 15* c15 *x14 + 13* c13 *x12 + 11* c11 *x10 + 9* c9 *x8 + 7* c7 *x6 + ...
            5 *c5 *x4 + 3 *c3 *x2 + 15*(1 - e);
        f2 = 210* c15 *x13 + 156* c13 *x11 + 110* c11 *x9 + 72* c9 *x7 + ...
            42* c7 *x5 + 20* c5 *x3 + 6* c3 *x;
        f3 = 2730*c15 *x12 + 1716* c13 *x10 + 990* c11 *x8 + 504* c9 *x6 +...
            210* c7 *x4 + 60* c5 *x2 + 6* c3;
        f4 = 32760* c15 *x11 + 17160* c13 *x9 + 7920* c11 *x7 + 3024* c9 *x5 +...
            840* c7 *x3 + 120* c5 *x;
        f5 = 360360* c15 *x10 + 154440* c13 *x8 + 55440* c11 *x6 + 15120* c9 *x4 ...
            + 2520* c7 *x2 + 120* c5;
        f6 = 3603600* c15 *x9 + 1235520* c13 *x7 + 332640* c11 *x5 + 60480* c9 *x3...
            + 5040* c7 *x;
        f7 = 32432400* c15 *x8 + 8648640* c13 *x6 + 1663200* c11 *x4 + 181440* c9 *x2 ...
            + 5040* c7;
        f8 = 25945920* c15 *x7 + 51891840* c13 *x5 + 6652800* c11 *x3 + ...
            362880* c9 *x;
        f9 = 1.8162144e9* c15 *x6 + 259459200* c13 *x4 + 19958400* c11 *x2 + ...
            362880* c9;
        f10 = 1.08972864e10* c15 *x5 + 1.0378368e9*c13 *x3 + 39916800* c11 *x;
        f11 = 5.4486432e10* c15 *x4 + 3.1135104e9* c13 *x2 + 39916800* c11;
        f12 = 2.17945728e11* c15 *x3 + 6.2270208e9* c13 *x;
        f13 = 6.53837184e11* c15 *x2 + 6.2270208e9* c13;
        f14 = 1.307674368e13* c15 *x;
        f15 = 1.307674368e13* c15;

        g1 = 1/2; g2 = 1/6; g3 = 1/24; g4 = 1/120;
        g5 = 1/720; g6 = 1/5040; g7 = 1/40320; g8 = 1/362880;
        g9 = 1/3628800; g10 = 1/39916800; g11 = 1/479001600;
        g12 =1/6.2270208e9 ; g13 = 1/8.71782912e10; g14 = 1/1.307674368e12;

        u1 = -f/f1;
        h2 = f1 + g1* u1* f2;
        u2 = -f/h2;
        h3 = f1 + g1* u2* f2 + g2* u2^2*f3;
        u3 = -f/h3;
        h4 = f1 + g1* u3* f2 + g2* u3^2*f3 + g3* u3^3*f4;
        u4 = -f/h4;
        h5 = f1 + g1* u4* f2 + g2* u4^2*f3 + g3* u4^3*f4 + g4* u4^4*f5;
        u5 = -f/h5;
        h6 = f1 + g1* u5* f2 + g2* u5^2*f3 + g3* u5^3*f4 + g4* u5^4*f5 +...
            g5* u5^5*f6;
        u6 = -f/h6;
        h7 = f1 + g1* u6* f2 + g2* u6^2*f3 + g3* u6^3*f4 + g4* u6^4*f5 +...
            g5* u6^5*f6 + g6* u6^6*f7;
        u7 = -f/h7;
        h8 = f1 + g1* u7* f2 + g2* u7^2*f3 + g3* u7^3*f4 + g4* u7^4*f5 +...
            g5* u7^5*f6 + g6* u7^6*f7 + g7* u7^7*f8;
        u8 = -f/h8;
        h9 = f1 + g1* u8* f2 + g2* u8^2*f3 + g3* u8^3*f4 + g4* u8^4*f5 +...
            g5* u8^5*f6 + g6* u8^6*f7+ g7* u8^7*f8 + g8* u8^8*f9;
        u9 = - f/h9;
        h10 = f1 + g1* u9* f2 + g2* u9^2*f3 + g3* u9^3*f4 + g4* u9^4*f5 +...
            g5* u9^5*f6 + g6* u9^6*f7+ g7* u9^7*f8 + g8* u9^8*f9 + ...
            g9* u9^9*f10;
        u10 = -f/h10;
        h11 = f1 + g1* u10* f2 + g2* u10^2* f3 + g3* u10^3* f4 + g4* u10^4*f5...
            + g5* u10^5* f6 + g6* u10^6*f7+ g7* u10^7* f8 + g8* u10^8* f9...
            + g9* u10^9* f10 + g10* u10^10* f11;
        u11 = -f/h11;
        h12 = f1 + g1* u11* f2 + g2* u11^2* f3 + g3* u11^3* f4 + ...
            g4* u11^4* f5 + g5* u11^5* f6 + g6* u11^6* f7...
        + g7* u11^7* f8 + g8* u11^8* f9 + g9* u11^9* f10 + g10* u11^10* f11...
        + g11* u11^11* f12;
        u12 = -f/h12;
        h13 = f1 + g1* u12* f2 + g2* u12^2* f3 + g3* u12^3* f4 +...
            g4* u12^4* f5 + g5* u12^5* f6 + g6* u12^6* f7...
        + g7* u12^7* f8 + g8* u12^8* f9 + g9* u12^9* f10 +...
        g10* u12^10* f11 + g11* u12^11* f12 + g12* u12^12* f13;
        u13 = -f/h13;
        h14 = f1 + g1* u13* f2 + g2* u13^2* f3 + g3* u13^3* f4 + ...
            g4* u13^4* f5 + g5* u13^5* f6 + g6* u13^6* f7...
        + g7* u13^7* f8 + g8* u13^8* f9 + g9* u13^9* f10 + g10* u13^10*f11...
        + g11* u13^11* f12+ g12* u13^12* f13 + g13* u13^13* f14;
        u14 = -f/h14;
        h15 = f1 + g1* u14* f2 + g2* u14^2* f3 + g3* u14^3* f4 + ...
            g4* u14^4* f5 + g5* u14^5* f6 + g6* u14^6* f7...
        + g7* u14^7* f8 + g8* u14^8* f9 + g9* u14^9* f10 +...
        g10* u14^10* f11 + g11*u14^11* f12+ g12* u14^12* f13 + ...
        g13* u14^13* f14 + g14* u14^14* f15;
        u15 = -f/h15;

        x =  x + u15;
        w = x - 0.01171875*x^17/(1 + e);
        E = M + e*(-16384*w^15 + 61440*w^13 - 92160*w^11+ 70400*w^9 -...
            28800*w^7 + 6048*w^5 - 560*w^3 + 15*w);
                  E_nit(i,j) = E;
                i = i + 1;
    end
    j = j + 1;
end
t_nit = toc;
fprintf('Calculation time = %4.6f [s] \n',t_nit);
Calculation time = 2.591335 [s] 
kepler

Newton Iterative Solver

tic;
j = 1;
eps = 1E-15;        % Tolerance
for M = 0:0.005:pi
    i = 1;
    for e = 0:0.005:1
        En  = M;
        Ens = En - (En-e*sin(En)- M)/(1 - e*cos(En));
        while ( abs(Ens-En) > eps )
            En = Ens;
            Ens = En - (En - e*sin(En) - M)/(1 - e*cos(En));
        end;
        E_it(i,j) = Ens;
        i = i + 1;
    end
    j = j + 1;
end
t_it = toc;
fprintf('Calculation time = %4.6f [s] \n',t_it);
Calculation time = 3.520123 [s]
Error = log10(abs(E_it - E_nit));
e = 0:0.005:1;
M = (0:0.005:pi);
contourf (M,e,Error);
xlabel('M - Mean anomaly [rad]');
ylabel('e - Eccentricity');
title('Error in E calculation on Log scale [rad] ');
colorbar;

 REFERENCES
 Adonis Pimienta-Pe˜nalver and John L. Crassidis†,ACCURATE KEPLER EQUATION
 SOLVER WITHOUT TRANSCENDENTAL FUNCTION EVALUATIONS,
 University at Buffalo, State University of New York, Amherst, NY, 14260-4400

 

%d bloggers like this: