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)
Small Satellites 
What could I do if I wanted to add the J2 effect to those orbits (only in the omega, of course)
Look at J2 propagator….all the code you need is here
https://smallsats.org/2013/01/20/j2-propagator/
Let me know if you need more help.
Can you explain what constitutes orbital2state?
orbital2state function computes state vectors R and V in the geo-centric equatorial frame of reference using orbital elements. Please look at the code here.
https://smallsats.org/2013/01/17/state-vectors-r-v-from-orbital-elements/