function surfaceplot (data, xblocks, yblocks, interpmethod)
%
%surfaceplot
%
%This plots the data as a surface (variables like shadingand colormap
%change the surface color and contrast, etc.
%
%EXAMPLE:
%surface
%figure(7)
%clf
%surfl(ee,nn,zi)
%shading interp
%colormap(gray)
%hold on
%plot3(e,n,hh,'b.')
%view(-30,60)
%xlabel('easting')    
%ylabel('northing')
%zlabel('elevation')   
%title('Cubic interpolation')
%print -depsc vmsurface.eps

interpmethod = 'linear'

%this brings up "help" and gives lots of instructions on griddata
help griddata

%tells griddata what the x and y grid looks like on which 
%it will operate-it's the easting/northing stepsize.
%It needs to be positive, so if it ends up being negative 
%you need to take the absolute value, which is the "abs" command
numsteps=100
estep = abs((min(e)-max(e))/numsteps);
nstep = abs((min(n)-max(n))/numsteps);
ee = [min(e):estep:max(e)];
nn = [min(n):nstep:max(n)];

%This grids the data.
%the 'n means transpose N (switch it from a row vector to a column vector)
% This was orginial  [xi,yi,zi] = griddata(e,n,hh,ee,nn');
[xi,yi,zi] = griddata(e,n,h,ee,nn');

%This prepares the data as a surface using the surfl function from matlab
surface
figure(4)
clf
surfl(ee,nn,zi)
shading interp
colormap(gray)
hold on

%This plots the data and labels the axes, etc.
plot3(e,n,h,'b.')
view(-30,60)
xlabel('easting')    
ylabel('northing')
zlabel('elevation')   
title('Surface Map')
axis equal
%This save the plot as a jpg
print('surfacemaps.jpg','-djpeg');

end