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


clear
close all


%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Declare variables  %%
%%%%%%%%%%%%%%%%%%%%%%%%%
grid_size                                           = 50;



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



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Loop over each line and read in the string   %%
%%  I think 2142 is the correct number of lines  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for cnt = 1 : 2142
   data_string                                      = fgetl(grav_file_stream);
   lat_deg(cnt)                                     = str2double(data_string(10:11));
   lat_deg_dec_min(cnt)                             = str2double(data_string(13:17));
   lon_deg(cnt)                                     = str2double(data_string(19:21));
   lon_deg_dec_min(cnt)                             = str2double(data_string(22:27));
   obs_grav(cnt)                                    = str2double(data_string(37:45))-979500.;
   obs_grav(cnt)                                    = str2double(data_string(37:45))-979500.;
   faa_grav(cnt)                                    = str2double(data_string(53:58));
   cba_grav(cnt)                                    = str2double(data_string(84:90));
   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);
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, faa_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
hold on
plot3(lon_deg_dec, lat_deg_dec, 100 * ones(size(lat_deg_dec)), 'k.')
axis equal
axis tight
colorbar
xlabel('West longitude')
ylabel('Latitude')
title('Free air gravity anomaly at CA-NV border [milligal]')    
%title('Observed gravity CA-NV border [milligal-979500]')