Home » 2013

# Yearly Archives: 2013

## Experimental Determination of a Geometric Form Factor in a Lidar Equation

In this example we experimentally determine a geometric form factor for an in-homogeneo atmosphere by using the technique described in the article  Sang Whoe Dho, Young Je Park, and Hong Jin Kong,”Experimental determination of a geometric form factor in a lidar equation for an inhomogeneous atmosphere. Link to the article http://lsrl.kaist.ac.kr/homepage/Publications/Papers_Files/0043.pdf

clear all;close all; clc;

```% Load Lidar observation data
R_site      = 0.418;            % [km]
dR          = 0.03;             % [km] Altitude step
n_pow       = 141;
pow         = 0;
Rmin        = 2.3;
Rmax        = 5.9;
R           = R_site:dR:15;      % Alt [km]
eps         = 0.01;              % Error
Gr_alt      = 0;                 % Overlap Altitude [km]
figN        = 1;                 % Figure Count
n           = 6;                 % Polynomial degree
indm = round((Rmin - R_site)/dR + 1);
indx = round((Rmax - R_site)/dR + 1);
Pow_vec = {raw_mat1 raw_mat2 raw_mat3 raw_mat4};

clear raw_mat1  raw_mat2 raw_mat3 raw_mat4
labels = {'Parallel - Analog';'Parallel - Digital';...
'Perpendicular - Analog';'Perpendicular - Digital'};

for f_ind = 1:1```
```    raw_mat = Pow_vec{1,f_ind};
%Calculating average power profile
pow = 0;
for i=1:n_pow
pow = pow + raw_mat(1:indx,i);
end
pow = pow'/n_pow;                 % Average power profile
plog = log(pow.*R(1:indx).^2);    % Average log power profile
%plog = log(R(1:indx).^2);
% plog = log(R(1:indx).^2)
figure(figN);
plot(pow,R(1:indx),'.');
hold on; grid on;
ylabel('Altitude [km] ');
xlabel('Power');
title(['Average Power Profile(N = 141)' labels(f_ind,1)]);
hold off;
figN = figN + 1;
%     % Average log power profile
%     figure(figN);
%     plot(plog,R(1:indx),'.g');
%     hold on; grid on;
%     ylabel('Altitude [km] ');
%     xlabel('log(P*R^2)');
%     title(['Average Log Power Profile(N = 141) log(P*R^2)' labels(f_ind,1)]);
%     hold off;
%     figN = figN + 1;``` Polynomial fit to the power profile

```    A = polyfit(R(indm:indx),pow(indm:indx),n);
Pint = 0;
for i =1:n+1
Pint  = Pint + A(i)*R(1:indx).^(n+1-i);
end
figure(figN);hold on; grid on;
plot(pow(indm:indx),R(indm:indx),'*b');
plot(Pint,R(1:indx),'.r');
ylabel('Altitude [km] ');
xlabel('Power');
title(['Polynomial fit to average power profile' labels(f_ind,1)]);
hold off;
legend('Actual', 'Fit')
figN = figN + 1;``` Geometric factor derived from average power profile

```    S   =  log(R(1:indx).^2.*pow);
Sint = log(R(1:indx).^2.*Pint);
GR_int = exp(S - Sint);         %Geometric factor,Overlap factor

figure(figN);hold on; grid on;
plot(R(1:indx),GR_int,'.r');
xlabel('Altitude [km] ');
ylabel('Geometric factor G(R)');
title(['Geometric factor derived from average power profile'...
labels(f_ind,1)]);
grid on;
hold off;
figN = figN + 1;``` Calculating geometric form factor for each power profile

```     for j = 1:n_pow
pow = raw_mat(1:indx,j);
pow = pow';
A = polyfit(R(indm:indx),pow(indm:indx),n);
Pint = 0;
for i = 1:n+1
Pint  = Pint + A(i)*R(1:indx).^(n+1-i);
end
S   =  log(R(1:indx).^2.*pow);
Sint = log(R(1:indx).^2.*Pint);
GR_int = exp(S - Sint);     %Geometric factor
for k = 1:indx
if (abs(GR_int(k) - 1.0) < eps )
Gr_alt(j) = R(k);
break;
end
end
end
figure(figN);
hist(Gr_alt,20);
xlabel('Altitude [km]');
ylabel('Number of Power profiles');
title(['G(R) altitude variation for various power profiles',...
labels(f_ind,1)]);
figN = figN + 1;
%     %% Mean altitude overlapping occurs, Sea Level
%     mean_Gr_alt = mean(Gr_alt);
%     std_Gr_alt =  std(Gr_alt);
%     display(labels(f_ind,1));
%     fprintf('Overlap Alt = %4.2f +- %4.2f [km] \n\n',mean_Gr_alt,std_Gr_alt);```
`end` 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};
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);
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);
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
for id = 1:3
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,:));
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
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]);
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

## 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
colormap(topomap1);
% Create the surface.
[x,y,z] = sphere(50);
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```