Small Satellites

Home » Featured (Page 2)

Category Archives: Featured

Advertisements

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

Phased Array Antenna, Radiation Pattern and Array Configuration


Phased array consisting of N lined up individual isotropic antennas with the composite main antenna beam into vertical direction. This code based on the governing equation (J. Röttger, The Instrumental Principles of MST Radars, Ch.2.1) and shows how the radiation pattern changes as the parameters are modified.

clc; clear all;
close all;

Ratio between wavelength and distance between individual elements Linear Array of N isotropic radiators

N    = 64;
alfa = -90:0.1:90;       % Zenith angle
d_l  = [0.1,0.5,1,2,5] ; % Ratio
for i = 1:5
    Ed = abs(sin(N*pi*d_l(i)*sind(alfa))./sin(pi*d_l(i)*sind(alfa)))/N;
    fig = figure(i);
    set(fig,'Position', [100, 100, 1049, 400]);
    subplot(2,2,[1 3]);
    plot(alfa,Ed);
    grid on;
    axis([-90 90 0 1]);
    xlabel('Zenith angle [deg]');
    ylabel('Normalized radiation pattern');
    title(['Ratio between distance and wavelength d/\lambda = ',num2str(d_l(i))]);

    subplot(2,2,[2 4]);
    polar(alfa*pi/180,Ed);
    hold on;
    polar((alfa+180)*pi/180,Ed);
    xlabel(['Polar plot for the radiation pattern d/\lambda = ',num2str(d_l(i))]);
    hold off;
end

PhasedArray_01 PhasedArray_02 PhasedArray_03 PhasedArray_04 PhasedArray_05

Distance between individual elements

f_EISCAT =  233E6;   % Centre frequency 233 Hz
c = 3*10^8;          % Light Speed m/s
l = c/f_EISCAT;
d = [0.1,0.5,1,2,3];  % m
for i = 1:5
    % Normalizied gain
    Ed = abs(sin(N*pi*d(i)/l*sind(alfa))./sin(pi*d(i)/l*sind(alfa)))/N;
    % Ed - Array Factor, Radiation patern for isotropic radiators
    fig = figure(i+5);
    set(fig,'Position', [100, 100, 1049, 400]);
    subplot(2,2,[1 3]);
    plot(alfa,Ed);
    grid on;
    axis([-90 90 0 1]);
    xlabel('Zenith angle [deg]');
    ylabel('Normalized radiation pattern');
    title([' Distance between individual elements (space weighting) = ',...
        num2str(d (i)),' m,  F = 233Mhz']);
        subplot(2,2,[2 4]);
    polar(alfa*pi/180,Ed);
    hold on;
    polar((alfa+180)*pi/180,Ed);
    xlabel('Polar plot for the radiation pattern');
    hold off;
 end

PhasedArray_06 PhasedArray_07 PhasedArray_08 PhasedArray_09 PhasedArray_10

Number of antenna elements The distance between individual elements equals to the half wavelength of the signal

d_l = 0.5;
N   = [4 16 64 100 169]; % Number of elements was varied from 4 to 256.
for i = 1:5
    Ed = abs(sin(pi*d_l*N(i)*sind(alfa))./sin(pi*d_l*sind(alfa)))/N(i);
    fig = figure(i+10);
    set(fig,'Position', [100, 100, 1049, 400]);
    subplot(2,2,[1 3]);
    plot(alfa,Ed);
    grid on;
    axis([-90 90 0 1]);
    xlabel('Zenith angle [deg]');
    ylabel('Normalized array factor,Radiation pattern');
    title(['Number of elements = ',num2str(N(i)),', d/\lambda = 0.5']);
    subplot(2,2,[2 4]);
    polar(alfa*pi/180,Ed);
    hold on;
    polar((alfa+180)*pi/180,Ed);
    xlabel('Polar plot for the radiation pattern');
    hold off;
end

PhasedArray_11 PhasedArray_12 PhasedArray_13 PhasedArray_14 PhasedArray_15

