%  1.  NAME
%
%  gps - Compare position from different navigation devices
%
%  2.  USAGE
%
%  gps
%
%  3.  DESCRIPTION
%
%  This Matlab script reads navigation data from two different files,
%  assuming that these result from two different navigation devices.
%  Typically, this could be the navigation positions stored in a CODAS
%  database, and an external GPS device.
%
%  Input files must have the following columns:
%
%  1. Time in decimal days
%
%  2. Longitude in degrees
%
%  3. Latitude in degrees
%
%  This script is quite crude and could be improved easily.
%
%  4.  PARAMETERS
%
%  The following parameters must be edited:
%
%     year_base
%        Base year used in input files.
%
%     nav_file
%        Name of first navigation file.
%
%     gps_file
%        Name of second GPS data file.
%
%     avg_time
%        Ensemble average time in seconds. This parameter is used to find
%        common positions and gaps between both input files.
%
%     plot_title
%        Title for plot.
%
%     note
%        Additionnal string to print on the plot.
%
%  5.  FILES
%
%  o  gps.ps
%
%  6.  AUTHORS
%
%  Pierre Jaccard, Geophysical Institute, University of Bergen, 1999
%

year_base = 1997;
nav_file  = 'getnav.nav';
gps_file  = 'azhtech.gps';
avg_time  = 300;
plot_title = 'AGMASCO';
note       = 'P. Jaccard';

%------------------------------------------------------------------
% No more editing!!!

%--- Read navigation data:
nav = read_asc(nav_file);
if size(nav,2) < 3,
	error([ 'Insuffiscient number of columns in file ' nav_file ]);
end

%--- Read gps data:
gps = read_asc(gps_file);
if size(gps,2) < 3,
	error([ 'Insuffiscient number of columns in file ' gps_file ]);
end

%--- Initialize indexes:
nav_t = 1; nav_lon = 2; nav_lat = 3;
gps_t = 1; gps_lon = 2; gps_lat = 3;
if size(gps,2) > 3,
	gps_fix = 4;
else
	gps_fix = -1;
end

%--- Get common positions and missing points:
nav = nav';
gps = gps';
eps = (avg_time * 0.2)/86400.0;
ic = [];
ii  = 1:length(gps(gps_t,:));
for i=ii,
	ic = [ ic, find(abs(nav(nav_t,:) - gps(gps_t,i)) < eps) ];
end
ibad = ones(size(nav(nav_t,:)));
ibad(ic) = 0;
ibad = find(ibad == 1);
if gps_fix > 0,
	ibad = [ ibad, find(gps(gps_fix,:) < 0) ];
end

%--- Calculate position difference in meters:
mn_lat = (nav(nav_lat,ic) + gps(gps_lat,:))/2.0;
d_lat  = gps(gps_lat,:) - nav(nav_lat,ic);  
d_lon  = gps(gps_lon,:) - nav(nav_lon,ic);  
dx     = lon_to_m(d_lon, mn_lat);
dy     = lat_to_m(d_lat, mn_lat);

%--- Prepare plot:
figure;
orient tall;

%--- Plot position differences:
subplot(2,1,1);
plot(gps(gps_t,:), dx, '-', gps(gps_t,:), dy, ':');
grid on;
xlabel([ 'Decimal days (' int2str(year_base), ')' ]);
ylabel('GPS - NAV (m) [-]=dx [:]=dy');
if length(ibad) > 0,
	hold on;
	ymin = min([ dx, dy ]);
	bad  = ones(size(ibad)) * ymin;
	plot(nav(nav_t,ibad), bad, 'x');
end

%--- Plot cruise track:
subplot(2,1,2);
plot(gps(gps_lon, :), gps(gps_lat, :), '-', ...
		nav(nav_lon,:), nav(nav_lat,:), ':');
grid on;
xlabel([ 'Decimal days (' int2str(year_base), ')' ]);
ylabel('Cruise Track [-]=GPS [:]=NAV');
if length(ibad) > 0,
	hold on;
	plot(nav(nav_lon,ibad), nav(nav_lat, ibad), 'x');
end

%--- Additionnal stuff:
date_it([ '- NAV File: ' nav_file ' - GPS File: ' gps_file ...
			' - PS File: gps.ps' ]);
old_axes = gca;
aa = axes('position', [0 0 1 1]);
axis('off');

str = plot_title;
text(0.5, 0.99, str, 'units', 'normalized', ...
		'FontAngle', 'italic', 'HorizontalAlignment', 'center', ...
		'VerticalAlignment', 'top', 'FontSize', 12, 'FontWeight', 'bold');

t1 = min(nav(nav_t,:));
t2 = max(nav(nav_t,:));
str = [ to_date(year_base, t1) ' to ' to_date(year_base, t2) ];
text(0.5, 0.96, str, 'units', 'normalized', ...
		'FontAngle', 'italic', 'HorizontalAlignment', 'center', ...
		'VerticalAlignment', 'top', 'FontSize', 10, 'FontWeight', 'bold');
	
str = note;
text(0.05, 0.02, str, 'units', 'normalized', ...
		'FontAngle', 'italic', 'HorizontalAlignment', 'left', ...
		'VerticalAlignment', 'bottom', 'FontSize', 8, 'FontWeight', 'normal');

str = sprintf('dx: mean=%-8.2f sdev=%-8.2f min=%-8.2f max=%8.2f\n', ...
		mean(dx), std(dx), min(dx), max(dx));
str = sprintf('%sdx: mean=%-8.2f sdev=%-8.2f min=%-8.2f max=%8.2f', ...
		str, mean(dy), std(dy), min(dy), max(dy));
text(0.95, 0.02, str, 'units', 'normalized', ...
		'FontAngle', 'normal', 'HorizontalAlignment', 'right', ...
		'VerticalAlignment', 'bottom', 'FontSize', 8, 'FontWeight', 'normal');

print -dps gps.ps
