Home » Posts tagged 'Small Satellites'
Tag Archives: Small Satellites
Mesosphere-Stratosphere-Troposphere(MST) Radar Data Analysis
In this example we analyse data from MST (mesosphere-stratosphere-troposphere) radar observations. MST radars used to do observation of the dynamics of the lower and middle atmosphere to study winds, waves, turbulence and instabilities generate irregularities in the atmosphere. The reflected radar signals from the random irregularities are collected by the receiver antenna. The return signal strength is highly depending on the refractive index which is a function of atmospheric parameters such as humidity, temperature, and pressure and electron density. Hence those parameters will highly affect the signal to noise ratio (SNR). Input data file structure {UT, Altitude,Signal amplitude (linear),Signal-to-noise ratio (SNR), dB, Zonal wind, m/s, Meridional wind, m/s}
clear all; clc; close all; % Data files name ID fl = {'150 m height resolutions'; '1200 m height resolutions'; 'Barker Coding'; 'Complementary Coding'; 'Uncoded Data'; }; % Files for three different observation days fid = { 'TXT_20061211_test1.fca','TXT_20071010_test1.fca',... 'TXT_20081014_test1.fca'; % 150 m height resolutions 'TXT_20061211_test2.fca','TXT_20071010_test2.fca',... 'TXT_20081014_test2.fca'; % 1200 m height resolutions 'TXT_20061211_test3.fca','TXT_20071010_test3.fca',... 'TXT_20081014_test3.fca'; % Barker coding 'TXT_20061211_test4.fca','TXT_20071010_test4.fca',... 'TXT_20081014_test4.fca'; % Complementary coding 'TXT_20061211_test5.fca','TXT_20071010_test5.fca',... 'TXT_20081014_test5.fca'; % Uncoded data }; dates = ['2006 Dec 11'; '2007 Oct 10'; '2008 Dec 14' ];
SNR value as a function of universal time and altitude for three different observation days
fs = size(fid,1); fd = size(fid,2); for jd = 1:fd for id = 1:fs clear SNR; fname = fid{id,jd}; data = load(fname); time = unique(data(:,1)); alt = unique(data(:,2)); size_t = size(time,1); size_a = size(alt,1); % Signal-to-noise ratio (SNR), dB for i = 1:size_a for j = 1:size_t SNR(i,j)= data(i+((j-1)*size_a),4); end end % Plot if(id > 2 ) FigHandle = figure(2); hold on; ii = id - 2; subplot(3,3,ii+(jd-1)*3); colormap(jet); pcolor(time,alt,SNR); shading flat; caxis([-25,30]); if(id == 5) colorbar; end if(id == 3) ylabel('Altitude [km]'); end if((id == 4)) xlabel(['UT/Date: ',dates(jd,:) ]); end if((jd == 1)&&(id == 4) ) title(['Signal-to-Noise Ratio(SNR) [dB]', fl(id)]); else if(jd == 1) title(fl(id)); end end set(FigHandle, 'Position', [100, 0, 800, 800]); else FigHandle = figure(1); hold on; subplot(3,2,id+(jd-1)*2); colormap(jet); pcolor(time,alt,SNR); shading flat; caxis([-25,30]); if(id == 2) colorbar; end if(id == 1) ylabel('Altitude [km]'); end xlabel(['UT/Date: ',dates(jd,:) ]); if((jd == 1)) title(['Signal-to-Noise Ratio(SNR) [dB]', fl(id)]); end set(FigHandle, 'Position', [100, 0, 600, 800]); end end end![]()
Magnitude and direction of horizontal winds
The figures below provide the horizontal wind variation over the altitude and universal time for three different days. The direction of the wind is mostly from west to east. The results are expected as the wind in this layer of atmosphere moves west to east because of the Coriolis acceleration due to force caused by the rotation of the earth.
clear all; clc; close all; dates = ['2006 Dec 11'; '2007 Oct 10'; '2008 Dec 14' ]; % Data files name ID fid = {'TXT_20061211_test4.fca','TXT_20071010_test4.fca',... 'TXT_20081014_test4.fca' }; % Complementary coding % Radar with Complementary coding for id = 1:3 data = load(fid{id}); time = unique(data(:,1)); alt = unique(data(:,2)); size_t = size(time,1); size_a = size(alt,1); % Signal-to-noise ratio (SNR), dB for i = 1:size_a for j = 1:size_t zWind(i,j)= data(i+((j-1)*size_a),5); % Zonal Wind, East mWind(i,j)= data(i+((j-1)*size_a),6); % Meridional Wind, North end end hWind = sqrt(zWind.^2 + zWind.^2); % Horizontal Wind [m/s] %dWind= atan2(zWind,mWind)*180/pi + 180; % Wind Direction angle, Zero in North direction a_max = 60; a_min = 30; FigHandle = figure(2 + id); set(FigHandle, 'Position', [100, 0, 800, 400]); subplot(2,2,[1 3]); colormap(cool); pcolor(time,alt(1:a_max),hWind(1:a_max,:)); shading flat; colorbar; xlabel(['UT/Date:',dates(id,:)]); ylabel('Altitude [km]'); title('Horizontal wind speed [m/s], Complementary coded signal'); subplot(2,2,[2 4]); hold on; whitebg([0.0 .0 .2]); quiver(time,alt(1:a_max),zWind(1:a_max,:),mWind(1:a_max,:),'r'); %contour(time,alt(1:a_max),hWind(1:a_max,:)); xlabel(['UT/Date:',dates(id,:)]); title('Horizontal wind direction(top - North,right - East)'); axis([min(time),max(time),min(alt(1:a_max)),max(alt(1:a_max))]); hold off; end
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
Relative Motion of Satellites
clc; clear all; close all;
% 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. 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
Input
State vectors of Satellite A and B
RA = [-266.77, 3865.8, 5426.2]; % [km] VA = [-6.4836, -3.6198, 2.4156]; % [km/s] RB = [-5890.7, -2979.8, 1792.2]; % [km] VB = [0.93583, -5.2403, -5.5009]; % [km/s] % Earth gravitational parameter mu = 398600; % [km^3/s^2]
Algorithm
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 = QXx*Rr'; % Relative velocity,inertial frame Vr = VB - VA - cross(Om,Rr); % Relative velocity,LVLH frame attached to A V_BA = QXx*Vr'; % Relative acceleration, inertial frame ar = aB - aA - cross(Om_dt,Rr) - cross(Om,cross(Om,Rr))- 2*cross(Om,Vr); % Relative acceleration,LVLH frame attached to A a_BA = QXx*ar'; fprintf('Position of B relative to A in LVLH frame attached to A \n'); fprintf('R_BA = [%4.2f %4.2f %4.2f] km \n\n', R_BA); fprintf('Velocity of B relative to A in LVLH frame attached to A \n'); fprintf('V_BA = [%6.4f %6.4f %6.4f] km/s \n\n', V_BA); fprintf('Acceleration of B relative to A in LVLH frame attached to A \n'); fprintf('a_BA = [%8.8f %8.8f %8.8f] km/s^2 \n', a_BA);
Position of B relative to A in LVLH frame attached to A R_BA = [-6701.22 6828.28 -406.24] km Velocity of B relative to A in LVLH frame attached to A V_BA = [0.3168 0.1120 1.2470] km/s Acceleration of B relative to A in LVLH frame attached to A a_BA = [-0.00022213 -0.00018083 0.00050590] km/s^2
Orbital Inclination Change
Transfer from a LEO 350 km circular orbit with 53.4 deg inclination to a Geostationary Equatorial Orbit(GEO)
clear all; clc; close all;
Input
R_LEO = 6378 + 350; % km R_GEO = 42164; % km mu = 398600; % km^3/s^2 incl = 53.4; % deg
Method 1: Hohman transfer from LEO to GEO and after inclination change
Rp = R_LEO; Ra = R_GEO; e = (Ra - Rp)/(Ra + Rp); % transfer orbit eccentricity a = (Ra + Rp)/2; % transfer orbit semimajor axis V_LEO = (mu/R_LEO)^0.5; Vp = (2*mu/R_LEO - mu/a)^0.5; Va = (2*mu/R_GEO - mu/a)^0.5; V_GEO = (mu/R_GEO)^0.5; dV_LEO = abs(Vp - V_LEO); dV_GEO = abs(V_GEO - Va); dV_Hoff = dV_GEO + dV_LEO; % Inclination change at GEO dV_incl = 2*V_GEO*sind(incl/2); dV_total1 = dV_incl + dV_Hoff; fprintf('Method 1: Hohman transfer from LEO to GEO and after inclination change\n'); fprintf('dV_Hoff = %6.2f [km/s] \n',dV_Hoff); fprintf('dV_incl = %6.2f [km/s] \n',dV_incl); fprintf('dV_total = %6.2f [km/s] \n',dV_total1);
Method 1: Hohman transfer from LEO to GEO and after inclination change dV_Hoff = 3.87 [km/s] dV_incl = 2.76 [km/s] dV_total = 6.64 [km/s]
Method 2: Inclination change in LEO and after Hohman transfer to GEO
dV_incl = 2*V_LEO*sind(incl/2); dV_total2 = dV_incl + dV_Hoff; fprintf('\nMethod 2: Inclination change in LEO and after Hohman transfer to GEO\n'); fprintf('dV_incl = %6.2f [km/s] \n',dV_incl); fprintf('dV_Hoff = %6.2f [km/s] \n',dV_Hoff); fprintf('dV_total = %6.2f [km/s] \n\n',dV_total2); fprintf('Method 2/Method 1 = %6.2f \n',dV_total2/dV_total1);
Method 2: Inclination change in LEO and after Hohman transfer to GEO dV_incl = 6.92 [km/s] dV_Hoff = 3.87 [km/s] dV_total = 10.79 [km/s] Method 2/Method 1 = 1.63