Radiation pattern change when the beam is not pointed vertically

clc;clear all;
N = 8;
d_l = 0.5;
alfa0 = 30;
alfa = -180:0.1:180;       % Zenith angle
Ed = abs(sin(pi*d_l*N*(sind(alfa) - sind(alfa0) ))./...
                        sin(pi*d_l*(sind(alfa) - sind(alfa0))))/N;
Ed0 = abs(sin(pi*d_l*N*(sind(alfa)))./...
                        sin(pi*d_l*sind(alfa)))/N;
Edl = 20*log10(Ed);
Edl0= 20*log10(Ed0);
    fig = figure(16);
    set(fig,'Position', [100, 100, 1049, 400]);
    subplot(2,2,[1 3]);
    plot(alfa,Ed);
    hold on;
    plot(alfa,Ed0,'r');
    grid on;
    axis([-90 90 0 1]);
    xlabel('Zenith angle [deg]');
    ylabel('Normalized array factor,Radiation pattern');
    title(['Number of elements = ',num2str(N),', d/\lambda = 0.5']);
    legend('\beta = 30 deg','\beta = 0 deg','Location','NorthWest')
    subplot(2,2,[2 4]);
    polar(alfa*pi/180,Ed);
    hold on;
    polar(alfa*pi/180,Ed0,'r');
    hold on;
    xlabel('Polar plot for the radiation pattern');
    hold off;
    figure(17);
    plot(alfa,Edl);
    hold on;
    plot(alfa,Edl0,'r');
    grid on;
    legend('\beta = 30 deg','\beta = 0 deg','Location','NorthWest')
    axis([-90 90 -100 0]);
    xlabel('Zenith angle [deg]');
    ylabel('Radiation pattern dB');
    title(['Number of elements = ',num2str(N),', d/\lambda = 0.5']);
PhasedArray_16

PhasedArray_17

Electrical weighting techniques

N     = 64;
d_l   = 0.1;
Ed = abs(sin(pi*d_l*N*(sind(alfa) ))./...
                        sin(pi*d_l*sind(alfa)))/N;
Edl = 20*log10(Ed);
% Try to lower the first side lobe relative to main lobe
% To be able to reduce the side lobe levels, we design the array so its radiate
% more power towards the center, and less at the edges.
w  = zeros(1,N);
ws = zeros(1,N);
n = 1:N;
E = 0;Es = 0;
for n = 1:N;
    if n <= N/2
        w(n) = (n-1)/30.5;
    else
        w(n) = w(65 - n);
    end
    ws(n) = sin((n-1)*pi/(N-1));
    E  = E + (w(n)*exp(complex(0,((2*pi*(n-1)*d_l)*sind(alfa)))));
    Es  = Es + (ws(n)*exp(complex(0,((2*pi*(n-1)*d_l)*sind(alfa)))));
end
figure(19);
hold on; grid on;
plot(w,'g');
plot(ws,'r');
axis([1 64 0 1]);
ylabel('Normalized tapering weight');
xlabel('Array elements index');
legend('polynomial','sin');
hold off;
E = abs(E);
E = 20*log10(E/max(E));
Es = abs(Es);
Es = 20*log10(Es/max(Es));
    fig = figure(18);
    set(fig,'Position', [100, 100, 600, 700]);
    hold on;
    subplot(3,1,1);
    plot(alfa,Edl);
    grid on;
    axis([-90 90 -100 0]);
    legend('Isotropic');
    ylabel('Radiation pattern dB');
    title(['Number of elements = ',num2str(N),', d/\lambda = ',num2str(d_l)]);
    subplot(3,1,2);
    plot(alfa,Es,'r');
    axis([-90 90 -100 0]);
    ylabel('Radiation pattern dB');
    grid on;
    legend('Weighted Sin');
    subplot(3,1,3);
    plot(alfa,E,'g');
    axis([-90 90 -100 0]);
    grid on;
    xlabel('Zenith angle [deg]');
    ylabel('Radiation pattern dB');
    legend('Weighted Poly');
    hold off;

