function [xg, yg, zg] = potgrid(x, y, z, size, weight)

% Function potential = potgrid(x, y, z, size, weight)
%    x = array of x coordinates
%    y = array of y coordinates
%    z = array of z coordinates
%    size = [nx, ny] = size of desired potential 
%    weight = [wi, wj] = weighting of 

% Get size of desired output array...
max_i = size(2);
max_j = size(1);

% Get weights for Laplacian solution...
wi = weight(1);
wj = weight(2);

% Eliminate NaNs to simplify processing...
ok = find((~isnan(x) & ~isnan(y)) & ~isnan(z));
x = x(ok);
y = y(ok);
z = z(ok);

% -----------------------------------------------------------------------

x1 = min(min(x)); x2 = max(max(x)); 
xs = (x2-x1)/(max_j-2);
xi = (x1-xs/2):xs:(x2+xs/2);

y1 = min(min(y)); y2 = max(max(y));
ys = (y2-y1)/(max_i-2);
yi = (y1-ys/2):ys:(y2+ys/2);

[xi, yi] = meshgrid(xi', yi);

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

pf = zeros(max_i,max_j);
nn = zeros(max_i,max_j);

% Compute location of data in grid...
ii = round((max_j-2)*(x - x1)/(x2 - x1))+1;
jj = round((max_i-2)*(y - y1)/(y2 - y1))+1;
ij = (ii-1)*max_i + jj;

% Order data...
% [ij, ok] = sort(ij);
% z = z(ok);

for i = 1:length(ij), 
	pf(ij(i)) = pf(ij(i)) + z(i);
	nn(ij(i)) = nn(ij(i)) + 1;
end;

pf(find(nn)) = pf(find(nn))./nn(find(nn));
nn(find(nn)) = nn(find(nn))./nn(find(nn));
memory = pf(find(nn));
clear ii; clear jj; clear ij; clear ok;

% pcolor(xi,yi,pf); 
% set( gca, 'YDir', 'rev'); 
% caxis([2.8 3.5]); 
% shading flat; colorbar

% -----------------------------------------------------------------------
% Initialize pf assuming little horizontal variability...
%   make sure there are no NaNs in the initial array....

a = mean(pf'.*nn')./mean(nn');
a1 = min(find(~isnan(a)));
if a1 > 1
   a(1:(a1-1)) = a(a1);
end;

a2 = max(find(~isnan(a)));
if a2 < (length(a)+1)
   a((a2+1):length(a)) = a(a2);
end;

ind = 1:max_i; a = interp1(ind(find(~isnan(a))), a(find(~isnan(a))), ind, 'nearest');

b = ones(max_j,1);
pf = a' * b';
pf(find(nn))=memory;

clear a; clear b; clear a1; clear a2; clear ind;

% -----------------------------------------------------------------------
% Iteratively solve Laplaces equation for the potential defined above...

for i = 1:1000,
	% Compute field in center of region...
	pfxp = pf(3:max_i, 2:max_j-1);
	pfxm = pf(1:max_i-2, 2:max_j-1);
	pfyp = pf(2:max_i-1,3:max_j);
	pfym = pf(2:max_i-1,1:max_j-2);
	pf(2:max_i-1,2:max_j-1) = (wi*(pfxp + pfxm) + wj*(pfyp + pfym))/(2*wi+2*wj);

	% Handle edges...
	pf(1,2:max_j-1) = (wi*pf(2,2:max_j-1)+wj*pf(1,1:max_j-2)+wj*pf(1,3:max_j))/(wi+2*wj);	
	pf(max_i,2:max_j-1) = (wi*pf(max_i-1,2:max_j-1)+wj*pf(max_i,1:max_j-2)+wj*pf(max_i,3:max_j))/(wi+2*wj);
	pf(2:max_i-1,1) = (wj*pf(2:max_i-1,2)+wi*pf(1:max_i-2,1)+wi*pf(3:max_i,1))/(wj+2*wi);
	pf(2:max_i-1,max_j) = (wj*pf(2:max_i-1,max_j-1)+wi*pf(1:max_i-2,max_j)+wi*pf(3:max_i,max_j))/(wj+2*wi);
	
	% Handle corners...
	pf(1,1) = (wj*pf(1,2)+wi*pf(2,1))/(wi+wj);
	pf(max_i,1) = (wj*pf(max_i,2)+wi*pf(max_i-1,1))/(wi+wj);
   pf(1,max_j) = (wj*pf(1,max_j-1)+wi*pf(2,max_j))/(wi+wj);
	pf(max_i,max_j) = (wj*pf(max_i,max_j-1)+wi*pf(max_i-1,max_j))/(wi+wj);

	% Reset initial potential...
	pf(find(nn))=memory;
end;

% -----------------------------------------------------------------------
% Make sure we only plot for regions with good data...
% Compute mask...
mask = (nn > 0);
memory = mask(find(nn));
for i = 1:50,
	% Compute field in center of region...
	maskxp = mask(3:max_i, 2:max_j-1);
	maskxm = mask(1:max_i-2, 2:max_j-1);
	maskyp = mask(2:max_i-1,3:max_j);
	maskym = mask(2:max_i-1,1:max_j-2);
	mask(2:max_i-1,2:max_j-1) = (wi*(maskxp + maskxm) + wj*(maskyp + maskym))/(2*wi+2*wj);

	% Handle edges...
	mask(1,2:max_j-1) = (wi*mask(2,2:max_j-1)+wj*mask(1,1:max_j-2)+wj*mask(1,3:max_j))/(wi+2*wj);	
	mask(max_i,2:max_j-1) = (wi*mask(max_i-1,2:max_j-1)+wj*mask(max_i,1:max_j-2)+wj*mask(max_i,3:max_j))/(wi+2*wj);
	mask(2:max_i-1,1) = (wj*mask(2:max_i-1,2)+wi*mask(1:max_i-2,1)+wi*mask(3:max_i,1))/(wj+2*wi);
	mask(2:max_i-1,max_j) = (wj*mask(2:max_i-1,max_j-1)+wi*mask(1:max_i-2,max_j)+wi*mask(3:max_i,max_j))/(wj+2*wi);
	
	% Handle corners...
	mask(1,1) = (wj*mask(1,2)+wi*mask(2,1))/(wi+wj);
	mask(max_i,1) = (wj*mask(max_i,2)+wi*mask(max_i-1,1))/(wi+wj);
   mask(1,max_j) = (wj*mask(1,max_j-1)+wi*mask(2,max_j))/(wi+wj);
	mask(max_i,max_j) = (wj*mask(max_i,max_j-1)+wi*mask(max_i-1,max_j))/(wi+wj);

	% Reset initial potential...
	mask(find(nn))=memory;
end;
not_ok = find(mask < 0.1);
pf(not_ok) = NaN;

% figure; pcolor(pf); shading flat; colormap jet; colorbar;
% hold on;
% plot((max_j-1)*(x - mnx)/(mxx - mnx)+1, (max_i-1)*(y - mny)/(mxy - mny)+1, 'k');

% figure;

pcolor(xi,yi,pf); 
set( gca, 'YDir', 'rev'); 
% caxis([2.8 3.5]);
shading flat;
colorbar;
hold on
plot(x, y, 'w');
% ylabel( 'depth (m)' );
% xlabel( 'x (m)' );
hold off;

xg = x;
yg = y;
zg = pf;

return;