Home » Posts tagged 'Molniya orbit'
Tag Archives: Molniya orbit
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)