PhasedArray_18 PhasedArray_19

 

Universal anomaly, Lagrange f and g coefficients

Contents

Universal anomaly, Lagrange f and g coefficients

Given the initial position R0 and velocity V0 vector of a satelite at time t0 in earth centered earth fixed reference frame. In this example we show how to calculate position and velocity vectors of the satellite dt minutes later using Lagrange f and g coefficients in terms of the universal anomaly

clc; clear all;
R0 = [7200, -13200 0];     %[km]   Initial position
V0 = [3.500, 2.500 1.2];   %[km/s] Initial velocity
dt = 36000;                %[s] Time interval

mu = 398600;               %[km^3/s^2] Earth’s gravitational parameter
r0  = norm(R0);
v0  = norm(V0);
vr0 = R0*V0'/r0;
alfa   = 2/r0 - v0^2/mu;

if     (alfa > 0 )
    fprintf('The trajectory is an ellipse\n');
elseif (alfa < 0)
    fprintf('The trajectory is an hyperbola\n');
elseif (alfa == 0)
    fprintf('The trajectory is an parabola\n');
end

X0 = mu^0.5*abs(alfa)*dt;  %[km^0.5]Initial estimate of X0
Xi = X0;
tol = 1E-10;              % Tolerance
while(1)
    zi = alfa*Xi^2;
    [ Cz,Sz] = Stumpff( zi );
    fX  = r0*vr0/(mu)^0.5*Xi^2*Cz + (1 - alfa*r0)*Xi^3*Sz + r0*Xi -(mu)^0.5*dt;
    fdX = r0*vr0/(mu)^0.5*Xi*(1 - alfa*Xi^2*Sz) + (1 - alfa*r0)*Xi^2*Cz + r0;
    eps = fX/fdX;
    Xi = Xi - eps;
    if(abs(eps) < tol )
        break
    end
end
fprintf('Universal anomaly X = %4.3f [km^0.5] \n\n',Xi)
% Lagrange f and g coefficients in terms of the universal anomaly

f  =  1 - Xi^2/r0*Cz;
g  = dt - 1/mu^0.5*Xi^3*Sz;
R  = f*R0 + g*V0;
r = norm(R);
df = mu^0.5/(r*r0)*(alfa*Xi^3*Sz - Xi);
dg = 1 - Xi^2/r*Cz;
V = df*R0 + dg*V0;
fprintf('Position after %4.2f [min] R = %4.2f*i + %4.2f*j + %4.2f*k[km] \n',dt/60,R);
fprintf('Velocity after %4.2f [min] V = %4.3f*i + %4.3f*j + %4.2f*k[km/s] \n',dt/60,V);
The trajectory is an ellipse
Universal anomaly X = 1922.210 [km^0.5] 

Position after 600.00 [min] R = -6781.27*i + -11870.72*j + -3270.69*k[km] 
Velocity after 600.00 [min] V = 3.488*i + -3.362*j + 0.41*k[km/s]

Ploting trajectory

t      = 0;
step   = 100;     %[s] Time step
ind    = 1;
while (t < dt)
    [R V] = UniversalLagrange(R0, V0, t);
    Ri(ind,:) = R;
    [Long(ind,:), Lat(ind,:)] = R2RA_Dec(R);
    ind = ind + 1;
    t = t + step;
end
% Earth 3D Plot
Earth3DPlot(1);
scatter3(Ri(:,1),Ri(:,2),Ri(:,3),'.r');
hold off;
% Earth topographic map
EarthTopographicMap(2,820,420);
scatter(Long ,Lat,'.r');
text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
 title('Satellite ground track');
 hold off;
%

Lagrange Universal Anomaly_02

Lagrange Universal Anomaly_01

Hyperbolic trajectory

clear Ri Long Lat
R0 = [7200, -6200 0];      %[km]   Initial position
V0 = [5.500, 7.500 1.2];   %[km/s] Initial velocity
dt = 12000;                %[s] Time interval

