%This function finds the error magnification factor of a specified number of 
%GPS satellites. The EMF measures how much the error of the approximate location 
%of a person using GPS will increase due to inaccuracies of their GPS reciever's
%clock.

%This function requires the other functions:
    %MVNewtonsInputs
    %FCreator
    %DCreator

%The input to this function is the number of GPS satellites.

function Result=SatelliteEMF(NumberOfSatellites)
%Rho is the distance of a satellite to the surface of the Earth. It is fixed at 
%26570 km and the following line creates a constant array of that value. The 
%size of the Rho vector is equal to the amount of satellites.
Rho=26570.*ones(1,NumberOfSatellites);

%Theta is the angle around the equator that a satellite is at. It is between 
%0 and 2pi
%Phi is the angle from the north pole that a satellie is at. It is between 0 and
%pi/2

%The following two lines create initial value for Theta and Phi. The size of the
%Theta and Phi vectors is equal to the number of satellites.
Theta=rand(1,NumberOfSatellites)*2*pi;
Phi=rand(1,NumberOfSatellites)/2*pi;

%The following loop checks the values of Theta to see that the satellites are 
%not close together. The value ThetaFarness measures the percent difference
%between the lowest and highest Theta.
ThetaFarness=0;
while and(ThetaFarness<.7, std(Theta)<1.7)
    Theta=rand(1,NumberOfSatellites)*2*pi;
    ThetaFarness=(max(Theta)-min(Theta))/max(Theta);
end

%This is the same as the above loop except for the value Phi. It makes sure the 
%values for Phi are far apart.
PhiFarness=0;
while and(PhiFarness<.7 ,std(Phi)<.45)
    Phi=rand(1,NumberOfSatellites)/2*pi;
    PhiFarness=(max(Phi)-min(Phi))/max(Phi);
end

%The following three lines convert the vectors for Rho, Theta, and Phi into the
%rectangular coordinates for the satellites. These 3 vectors have the 
%coordinates for the specified amount of satellites. The coordinates for the nth
%satellite are the nth entries in each vector.
A=Rho.*cos(Phi).*cos(Theta);
B=Rho.*cos(Phi).*sin(Theta);
C=Rho.*sin(Phi);

%c is light speed
c=299792.458;

%This is the initial guess that will be put into the Newtons or Gauss Newtons
%method to solve the GPS equations for the satellite positions created earlier 
%in this funciton.
%The initial guess is the position of a hypothetical person on the surface of 
%the Earth.
InitialGuess=[0,0,6370,0.0001];

%R is a vector that calculates the distance of each satellite from the initial 
%guess.
R=sqrt(A.^2+B.^2+(C-6370).^2);

%T calculates the travel times for a signal from a satellite to the person.
T=InitialGuess(4)-R./c;

%The following line adds error to the travel times ranging from -10^-8 to 10^-8
%seconds by multiplying 10^-8 by a random number between -1 and 1 and then 
%a vector of the random values to the travel time vector T.
T=T+10^(-8).*(2*rand(1,NumberOfSatellites)-.5);

%The approximation for the persons position is found using Newton's or Gauss
%Newtons method depending on which is the correct method.
Approximation=MVNewtonsInputs(InitialGuess,A,B,C,T,NumberOfSatellites);

%The forward error is the greatest distance between the approximation for the
%person's location and the initial guess (their actual location). This value is
%in kilometers.
ForwardError=norm((Approximation'-InitialGuess),inf);

%Backward error is the maximum difference in distance caused by the maximum 
%variation in time. This is approximately .003 kilometers.
BackwardError=c*10^(-8);

%The error magnification factor is forward error divided by the backward error.
%The error magnification factor measures how much the position error will 
%increase for an increased time error.
ErrorMagnificationFactor=ForwardError/BackwardError;

%The output of the equation is the forward error and the EMF
Result=[ForwardError, ErrorMagnificationFactor];