%script to import a .fld file from LISCAD
%and make a contour map.  It will not as of yet do breaklines.

%CHANGE things between here:
%file to load is NETR.dat.txt or whatever
load camel.fld

%The loaded file is in the array  NETR
data = camel;
%refer to columns as:
%data(:,1) is block number
%data(:,2) is 116
%data(:,3) is point number
%data(:,4) is Easting
%data(:,5) is Northing
%data(:,6) is Elevation
%data(:,7) is code

load trench.m

trenche=trench(:,4);
trenchn=trench(:,5);
trenchh=trench(:,6);

%set parameters:
%number of x and y nodes for grid
%make these smaller if you get an error about excceding the max array size
nx = 70;
ny = nx;
%contour interval:
ci = 0.5;

%AND HERE





%Here is the map of points:
figure(1)
clf
plot(data(:,4), data(:,5), 'k+')

%Here is contour map:
figure(2)
clf
e=data(:,4);
n=data(:,5);
h=data(:,6);
estep = abs((min(e)-max(e))/nx);
nstep = abs((min(n)-max(n))/ny);
ee = [min(e):estep:max(e)];
nn = [min(n):nstep:max(n)];
[xi,yi,zi] = griddata(e,n,h,ee,nn');

v = floor(min(h)):ci:max(h);


[c,h] = contour(xi,yi,zi,v);

%label contours
clabel(c,h);
colormap(jet)
hold on

%show data:
plot(e,n,'o')

%show grid:
plot(xi,yi,'+')

xlabel('easting (m)')    
ylabel('northing (m)')   
title('Elevation and observation points')
axis equal





%nice contour map
figure(3)
clf
[c,h] = contour(xi,yi,zi,v);

%label contours
%clabel(c,'manual');
colormap(jet)
hold on

%show data:
plot(e,n,'.')

%%show trench
plot(trenche,trenchn,'k-')

xlabel('easting (m)')    
ylabel('northing (m)')   
title('Elevation and observation points')
axis equal
axis([940 1100 875 1030])
print -depsc contours.eps

%nice view
figure(4)
clf
surfl(xi,yi,zi)
colormap(pink)
view(-135,54)
LIGHTANGLE(100,45)
hold on
%%show trench
plot3(trenche,trenchn,trenchh,'k-')
zlabel('elevation (m)')
xlabel('easting (m)')    
ylabel('northing (m)')

%nice view
figure(5)
clf
surfl(xi,yi,zi)
shading interp
colormap(pink)
view(-135,54)
LIGHTANGLE(100,45)
hold on
%%show trench
plot3(trenche,trenchn,trenchh,'k-')
zlabel('elevation (m)')
xlabel('easting (m)')    
ylabel('northing (m)')