% Ploting trajectory
t      = 0;
step   = 100;     %[s] Time step
ind    = 1;
while (t < dt)
    [R V] = UniversalLagrange(R0, V0, t);
    Ri(ind,:) = R;
    [Long(ind,:), Lat(ind,:)] = R2RA_Dec(R);
    ind = ind + 1;
    t = t + step;
end
% Earth 3D Plot
Earth3DPlot(3);
scatter3(Ri(:,1),Ri(:,2),Ri(:,3),'.r');
% Earth topographic map
EarthTopographicMap(4,820,420);
scatter(Long ,Lat,'.r');
text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
'HorizontalAlignment','left','FontSize',14 );
 title('Satellite ground track');

Lagrange Universal Anomaly_03Lagrange Universal Anomaly_04

Equations from the book Orbital Mechanics for Engineering Students, Second Edition Aerospace_Engineering

 

Satellite Ground Track, MOLNIYA 1-93

Molniya satellite systems were military communications satellites used by the Soviet Union. The satellites used highly eccentric elliptical orbits of 63.4 degrees inclination and orbital period of about 12 hours. This type of orbits allowed satellites to be visible to polar regions for long periods.

In this example we will show how to plot ground track for satellites.

clear all;
clc;
% Earth topographic map
figure(1);
xwidth = 820;
ywidth = 420;
hFig = figure(1);
 set(gcf,'PaperPositionMode','auto')
 set(hFig, 'Position', [100 100 xwidth ywidth])
hold on;
grid 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',[0 360],'YLim',[-90 90], ...
    'XTick',[0 60 120 180 240 300 360], ...
    'Ytick',[-90 -60 -30 0 30 60 90]);
image([0 360],[-90 90],topo,'CDataMapping', 'scaled');
colormap(topomap1);
ylabel('Latitude [deg]');
xlabel('Longitude [deg]');
title('MOLNIYA 1-93 satellite ground track');

R_e = 6378;        % Earth's radius
mu = 398600;       % Earth’s gravitational parameter [km^3/s^2]
J2 = 0.0010836;
we = 360*(1 + 1/365.25)/(3600*24);      % Earth's rotation [deg/s]
% MOLNIYA 1-93 Orbital Elements /Source: AFSPC
rp    =  1523.9 + R_e;       % [km] Perigee Radius
ra    =  38843.1 + R_e;      % [km] Apogee Radius
theta =  0;                  % [deg] True anomaly
RAAN  =  141.4992;           % [deg] Right ascension of the ascending node
i     =  64.7;               % [deg] Inclination
omega =  253.3915 ;          % [deg] Argument of perigee

a = (ra+rp)/2;               % Semimajor axis
e = (ra -rp)/(ra+rp) ;       % Eccentricity
h = (mu*rp*(1 + e))^0.5;     % Angular momentum
T = 2*pi*a^1.5/mu^0.5;       % Period
dRAAN = -(1.5*mu^0.5*J2*R_e^2/((1-e^2)*a^3.5))*cosd(i)*180/pi;
domega = dRAAN*(2.5*sind(i)^2 - 2)/cosd(i);
% Initial state
[R0 V0] = Orbital2State( h, i, RAAN, e,omega,theta);
[ alfa0 ,delta0 ] = R2RA_Dec( R0 );
scatter(alfa0,delta0,'*k');

ind = 1;
eps = 1E-9;
dt = 20;        % time step [sec]
ti = 0;

while(ti <= 2.5*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;
    [ alfa(ind) ,delta(ind) ] = R2RA_Dec( R );
    ti = ti+dt;
    ind = ind + 1;
end
scatter(alfa,delta,'.r');
text(280,-80,'smallsats.org','Color',[1 1 1], 'VerticalAlignment','middle',...
	'HorizontalAlignment','left','FontSize',14 );

Molniya Ground Track

%d bloggers like this: