clear;
clf;

fname = input('Which topo file?  ','s');

z = double(imread(sprintf('/morgan/topo/tif/%s.tif',fname),'tif'));

loadstr = sprintf('/home/kana/morgan/topo/elev/%s.elev',fname);
elev = dlmread(loadstr,' ');
xcoord = elev(2,:);   ycoord = elev(1,:);
minelev = elev(3,1);  maxelev = elev(3,2);
z = (z/255)*(maxelev-minelev) + minelev;

[n1 n2] = size( z);
x0 = 0.5*(xcoord(3)+xcoord(4));  x1 = 0.5*(xcoord(1)+xcoord(1));
y0 = 0.5*(ycoord(1)+ycoord(3));  y1 = 0.5*(ycoord(2)+ycoord(4));
dx = (x1-x0)/n2;  dy = (y1-y0)/n1;
x = [x0:dx:x1-dx];
y = [y0:dy:y1-dy];

raysel = input('Do you want to do raytracing from a point?  (1) = yes:  ');
if( raysel == 1)
	raytype = input('shooting angle select: (0)=down only, (1)=up only (2)=both:  ');
	nbhdtype = input('Neighborhood choice: (0)=exact pt, (1)=highest in nbhd., (2)=lowest in nbhd:  ');
	visible = raytrace(z,raytype,nbhdtype);
end

plottype = input('Choose: (1) color = a*z, or (2) contour map: ');

if( plottype == 1)
	surf(y,x,z');
	shading interp;  colorbar;  colormap('hot');
%	vertexag = input('Pick c, where plot is scaled by 1/c: ');
%	axis([1 n1 1 n2 min(min(z)) vertexag*max(max(z))]);
	axis equal
	title( sprintf('%s',fname) ); xlabel('east'); ylabel('north');
	if( raysel == 1)
		colormap('bone');
		[ivis,jvis] = find(visible==1.0);
		[i0,j0] = find(visible==2.0);
		yy=(i0-1)*dy+y0; xx=(j0-1)*dx+x0;
		hold on;
		for iv=1:length(ivis)
			yvis = (ivis(iv)-1)*dy+y0;  xvis = (jvis(iv)-1)*dx+x0;
			plot3(yvis,xvis,z(ivis(iv),jvis(iv)) + 2.0,'r.')
		end;
		plot3(yy,xx,z(i0,j0)+10,'r*')
		plot3([yy,yy],[xx,xx],[0,z(i0,j0)+10],'r-')
		hold off;
	end
%	view( -15, 60)
	rotate3d
else
	b = ones(3,3);  b(2,2) = 4.;  b = b/sum(sum(b));
	z = filter2(b,z);
	contourf(z',20)
	colorbar
	con_int = (maxelev-minelev)/20;
	title( sprintf('%s contour map - contour interval = %5.1f m',fname,con_int), 'FontSize',14 );
	xlabel('east'); ylabel('north');
end

frite = input('Write this figure to an image file? (1)=yes: ');
if( frite == 1)
	delete matlab-out.tif;
	disp('Writing file to matlab-out.tif');
	print -dtiff 'matlab-out.tif';
end
