function Leica(filepath)

% Reads a text file with Leica points and generates contour and surface
% plots of the surveyed dune.
%
% Michel Louge 04-26-2013
%
%	input variables
%           -	filepath    =	path of the text file
%
%	output variables
%           -   xx          =   xx
%
%   global variables
%           -   none        =   none
%
% Example: Leica('DunePoints.txt');
%

%
% 	other internal variables
%
%           -   fid         =   file index
%           -   data_cell   =   cell structure with data
%           -   dune        =   dune points [Northing Easting Height]
%           -   npoint      =   number of grid points
%           -   x_min       =   minimum Easting (m)
%           -   x_max       =   maximum Easting (m)
%           -   x_range     =   range of Easting (m)
%           -   x_step      =   step of Easting (m)
%           -   x_grid      =   regular Easting grid (m)
%           -   y_min       =   minimum Northing (m)
%           -   y_max       =   maximum Northing (m)
%           -   y_range     =   range of Northing (m)
%           -   y_step      =   step of Northing (m)
%           -   y_grid      =   regular Northing grid (m)
%           -   FF          =   interpolation function of dune profile
%           -   z_grid      =   interpolated dune profile at grid points
%           -   h           =   plot variable
%           -   N0          =   Northing of recovered probe (m)
%           -   E0          =   Easting of recovered probe (m)
%           -   H0          =   height of recovered probe (m)
%           -   alpha       =   bearing of dune (rad)
%           -   col         =   color line vector
%

%% Set parameters:
%
npoint      =   500     ; % number of grid points to describe the dune

%% Read input file:
%
fid         =   fopen(filepath)                             ;
data_cell   =   textscan(fid,'%9.4f %9.4f %9.4f',...
                    'delimiter','\t','HeaderLines',0)       ;
fclose(fid)                                                 ;

dune        =   [data_cell{1} data_cell{2} data_cell{3}]    ;

%% Set up the regular grid for plotting:
%
x_min       =   min(dune(:,2))  ;
x_max       =   max(dune(:,2))  ;
x_range     =   x_max-x_min     ;
x_step      =   x_range/npoint  ;

y_min       =   min(dune(:,1))  ;
y_max       =   max(dune(:,1))  ;
y_range     =   y_max-y_min     ;
y_step      =   y_range/npoint  ;

% setup the mesh grid with increasing x and y 
% (increasing is crucial to find outward normal):

[x_grid,y_grid] = meshgrid(x_min:x_step:x_max,y_min:y_step:y_max)   ;    

%% Interpolate profile

%
FF    	= TriScatteredInterp(dune(:,2),dune(:,1),dune(:,3))	; 
z_grid  = FF(x_grid,y_grid)	; % interpolated dune profile at grid points

%% Plot
%

scrsz   = get(0,'ScreenSize')       ; % finds screen size of machine
figure('Position',...
 [1 scrsz(4) scrsz(3) scrsz(4)]) ; % reserves screen space for figures
                                   % ScreenSize is a four-element vector:
                                   % [left, bottom, width, height]

contour3(x_grid,y_grid,z_grid,8) % superimpose contour on surface
h = findobj('Type','patch');
set(h,'LineWidth',2)
hold on
surf(x_grid,y_grid,z_grid,'FaceAlpha',.8) % for surface alone
daspect([1 1 0.3])  % set aspect ratio of surface plot for aethestics
xlabel('Easting Coord. (m)','fontsize',16)
ylabel('Northing Coord. (m)','fontsize',16)
zlabel('Height (m)','fontsize',16)
shading interp
colormap hsv
colorbar
hold on

%% Mark general direction of dune motion and position of recovered probe
%
% from Leica file M_DUNE0413.TXT of April 2013:
%
%BP     -73.4063	37.7983	4.2067
%
N0      =   -73.4063                ; % coordinates of the recovered probe
E0      =   37.7983                 ;
H0      =   4.2067                  ;
%
dd      =   25                      ; % estimated length of the vector
alpha   =   160/180*pi              ; % general bearing of the dune
col     =   [61/256 61/256 61/256]  ; % color for plot
%quiver3(E0,N0,H0,...
%        dd*sin(alpha),dd*cos(alpha),0,...
%        'LineWidth',2,...
%        'Color',col)
quiver3(E0+dd*sin(alpha),N0+dd*cos(alpha),H0,...
        -dd*sin(alpha)*1.5,-dd*cos(alpha)*1.5,0,...
        'LineWidth',2,...
        'Color',col) % reverse vector
% mark the location of the uncovered RH/T probe
plot3(E0,N0,H0,'o',...
               'MarkerEdgeColor',col,...
               'MarkerFaceColor',col,...
               'MarkerSize',6)



