Home » Featured (Page 2)
Category Archives: Featured
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');
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
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
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
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
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']);
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;
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 );