clear all;
close all;

%% graben2.txt should contain a list of the locations and UMASS gravimeter
%% readings (converted) for flag points 00 - 08n, in order. I created the
%% file by hand from graben2_utm.txt (the output of readgps.m) and a file
%% of gravity data.

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

%% make some empty arrays to read data into
name = []; e = []; n = []; d = []; h = []; meas = [];

%% read the data in a loop, appending rows to the arrays
while ~feof(fgra)
   name = str2mat(name, fscanf(fgra,'%s',1));
   n = [n; fscanf(fgra,'%d',1)];
   e = [e; fscanf(fgra,'%d',1)];
   d = [d; fscanf(fgra,'%lf',1)];
   h = [h; fscanf(fgra,'%lf',1)];
   meas = [meas; fscanf(fgra,'%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(n)),:);

fclose(fgra);

%% convert gravity measurements to mgals
%% NOTE! This is for the MIT instrument ONLY 
%% meas = (cntr - 3100)*1.05725 + 3272.0;


%% draw a top view of the flag path

figure; clf; subplot(2,1,1);
plot(e,n,'o'), hold on;
t = text(e+100, n ,name);
set(t, 'FontSize', [6]);

bfl = polyfit(e, n, 1);
plot(e, polyval(bfl, e), 'r'), axis('equal'), grid;
title('Map view of NE/SW graben traversal');
xlabel('Easting (m)');
ylabel('Northing (m)');

%% plot a side view to put best-fit lines on

subplot(2,1,2);
plot(d,h,'o'), hold on, grid;
t = text(d,h-5,name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', -80);

%% compute and draw the elevation best fit line

bbfl = polyfit(d, h, 1);
bbfx = [min(d)-20:5:max(d)+20];
bbfy = polyval(bbfl, bbfx);
plot(bbfx, bbfy, 'g');

title('Side view of NE/SW graben traversal');
xlabel('Horizontal distance (m) from 00');
ylabel('Elevation (m)');


%% fiddle with the basal sliding model
%% bouganom = model2(rho1, rho2, dip, dist, ht, meas, base);

%% try a few different angles

figure; clf; subplot(3,1,1);

for dip = 6:10,
   anom = model2(2.2, 2.7, dip, d, h, meas, 1);
   plot(d, anom, 'r'), hold on, grid;
   text(d(length(d))+20, anom(length(d)), sprintf('%d', dip));
end

title('Expected Bouguer anomaly with 6 to 10-degree rock dips (\rho_1 = 2.2, \rho_2 = 2.7)');
xlabel('Distance (m) from 00');
ylabel('Bouguer anomaly (mgal)');
 
%% now try varying the densities

subplot(3,1,2);

for rho1 = 2.0:0.1:2.4,
   anom = model2(rho1, 2.7, 8, d, h, meas, 1);
   plot(d, anom, 'r'), hold on, grid;
   text(d(length(d))+20, anom(length(d)), sprintf('%3.1f', rho1));
end

title('Expected Bouguer anomaly with 2.0 to 2.4 sediment density (\rho_2 = 2.7, dip = 8)');
xlabel('Distance (m) from 00');
ylabel('Bouguer anomaly (mgal)');

subplot(3,1,3);

for rho2 = 2.5:0.1:2.9,
   anom = model2(2.2, rho2, 8, d, h, meas, 1);
   plot(d, anom, 'r'), hold on, grid;
   text(d(length(d))+20, anom(length(d)), sprintf('%3.1f', rho2));
end

title('Expected Bouguer anomaly with 2.5 to 2.9 rock density (\rho_1 = 2.2, dip = 8)');
xlabel('Distance (m) from 00');
ylabel('Bouguer anomaly (mgal)');

%% let's see if we can use the Bouguer anomaly to determine the rock shape

rho1 = 2.2; rho2 = 2.7;
anom = model2(rho1, rho2, 0, d, h, meas, 1);
delta_h = anom / (2*pi*6.673e-8*(rho2-rho1)*1e5);
rockht = delta_h + 300;

figure;
plot(d, rockht, '+'), hold on;
plot(d, h, 'x'), axis('equal'), grid;
text(500,600,'surface');
text(500,500,'sediment');
text(500,300,'rock');

title(sprintf('Predicted rock profile for \\rho_1=%3.1f, \\rho_2=%3.1f', ...
              rho1, rho2));
xlabel('Distance (m) from 00');
ylabel('Elevation (m)');

%% now try to find the best angle for different density values

figure; clf;

for dip = 4:0.5:8,
   anom = model2(2.4, 3.0, dip, d, h, meas, 1);
   plot(d, anom, 'r:'), hold on;
   text(d(length(d))+20, anom(length(d)), sprintf('%4.1f', dip));
end

dip = bestdip(2.4, 3.0, 4:0.01:8, d, h, meas, 1);
anom = model2(2.4, 3.0, dip, d, h, meas, 1);
plot(d, anom, 'b'), grid;
text(d(length(d))+20, anom(length(d)), sprintf('%4.2f', dip));

title('Example of best (min squared anomaly) rock dip (\rho_1 = 2.4, \rho_2 = 3.0)');
xlabel('Distance (m) from 00');
ylabel('Bouguer anomaly (mgal)');

figure; clf; subplot(2,1,1);

rhostep = 0.05;
alldips = [];
for rho2 = 2.6:rhostep:2.9,
   dip = [];
   for rho1 = 2.0:rhostep:2.4,
      dip = [dip; bestdip(rho1, rho2, 4:0.1:18, d, h, meas, 1)];
   end
   alldips = [alldips dip];
   plot(2.0:rhostep:2.4, dip, 'r'), grid, hold on;
   text(2.41, dip(length(dip)), sprintf('%4.2f',rho2));
end

title('Best rock dip for different densities');
xlabel('sediment density \rho_1  (g/cm^3)');
ylabel('best rock dip (degrees)');

subplot(2,1,2);
[foo, bar] = contour(2.6:rhostep:2.9, 2.0:rhostep:2.4, alldips, 4:18);
clabel(foo, bar), grid;
title('Best rock dip for different densities');
xlabel('rock density \rho_2  (g/cm^3)');
ylabel('sediment density \rho_1  (g/cm^3)');