clear all;
close all;

%% graben1.txt should contain a list of the locations and MIT gravimeter
%% readings for flag points 02W - 36E, in order. I created the file by hand
%% from graben1_utm.txt (the output of readgps.m) and a file of gravimeter
%% readings.

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

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

%% 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)];
   cntr = [cntr; 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, 4020100*ones(length(n),1) ,name);
set(t, 'FontSize', [6]);
set(t, 'Rotation', -80);

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

%% plot a cross-section 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 west best fit line

wbfl = polyfit(d(1:3), h(1:3), 1);
wbfx = [d(1):5:d(3)+20];
wbfy = polyval(wbfl, wbfx);
plot(wbfx, wbfy, 'r');

%% compute and draw the basin best fit line

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

%% compute and draw the east best fit line

ebfl = polyfit(d(26:32), h(26:32), 1);
ebfx = [d(26)-20:5:d(32)+20];
ebfy = polyval(ebfl, ebfx);
plot(ebfx, ebfy, 'b');

title('Side view of main E/W graben traversal');
xlabel('Horizontal distance (m) from 02W');
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 = 10:14,
   anom = model2(2.2, 2.7, dip, d, h, meas, 3);
   plot(d, anom, 'r'), hold on, grid;
   text(d(size(d))+20, anom(size(d)), sprintf('%d', dip));
end

title('Expected Bouguer anomaly with 10 to 14-degree rock dips (\rho_1 = 2.2, \rho_2 = 2.7)');
xlabel('Distance (m) from 02W');
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, 12, d, h, meas, 3);
   plot(d, anom, 'r'), hold on, grid;
   text(d(size(d))+20, anom(size(d)), sprintf('%3.1f', rho1));
end

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

subplot(3,1,3);

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

title('Expected Bouguer anomaly with 2.5 to 2.9 rock density (\rho_1 = 2.2, dip = 12)');
xlabel('Distance (m) from 02W');
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, 3);
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 02W');
ylabel('Elevation (m)');

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

figure; clf;

for dip = 7:0.5:10,
   anom = model2(2.4, 3.0, dip, d, h, meas, 3);
   plot(d, anom, 'r:'), hold on;
   text(d(size(d))+20, anom(size(d)), sprintf('%4.1f', dip));
end

dip = bestdip(2.4, 3.0, 8:0.01:9, d, h, meas, 3);
anom = model2(2.4, 3.0, dip, d, h, meas, 3);
plot(d, anom, 'b'), grid;
text(d(size(d))+20, anom(size(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 02W');
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, 7:0.1:24, d, h, meas, 3)];
   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, 8:23);
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)');