sample_frequency=1;
topo_file_name = 'VIDAL_CA-24000.dem';
%function [x_mat, y_mat, topo_mat, pixel_width] = show_24k_topo(topo_file_name, sample_frequency)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%                                                     %%
%%  show_24k_topo.m                                    %%
%%                                                     %%
%%  This will plot a DEM format 1:24000 USGS topo map  %%
%%                                                     %%
%%  Arguments:                                         %%
%%    topo_file_name  : full name of the topo file     %%
%%    sample_frequency: take every n pixel             %%
%%                                                     %%
%%  Returned variables:                                %%
%%    x_mat           : x coordinate matrix (meters)   %%
%%    y_mat           : y coordinate matrix (meters)   %%
%%    topo_mat        : z coordinate matrix (meters)   %%
%%    pixel_width     : sample size of pixes (meters)  %%
%%                                                     %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Declare variables  %%
%%%%%%%%%%%%%%%%%%%%%%%%%
pixel_spacing                  = 30;
pixel_width                    = sample_frequency * pixel_spacing;

%%%%%%%%%%%%%%%%%%%%%%%%
%%  Read in the data  %%
%%%%%%%%%%%%%%%%%%%%%%%%
[latgrat, longgrat, topo_mat, dem_info]   = usgs24kdem(topo_file_name, sample_frequency);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Announce the internal name of the DEM that we just read in  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp(sprintf('Read DEM file: %s', getfield(dem_info, 'Quadranglename')))

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Build some position matrices assuming the pixel spacing, %%
%%  flipping the y axis                                      %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[topo_mat_rows, topo_mat_cols]            = size(topo_mat);

%%%%%%%%%%%%%%%%%%%%%%%%
%%  Plot the results  %%
%%%%%%%%%%%%%%%%%%%%%%%%
surf(longgrat(364:374,265:275), latgrat(364:374,265:275), topo_mat(364:374,265:275))
axis tight
shading interp
daspect([1 0.8285 1/5])
material shiny
%camlight
xlabel('long [deg]')
ylabel('lat [deg]')
zlabel('z [km]')
title(sprintf('%s, %d meter pixel spacing', getfield(dem_info, 'Quadranglename'), sample_frequency * pixel_spacing))
view(0, 90)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  Open Gravity Datafile  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load gravity_chris;
longs_wgs=gravity_chris(:,1);
lats_wgs=gravity_chris(:,2);
heights=gravity_chris(:,3);
g_boug=gravity_chris(:,4);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  convert WGS84 into NAD27  %%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

longs=longs_wgs+7.899917e-4;
lats=lats_wgs-1.606944e-5;

hold on
for i=1:length(g_boug)
	plot3(longs(i),lats(i),heights(i)+350,'ok','markersize',(g_boug(i)-14.2)*6);
end;
hold off
colorbar