Small Satellites

Home » Posts tagged 'Communications satellites'

Tag Archives: Communications satellites

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

Molniya Orbit, Satellite Orbits

Molniya orbit is a highly elliptical orbit with an inclination of 63.4 degrees, an argument of perigee of -90 degrees and an orbital period of one half of a sidereal day. Molniya orbits are named after a series of Soviet/Russian Molniya communications satellites which have been using this type of orbit since the mid 1960s. In this example we will show how to calculate and plot satellite orbit around earth based on its orbital elements.

clc;
clear all;
mu = 398600;            % Earth’s gravitational parameter [km^3/s^2]
R_earth = 6378;         % Earth radius [km]

% Molniya Orbit 1
a     = 26554;
e     = 0.72;
i     = 63.4;
omega = -90;
theta = 0;
RAAN  = 0;
h =(a*mu*(1 - e^2))^0.5;
% Calculating initial state vector
[ R0 V0 ] = Orbital2State( h, i, RAAN, e,omega,theta);
r = norm(R0);           % Initial radius  [km]
v = norm(V0);           % Initial speed   [km/s]
T = 2*pi*a^1.5/mu^0.5;  % Orbital period [s]
dt = T/10000;           % time step [s]
t = 0;                  % initial time
% Using fourth-order Runge–Kutta method to solve fundamental equation
% of relative two-body motion
F_r = @(R) -mu/(norm(R)^3)*R;
Rd = V0; R  = R0;
i = 1;
while (t <= T)
    Rv(i,:) = R;
    tv(i) =t;
    k_1 = dt*F_r(R);
    k_2 = dt*F_r(R+0.5*k_1);
    k_3 = dt*F_r(R+0.5*k_2);
    k_4 = dt*F_r(R+k_3);
    Rd  = Rd + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    R   = R + Rd*dt;
    t   = t+dt;
    i = i+1;
end

% Molniya Orbit 2
a     = 26554;
e     = 0.72;
i     = 63.4;
omega = -90;
theta = 0;
RAAN  = 250;
h =(a*mu*(1 - e^2))^0.5;
% Calculating initial state vector
[ R0 V0 ] = Orbital2State( h, i, RAAN, e,omega,theta);
r = norm(R0);           % Initial radius  [km]
v = norm(V0);           % Initial speed   [km/s]
T = 2*pi*a^1.5/mu^0.5;  % Orbital period [s]
dt = T/10000;           % time step [s]
t = 0;                  % initial time
% Using fourth-order Runge–Kutta method to solve fundamental equation
% of relative two-body motion
F_r = @(R) -mu/(norm(R)^3)*R;
Rd = V0; R  = R0;
i = 1;
while (t <= T)
    Rv2(i,:) = R;
    tv(i) =t;
    k_1 = dt*F_r(R);
    k_2 = dt*F_r(R+0.5*k_1);
    k_3 = dt*F_r(R+0.5*k_2);
    k_4 = dt*F_r(R+k_3);
    Rd  = Rd + (1/6)*(k_1+2*k_2+2*k_3+k_4);
    R   = R + Rd*dt;
    t   = t+dt;
    i = i+1;
end
% Plotting
figure('Color',[0 0 0]);
figure(1);
hold on;
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);

% Inertial Frame Axis
Xa  = [radius_earth:100:radius_earth+2500];
Z0  = Xa*0;
plot3(-Xa,Z0,Z0,'r')
plot3(Z0,-Xa,Z0,'y')
plot3(Z0,Z0,Xa,'g')
% Plotting Orbits
plot3(Rv(:,1),Rv(:,2),Rv(:,3));
plot3(Rv2(:,1),Rv2(:,2),Rv2(:,3),'y');
axis square off
view(3)
zoom(2)

Molniya Orbit

%d bloggers like this: