% -----------------------------------------------------------------------
% Merge navigation track to ADCP data....

load nav.out;
n = nav(1:length(nav)-1,1)-nav(2:length(nav),1);
which_n = find(n)+1;

nav_x = ones(size(Atim))*NaN;
nav_y = ones(size(Atim))*NaN;

nav_x(nav(which_n,1)) = nav(which_n,4);
nav_y(nav(which_n,1)) = nav(which_n,3);

% -----------------------------------------------------------------------
% Turn bad ADCP readings to NaNs...

n = size(EWv);
bad = EWv > 32700;
for i=1:n(1),
	temp = find(bad(i,:));
	EWv(i, temp) = NaN*zeros(size(EWv(i, temp)));
end;

bad = NSv > 32700;
for i=1:n(1),
	temp = find(bad(i,:));
	NSv(i, temp) = NaN*zeros(size(NSv(i, temp)));
end;

% -----------------------------------------------------------------------
% Eliminate portion of data prior to ship coming on circle...

start = 72; 
endit = length(nav_x)-1;

nav_x = nav_x(start:endit);
nav_y = nav_y(start:endit);
Vdir = Vdir(:,start:endit);
Vmag = Vmag(:,start:endit);
EWv = EWv(:,start:endit);
NSv = NSv(:,start:endit);
Vv = Vv(:,start:endit);
Er = Er(:,start:endit);
Vi = Vi (:,start:endit);
Sew = Sew(start:endit);
Sns = Sns(start:endit);
Svv = Svv(start:endit);
Atim = Atim(start:endit);

% Define each circle...

xc = mean(nav_x);
yc = mean(nav_y);
radial = atan2(nav_x-xc,nav_y-yc);
circ = [1,181,377,568,767,985]

% -----------------------------------------------------------------------
% Compute potential by integrating current along vehicle path...

potential = zeros(size(nav_x));
level = 20;
for j = 1:5,	% loop through the 5 turns accompished in the mission...
	potential(circ(j)+1) = 0;
	i = circ(j)+1;
	
	while (i<circ(j+1)),
		ii = 1;
		if isnan(EWv(level,i)+NSv(level,i)+nav_x(i)+nav_y(i)) == 0,
			while isnan(EWv(level,i+ii)+NSv(level,i+ii)+nav_x(i+ii)+nav_y(i+ii)) > 0,
				ii = ii + 1;
				if (ii + i) > circ(j+1)
					break;
				end;
			end;
			if (ii + i) > circ(j+1)
				break;
			end;
			
			temp = (nav_x(i+ii)-nav_x(i))*(EWv(level,i)+EWv(level,i+ii))/2;
			temp = temp+(nav_y(i+ii)-nav_y(i))*(NSv(level,i)+NSv(level,i+ii))/2;
		end;
		if (ii + i) > circ(j+1)
			break;
		end;
		
		potential(i+ii) = potential(i)+temp;
		i = i+ii;
	end;
	i;
	ii;
end;

% Plot each circle...

figure(1); clf; hold on;
color= ['.r', '.g', '.b', '.c', '.m'];
for j = 1:5
	a = circ(j)+1; b = circ(j+1);
	plot(a:b, potential(a:b), color(2*j-1:2*j));
end;

% Remove vorticity...

for j = 1:5,
	potential(circ(j)+1:circ(j+1))= potential(circ(j)+1:circ(j+1)) - potential(circ(j+1))*(0:(circ(j+1)-(circ(j)+1)))/(circ(j+1)-(circ(j)+1));
end;

% Plot each circle...

color= ['r', 'g', 'b', 'c', 'm'];
for j = 1:5
	a = circ(j)+1; b = circ(j+1);
	plot(a:b, potential(a:b), color(j));
end;
hold off;

% -----------------------------------------------------------------------
% Create an x-y array (horizontal plane) on which to solve Laplace's eqn...

max_size = 51;
pf = zeros(max_size,max_size);
nn = zeros(max_size,max_size);
circle = 1;

for i = circ(circle):circ(circle+1), 
	ii = round((max_size-1)*(nav_x(i) - min(nav_x))/(max(nav_x) - min(nav_x)))+1;
	jj = round((max_size-1)*(nav_y(i) - min(nav_y))/(max(nav_y) - min(nav_y)))+1;
	if isnan(potential(i))==0,
		pf(jj, ii) = pf(jj, ii) + potential(i);
		nn(jj, ii) = nn(jj, ii) + 1;
	end;
end;

pf(find(nn)) = pf(find(nn))./nn(find(nn));
nn(find(nn)) = nn(find(nn))./nn(find(nn));
memory = pf(find(nn));

for i = 1:500,
	% Compute field in center of region...
	pfxp = pf(3:max_size, 2:max_size-1);
	pfxm = pf(1:max_size-2, 2:max_size-1);
	pfyp = pf(2:max_size-1,3:max_size);
	pfym = pf(2:max_size-1,1:max_size-2);
	pf(2:max_size-1,2:max_size-1) = (pfxp + pfxm + pfyp + pfym)/4;

	% Handle boundary...
	pf(1,2:max_size-1) = (pf(2,2:max_size-1)+pf(1,1:max_size-2)+pf(1,3:max_size))/3;
	pf(max_size,2:max_size-1) = (pf(max_size-1,2:max_size-1)+pf(max_size,1:max_size-2)+pf(max_size,3:max_size))/3;
	pf(2:max_size-1,1) = (pf(2:max_size-1,2)+pf(1:max_size-2,1)+pf(3:max_size,1))/3;
	pf(2:max_size-1,max_size) = (pf(2:max_size-1,max_size-1)+pf(1:max_size-2,max_size)+pf(3:max_size,max_size))/3;
	pf(1,1) = (pf(1,2)+pf(2,1))/2;
	pf(max_size,1) = (pf(max_size,2)+pf(max_size-1,1))/2;
	pf(1,max_size) = (pf(1,max_size-1)+pf(2,max_size))/2;
	pf(max_size,max_size) = (pf(max_size,max_size-1)+pf(max_size-1,max_size))/2;

	pf(find(nn))=memory;
end;

image(pf)
% -----------------------------------------------------------------------
% Plot succesive radial, tangential currents around ship path...

for circle = 1:5,
	a = circ(circle);
	b = circ(circle+1);
	
	Cx = ones(38,b-a+1);
	Cy = ones(38,b-a+1);
	Cz = ones(38,b-a+1);
	CRv = ones(38,b-a+1);
	CVv = ones(38,b-a+1);
	CTv = ones(38,b-a+1);
	
	for j = 1:38,
		Cx(j,:) = unwrap(radial(a:b));
		Cy(j,:) = ones(1,b-a+1).*dpth(j);
		CVv(j,:) = Vv(j, a:b);
		CRv(j,:) = EWv(j, a:b).*sin(radial(a:b))+NSv(j, a:b).*cos(radial(a:b));
		CTv(j,:) = EWv(j, a:b).*cos(radial(a:b))-NSv(j, a:b).*sin(radial(a:b));
	end;

	%---  Do separate plots of radial, tangential, and vertical velocity...	
	cla;
	subplot(3,1,1);
	pcolor(Cx, Cy, CRv);
	shading flat;
	caxis([-101 101])
	colorbar;
	colormap jet;
    set(gca, 'ydir', 'reverse');

	subplot(3,1,2);
	pcolor(Cx, Cy, CTv);
	shading flat;
	caxis([-101 101])
	colorbar;
	colormap jet;
    set(gca, 'ydir', 'reverse');

	subplot(3,1,3);
	pcolor(Cx, Cy, CVv);
	shading flat;
	caxis([-20 20])
	colorbar;
	colormap jet;
    set(gca, 'ydir', 'reverse');

	pause;
end;

% -----------------------------------------------------------------------
% Plot cylinder w radial, tangential currents mapped onto it...

for circle = 1:5,
	a = circ(circle);
	b = circ(circle+1);
	
	Cx = ones(38,b-a+1);
	Cy = ones(38,b-a+1);
	Cz = ones(38,b-a+1);
	CRv = ones(38,b-a+1);
	CVv = ones(38,b-a+1);
	CTv = ones(38,b-a+1);
	
	for j = 1:38,
		Cx(j,:) = nav_x(a:b);
		Cy(j,:) = nav_y(a:b);
		Cz(j,:) = ones(1,b-a+1).*dpth(j);
		CVv(j,:) = Vv(j, a:b);
		CRv(j,:) = EWv(j, a:b).*sin(radial(a:b))+NSv(j, a:b).*cos(radial(a:b));
		CTv(j,:) = EWv(j, a:b).*cos(radial(a:b))-NSv(j, a:b).*sin(radial(a:b));
	end;

	%---  Do separate plots of radial, tangential, and vertical velocity...	
	clf;
	surf_handle = surf(Cx, Cy, Cz);
	set(surf_handle, 'FaceColor', 'texturemap', 'cdata', flipud(CTv));
	set(gca, 'box', 'on');
	shading flat;
	caxis([-101 101])
	colorbar;
	colormap jet;

	pause;
end;


% -----------------------------------------------------------------------
% Make decimated version of arrays for plotting...

frac = 10; level = 20;
small_u = EWv(level,1:frac:length(EWv));
small_v = NSv(level,1:frac:length(EWv));
small_x = nav_x(1:frac:length(EWv));
small_y = nav_y(1:frac:length(EWv));

clf;
hold on;
color= ['r', 'g', 'b', 'c', 'm'];
for j = 1:5
	a = circ(j)/frac+1; b = circ(j+1)/frac;
	quiver(small_x(a:b), small_y(a:b), small_u(a:b), small_v(a:b), color(j));
end;
hold off;
