Solutions
%ECE 6390 Group Project Fall 2006
%Team Perihelion
%this program uses a "brute force" method for getting the location from the
%data. It's not exactly the most elegant method in the world, but it's
%fairly fast and can be reused for other data points. The key to this
%program is exact knowledge of the radius of the earth at the receiver
%location. This effectively eliminates one of the four variables, making
%the system of equations have a single solution with only three data
%points.
%the speed of light
c = 299860000;
tic
%latitude/longitudes/pseudoranges of the GPS satellites
lon1 = -62.3256;
lat1 = 7.8520;
pr1 = 0*c;
lon2 = -85.2234;
lat2 = 45.1053;
pr2 = -0.003817986*c;
lon3 = -92.8388;
lat3 = 29.7834;
pr3 = -0.004077077*c;
%the earth's radius
rad = 6380000;
%the altitude of the satellite
satalt = 20200000;
%the radius of the satellite
rgps = rad+satalt;
%the positions in x-y-z space of the GPS satellites
x1 = rgps*cosd(lon1)*cosd(lat1);
y1 = rgps*sind(lon1)*cosd(lat1);
z1 = rgps*sind(lat1);
x2 = rgps*cosd(lon2)*cosd(lat2);
y2 = rgps*sind(lon2)*cosd(lat2);
z2 = rgps*sind(lat2);
x3 = rgps*cosd(lon3)*cosd(lat3);
y3 = rgps*sind(lon3)*cosd(lat3);
z3 = rgps*sind(lat3);
%number of iterations - 8 is way more than enough to converge to a point
iterations = 8;
%number of "partitions" to divide the area into
%as this increases, the time of execution increases with the square, so
%increase cautiously
steps = 50;
%vector of latitudes and longitudes, filled in as the answer converges
minlon = zeros(1,iterations+1);
maxlon = zeros(1,iterations+1);
minlat = zeros(1,iterations+1);
maxlat = zeros(1,iterations+1);
%initial values for the latitude-longitude range. If desired, this can be
%restricted to areas inside the perimeter, but it's unnecessary.
minlat(1) = 30;
minlon(1) = -90;
maxlat(1) = 40;
maxlon(1) = -80;
%run through this loop to "zoom in" on the correct location
for iter = 1:iterations
%reset the "error matrix" to steps X steps
errors = zeros(steps);
for i = 1:steps
for j = 1:steps
%get the test latitude and longitude
testlon = minlon(iter) + i*(maxlon(iter)-minlon(iter))/steps;
testlat = minlat(iter) + j*(maxlat(iter)-minlat(iter))/steps;
%get the x-y-z coordinates for the test point
xg = rad*cosd(testlon)*cosd(testlat);
yg = rad*sind(testlon)*cosd(testlat);
zg = rad*sind(testlat);
%assume the test point is correct, and get the c*tau values
%from each of the satellites. if the point really is correct,
%all these values will be identical.
ctau1 = -(sqrt((x1-xg)^2+(y1-yg)^2+(z1-zg)^2)-pr1);
ctau2 = -(sqrt((x2-xg)^2+(y2-yg)^2+(z2-zg)^2)-pr2);
ctau3 = -(sqrt((x3-xg)^2+(y3-yg)^2+(z3-zg)^2)-pr3);
%take the average of the c*tau values
ctauavg = mean([ctau1 ctau2 ctau3]);
%store the RMS error of the three values in the error matrix
errors(i,j) = sqrt((ctauavg-ctau1)^2+(ctauavg-ctau2)^2+(ctauavg-ctau3)^2);
end
end
%find the point with the smallest error
smallerror = find(errors-min(min(errors)) == 0);
%find the corresponding latitude and longitude
smalllon = mod(smallerror(1),steps);
smalllat = ceil(smallerror(1)/steps);
%get a range from that point of the minimum and maximum longitudes and
%latitudes, which can be used in the next iteration
minlon(iter+1) = minlon(iter) + (maxlon(iter)-minlon(iter))*(smalllon-1.1)/steps;
maxlon(iter+1) = minlon(iter) + (maxlon(iter)-minlon(iter))*(smalllon+1.1)/steps;
minlat(iter+1) = minlat(iter) + (maxlat(iter)-minlat(iter))*(smalllat-1.1)/steps;
maxlat(iter+1) = minlat(iter) + (maxlat(iter)-minlat(iter))*(smalllat+1.1)/steps;
end
toc
%after we're done iterating, show the last values we got, which should be
%the answer.
maxlat(end)
maxlon(end)
%plug these values into http://www.mapquest.com/maps/latlong.adp and see
%what you get!
%northspring
%=the varsity?!?!
%No, of course Prof. Durgin has changed the location to Mary Mac's :)
