%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                        %%
%%  brad_grav.m                                           %%
%%                                                        %%
%%  This will plot the gravity data that Brad sent to me  %%
%%                                                        %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



clear
close all


%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Declare variables  %%
%%%%%%%%%%%%%%%%%%%%%%%%%
grid_size                                           = 100;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Open file stream and get rid of the first line  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
grav_file_stream                                    = fopen('USGS.txt', 'r');
trash                                               = fgets(grav_file_stream);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Loop over each line and read in the string   %%
%%  I think 2224 is the correct number of lines  (2141 original)  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for cnt = 1 : 2141
   data_string                                      = fgetl(grav_file_stream);
   lat_deg(cnt)                                     = str2double(data_string(15:16));
   lat_deg_dec_min(cnt)                             = str2double(data_string(19:24));
   lon_deg(cnt)                                     = str2double(data_string(26:28));
   lon_deg_dec_min(cnt)                             = str2double(data_string(29:32));
   obs_grav(cnt)                                    = str2double(data_string(41:48))-979500.;
   faa_grav(cnt)                                    = str2double(data_string(49:56));
   sba_grav(cnt)                                    = str2double(data_string(57:64));
   %iso_grav(cnt)                                    = str2double(data_string(92:98));
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Close the file stream before we finish  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fclose(grav_file_stream);



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Convert to decimal degrees  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lat_deg_dec_per                                     = lat_deg_dec_min ./ 60;
lon_deg_dec_per                                     = lon_deg_dec_min ./ 60;
lat_deg_dec                                         = lat_deg + lat_deg_dec_per;
lon_deg_dec                                         = lon_deg + lon_deg_dec_per;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Convert west longitude to east longitude  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
lon_deg_dec                                         = 360 - lon_deg_dec;



%%%%%%%%%%%%%%%%%%%%
%%  Grid results  %%
%%%%%%%%%%%%%%%%%%%%
%min_x                                               = min(lon_deg_dec);
%max_x                                               = max(lon_deg_dec);
%min_y                                               = min(lat_deg_dec);
%max_y                                               = max(lat_deg_dec);
min_x = 244.45
max_x = 244.6;
min_y = 35.74;
max_y = 35.84;
xvec                                                = linspace(min_x, max_x, grid_size);
yvec                                                = linspace(min_y, max_y, grid_size);
[xmat, ymat]                                        = meshgrid(xvec, yvec);
[grav_mat]                                          = griddata(lon_deg_dec, lat_deg_dec, sba_grav, xmat, ymat, 'cubic');
%[grav_mat]                                          = griddata(lon_deg_dec, lat_deg_dec, iso_grav, xmat, ymat, 'cubic');



%%%%%%%%%%%%%%%%%%%%%%%%
%%  Plot the results  %%
%%%%%%%%%%%%%%%%%%%%%%%%
figure
surf(xmat, ymat, grav_mat)
view(2)
shading interp
%contour(xmat, ymat, grav_mat)
hold on
plot3(lon_deg_dec, lat_deg_dec, 100 * ones(size(lat_deg_dec)), 'k.')
axis equal
axis tight
axis ([244.45 244.6 35.74 35.84])
    colorbar
xlabel('West longitude')
ylabel('Latitude')
title('Simple Bouguer gravity anomaly, USGS, at CA-NV border [milligal]')    
%title('Observed gravity CA-NV border [milligal-979500]')