%  ECE 6390 radiolocation project 1
%  Team:  The Gibbs Phenomenon
clear
format long
% load data from file
load DopplerData.txt;
freqDat = DopplerData(2,:);
freqDat = MYLPF(freqDat);  %  low-pass filter allows for better
                           %  differentiation.
timeDat = DopplerData(1,:);

%--------- find center frequency---------------------------------------
[a,b] = size(freqDat);

%  Center frequency is halfway between the first and last frequencies
%  instead of as a hard-coded number because the oscillator on the 
%  EPIRB might not be perfect.
centerFreq = (freqDat(1)+freqDat(b))/2

%  The following two lines determine the point in the frequency data vector
%  that is closest to the center frequency.  The point at the closest
%  frequency is considered to mark the latitude of the transmitter because
%  that is where the doppler shift effects are zero.
latloc = abs(centerFreq - freqDat);

[Y,I] = min(latloc)         % Y is min value, I is index of value

hitTime = timeDat(I);   % This is the time at which the latitude of the
                        % satellite is that of the transmitter

T = 6119;                   % Orbital period (known for SARSAT system)
satVel = 360/T*cos(9*pi/180);   %  Satellite orbital velocity in
                                %  degrees per second.

%  Now, find latitude by multiplying orbital velocity in degrees per second
%  by the hitTime, which is the amount of time between the satellite
%  passing the equator and the doppler effect being zero.
lat = satVel * hitTime      % lattitude

%---------------------------Now to find the Longitude---------------------

%  Derivative of the received frequency at hit time.  Note that I use the
%  min operator because the derivative is negative.  I divide by 0.1 to put
%  the answer in Hz/second.
derSig = min(diff(freqDat))/0.1;                            
                            
%  wavelength of transmitted EPIRB signal
lambda = 3e8/centerFreq;

re = 6370e3;                % radius of the earth(m)
alt = 850e3;                % altitude of satellite(m) (measured surface)

%  Convert satellite velocity from degrees/second to meters/second
Vy = satVel*(2*pi*(re+alt)/360);       % Velocity of satellite in m/s

%  Really need a picture to explain this.  The following two lines
%  determine the straight-line distance from the transmitter on the earth's
%  surface with the satellite in orbit at the satellite's closest point to
%  the transmitter.  It comes from imagining a triangle in space with one
%  vertex at the transmitter, one at the point where the satellite is
%  closest to the transmitter, and one point at the location of the
%  satellite one second later.  Theta is the angle at the vertex where the
%  transmitter.

%                                      -----------------x<sat after 1 sec
%                    ----------------                   |
%       -theta-----                                     |  <  Vy
% x-----------------------------------------------------x
% ^transmitter               dist                       ^- sat at latitude

%  theta is found using the doppler shift formula
%  fc + deltaf = fc + cos(something with lambda -- look it up)

theta = acos((derSig*lambda)/Vy) - pi/2;

dist = Vy/tan(theta);      % distance from satellite to beacon

%  I just print these to make sure they're reasonable.
dist
re
alt

%  Another hard one.  Here, I imagine a giant triangle with vertices at the
%  transmitter, the satellite, and the core of the earth.  I then use the
%  law of cosines.

%                                       x<satellite
%                                     / |
%                              dist /   |  alt
%                                  / ---|-----<surface of earth
%                   transmitter>x----   |
%                                  \    |
%                                   \   |re
%                                 re \  |
%                                     \ |
%                                      \|
%                                       x<core of earth
%  re = radius of earth

deltaLong = acos((dist^2 - re^2 - (re + alt)^2)/(-2*re*(re+alt)))
%  deltaLong is in radians



deltaLong = deltaLong*180/pi;       %convert to degrees

long = -76-deltaLong