Home » Posts tagged 'Satellite Ground Track'
Tag Archives: Satellite Ground Track
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');
Ground track generated by AGI STK10
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; %
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');
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 );
Small Satellites 