%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Radiolocation Scavanger Hunt
% Assignment One - 'SARSAT Rescue!'
% Georgia Institute of Technology
% ECE 6390, Fall 2005
% Solution by Team Space Invaders
%	Melissa Amoros
%	Dean Lewis
%	Joel Prothro
%	Tushar Thrivikraman
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Purpose:
%	To locate the position of an EPIRB distress beacon based on a set
%	of dopplar shift measurements from a SARSAT satellite passing
%	overhead.
% Input:
%	Doppler shift data - expected in 'DopplerData.txt'
%	Longitudinal line the receiving SARSAT satellite is following
% Output:
%	Latitude/Longitude position of the EPIRB
% Requires:
%	'DopplerData.txt' - input file containing received frequencies and
%		the time each measurement was taken
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Initializa MatLab
clear
close all
%clc
format compact

% Define Constants
lon_sat = -76; % Longitudinal line the satellite is following (degrees)
Hsat = 850;	   % Height of the satellite above the Earth's surface (km)
Vx = 1.2;      % Longitudinal velocity of satellite (km/s)
Vy = 7.3;	   % Latitudinal velocity of satellite (km/s)
Fc = 406e6;	   % SATSAT carrier frequency (Hz)
lamda = 3e5/Fc;% SARSAT carrier wavelength (km)
Re = 6378.1;   % Radius of Earth (km)
dt = 0.1;	   % Time between frequency measurements (sec)

% Load doppler shift data
data = load('DopplerData.txt');


%%%%% Part 1 : Find EPIRB Latitude %%%%%

% Partition data into frequency and time data sets
max_n = 9650;
n = 1:(max_n-1);
F = data(2,n);
t = data(1,n);
figure(1), plot(t,F); % Plot received frequency vs time

figure(2), plot(t(2:numel(t)),-diff(F)/dt);

% The data suffers from high quantization error. Therefore, it must be
% sampled to produce a useful derivative plot
k = 7; % Sampling interval
i = 1:k:max_n; % Indexes of the samples to use

% Pull out samples from original data sets
sF = F(i); % sampled frequency
st = t(i); % sampled time

% Compute the derivate dF/dt
df = -diff(sF)/(dt*k);

% Plot the sampled derivative versus sampled time
figure(3), plot(st(2:numel(st)), df);

% Find the index of the maximum value of the derivative
max_df_ind = find(df==max(df),1);

% Lookup the time the maximum derivative occured at
max_time = st(max_df_ind);

% It's given that time t=0s occurs when the satellite crosses the equator
% The distance the satellite has traveled is computed from d=vt
dist_km = Vy*max_time; % Distance north of the equator that the satellite
	% has travelled when it encounters the maximum doppler shift

% Convert distance north to degrees latitude
lat_EPIRB = dist_km*360/(2*pi*Re)


%%%%% Part 2 : Generate Known Data Sets %%%%%

Lship = lat_EPIRB; % Latitude of EPIRB (taken from Part 1) (degrees)

i=40;	 % number of longitude points to generate a data set for
j=10000; % number of data points to generate for each data set

% Generate EPIRB Locations
Xship = 50*[0:i]; % generate a data set for i points left/right of the
	% satellite's longitude line, starting directly beneath the satellite,
	% spaced every 50 km
Dship = sqrt(Xship.^2+Hsat^2)'; % Calculate the distance to the satellite
	% from the EPIRB, accounting for satellite height

% Points to calculate known received frequency at
n = [-j/2:j/2-1];

% Calculate the latitudinal position of the satellite at each point
Ysat = n*dt*Vy;

% For each known position of the EPIRB
for k=1:length(Dship),
	% For each satellite position north of the equator
    for l=1:length(Ysat),
		% Calculate the recieved frequency
        theta = atan2(Dship(k),Ysat(l));
			% angle between satellites flight path and the EPIRB
        Fr(k,l)=Fc-Vy*cos(theta)/lamda; % receive frequency
    end
end

% Plot the received frequencies for known distances per vs time
figure(4), plot(n*dt,Fr);

% For each known EPIRB position, compute the derivative of the received
% frequency and the maximum of this derivative
for k=1:length(Dship),
    dFr(k,:) = -diff(Fr(k,:))/dt;
    MAXdFr(k) = max(dFr(k,:));
end

% Plot the derivative of the received frequencies vs time
figure(5), plot(n(2:length(n))*dt,dFr);

% Compute and plot the maximum derivative versus distance from the satellite
deltadf=-1*diff(MAXdFr)./diff(Xship);
figure(6), plot(Xship(2:numel(Xship)),deltadf);


%%%%% Part 3 : Find EPIRB Longitude %%%%%

% Calculate the km per degree longitude
km_degree_long = 2*pi*Re*cosd(lat_EPIRB)/360;

% Calculate the difference between the max amplitude of the derivative (when
% the EPIRB is known to to be directly under the satellite) and the max
% amplitude the derivative of the actual data (found in Part 1)
derv_amp_diff = (max(MAXdFr)-df(max_df_ind));

% MAXdFr from Part 1 was found to be 70, which fell between a distance of
% 600km and 650km. From the deltadf vs Xship curve, the appropriate deltadf
% value to use is 0.03837, which also happens to be the max of deltadf. So
% use this linear factor to convert the difference in amplitudes to a
% difference in distance
dist = derv_amp_diff/max(deltadf);

% Convert distance in km to distance in longitude
distdeg=dist./km_degree_long;

% Calculate the two possible longitudes
long_EPIRB1 = lon_sat-distdeg
long_EPIRB2 = lon_sat+distdeg

lat_EPIRB =
   25.6473
long_EPIRB1 =
  -79.4849
long_EPIRB2 =
  -72.5151

Published with MATLAB® 7.0.4