clear all;
close all;

%% gpsdata.txt should be somewhere in matlab's path, preferably the current
%% directory. It should have *only* lines with the following format:
%% Station	Lat(deg min sec)	Long(deg min sec)	Height
%% For example:
%% Base	36	17	54.40000	114	19	45.80000	545.000

fgps = fopen('gpsdata.txt','rt');
if fgps < 0
   error('Hmm, either gpsdata.txt doesn''t exist or you can''t read it...');
end

%% make some empty arrays to read data into
name = []; lad = []; lam = []; las = []; lod = []; lom = []; los = []; ht = [];

%% read the rest of the lines in a loop, appending rows to the arrays
while ~feof(fgps)
   name = str2mat(name, fscanf(fgps,'%s',1));
   lad = [lad; fscanf(fgps,'%d',1)];
   lam = [lam; fscanf(fgps,'%d',1)];
   las = [las; fscanf(fgps,'%lf',1)];
   lod = [lod; fscanf(fgps,'%d',1)];
   lom = [lom; fscanf(fgps,'%d',1)];
   los = [los; fscanf(fgps,'%lf',1)];
   ht = [ht; fscanf(fgps,'%lf',1)];
end

%% the name array will contain an extra blank line at the beginning from
%% converting the initial empty array into a string let's chop it to the
%% length of the other arrays.

name = name(2:(1 + size(lad)),:);

fclose(fgps);

%% now get the UTM coords

[n, e] = dms2utm(lad, lam, las,  lod, lom, los);

%% filter the main E/W graben cross-section points

graben1 = []; g1name = []; g1n = []; g1e = []; g1h = [];

for i = 1:size(name),
   if ~isempty(findstr(name(i,:),'w/')) ...
         | ~isempty(findstr(name(i,:),'e/')) ...
	 & isempty(findstr(name(i,:),'se/'))
      graben1 = [graben1; i];     % collect the point array indices
      g1name = str2mat(g1name, name(i,:));  % point names
      g1n = [g1n; n(i)];                    % northings
      g1e = [g1e; e(i)];                    % eastings
      g1h = [g1h; ht(i)];                   % heights
   end
end

%% remove the initial empty string
g1name = g1name(2:(size(g1n)+1),:);

%% compute and draw the horizontal best fit line

figure;
bfl = polyfit(g1e, g1n, 1);
bfx = [min(g1e):50:max(g1e) + 50];
bfy = polyval(bfl, bfx);
plot(bfx, bfy, 'r'), hold on;
plot(g1e, g1n, 'x'), axis('square'), axis('equal'), grid;
t = text(g1e, 4020210 + 0*g1n, g1name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', 90);
title('Top view of graben cross-section (to scale)');
xlabel('Easting (m)');
ylabel('Northing (m)');

%% project into vertical best-fit plane and show side view
%% note bfl(1) is slope of best fit line, bfl(2) is y-intercept

g1d = (polyval(bfl, g1e) - g1n) / (bfl(1) * sqrt(1 + bfl(1)^2)) + ...
      sqrt((g1e - (polyval(bfl, g1e) - g1n) / bfl(1)).^2 + (g1n - bfl(2)).^2);
g1d = g1d - min(g1d);

figure;
plot(g1d, g1h, 'x'), axis('equal'), hold on, grid;
t = text(g1d, g1h - 25, g1name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', -80);
title('Side view of graben cross-section (corrected to vertical best-fit plane)');
xlabel('Corrected horizontal distance (m) from 02W');
ylabel('Elevation (m)');

%% same graph, not to scale (to show slopes better)

figure;
plot(g1d, g1h, 'x'), hold on, grid;
t = text(g1d, g1h - 3, g1name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', -80);
title('Side view of graben cross-section (corrected to vertical best-fit plane)');
xlabel('Corrected horizontal distance (m) from 02W');
ylabel('Elevation (m)');

%% plot all the points in 2d
figure;
plot(e, n, '.'), axis('equal'), hold on, grid;
t = text(e+75, n-75, name);
set(t, 'FontSize', [4]);
set(t, 'Rotation', -45);
title('All points located with differential GPS (to scale)');
xlabel('Easting (m)');
ylabel('Northing (m)');

%% plot all the points in 3d

figure;
plot3(e, n, ht, '.'), hold on;
plot3(g1e, g1n, g1h, 'r.'), axis('equal'), grid;
t = text(e+75, n-75, ht, name);
set(t, 'FontSize', [4]);
title('All points located with differential GPS (height not to scale)');
xlabel('Easting (m)');
ylabel('Northing (m)');
zlabel('Elevation (m)');

%% output the points to a file

futm = fopen('utmdata.txt','wt');
if futm < 0
   error('Sorry, you can''t open utmdata.txt for writing');
end

for i = 1:size(n),
   fprintf(futm, '%10s', name(i,:));
   fprintf(futm, '\t%d', round(n(i)));
   fprintf(futm, '\t%d', round(e(i)));
   fprintf(futm, '\t%7.3f\n', ht(i));
end

fclose(futm);

%% output the points to a file

fgra = fopen('graben1_utm.txt','wt');
if fgra < 0
   error('Sorry, you can''t open graben1_utm.txt for writing');
end

for i = 1:size(g1n),
   fprintf(fgra, '%10s', g1name(i,:));
   fprintf(fgra, '\t%d', round(g1n(i)));
   fprintf(fgra, '\t%d', round(g1e(i)));
   fprintf(fgra, '\t%7.3f', g1d(i));
   fprintf(fgra, '\t%7.3f\n', g1h(i));
end

fclose(fgra);

%% filter the 00-08N graben cross-section points

graben2 = []; g2name = []; g2n = []; g2e = []; g2h = [];

for i = 1:size(name),
   if (~isempty(findstr(name(i,:),'n/0')) ...         % yuck, total hack
         | ~isempty(findstr(name(i,:),'N/0')) ...     % should do by hand
         | ~isempty(findstr(name(i,:),'n/2')) ...
         | ~isempty(findstr(name(i,:),'N/2')) ...
         | ~isempty(findstr(name(i,:),'00e/'))) ...
         & ~(~isempty(findstr(name(i,:),'01n')) ...
             | ~isempty(findstr(name(i,:),'01N')))

      graben2 = [graben2; i];     % collect the point array indices
      g2name = str2mat(g2name, name(i,:));  % point names
      g2n = [g2n; n(i)];                    % northings
      g2e = [g2e; e(i)];                    % eastings
      g2h = [g2h; ht(i)];                   % heights
   end
end

%% remove the initial empty string
g2name = g2name(2:(size(g2n)+1),:);

%% compute and draw the horizontal best fit line

figure;
bfl = polyfit(g2e, g2n, 1);
bfx = [min(g2e):50:max(g2e) + 50];
bfy = polyval(bfl, bfx);
plot(bfx, bfy, 'r'), hold on;
plot(g2e, g2n, 'x'), axis('square'), axis('equal'), grid;
t = text(g2e, 4020210 + 0*g2n, g2name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', 90);
title('Top view of graben cross-section (to scale)');
xlabel('Easting (m)');
ylabel('Northing (m)');

%% project into vertical best-fit plane and show side view

g2d = (polyval(bfl, g2e) - g2n) / (bfl(1) * sqrt(1 + bfl(1)^2)) + ...
      sqrt((g2e - (polyval(bfl, g2e) - g2n) / bfl(1)).^2 + (g2n - bfl(2)).^2);
g2d = g2d - min(g2d);

figure;
plot(g2d, g2h, 'x'), axis('equal'), hold on, grid;
t = text(g2d, g2h - 25, g2name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', -80);
title('Side view of graben cross-section (corrected to vertical best-fit plane)');
xlabel('Corrected horizontal distance (m) from 00');
ylabel('Elevation (m)');

%% same graph, not to scale (to show slopes better)

figure;
plot(g2d, g2h, 'x'), hold on, grid;
t = text(g2d, g2h - 3, g2name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', -80);
title('Side view of graben cross-section (corrected to vertical best-fit plane)');
xlabel('Corrected horizontal distance (m) from 00');
ylabel('Elevation (m)');

%% output the points to a file

fgra = fopen('graben2_utm.txt','wt');
if fgra < 0
   error('Sorry, you can''t open graben2_utm.txt for writing');
end

for i = 1:size(g2n),
   fprintf(fgra, '%10s', g2name(i,:));
   fprintf(fgra, '\t%d', round(g2n(i)));
   fprintf(fgra, '\t%d', round(g2e(i)));
   fprintf(fgra, '\t%7.3f', g2d(i));
   fprintf(fgra, '\t%7.3f\n', g2h(i));
end

fclose(fgra);