function [xp, yp, pf] = potgridy(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 
% In contrast to potgrid, potgridy initializes by interpolating
% along the x axis (assumes y variability uncoupled).

% 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); 

% Generate a mesh grid....
x1 = min(x); x2 = max(x); xs = (x2-x1)/max_j; xi = (x1-xs):xs:(x2+xs);
y1 = min(y); y2 = max(y); ys = (y2-y1)/max_i; yi = (y1-ys):ys:(y2+ys);
[xi, yi] = meshgrid(xi', yi);

% grid for plotting (JR)
xs = (x2-x1)/(max_j-2); xp = (x1-xs/2):xs:(x2+xs/2);
y1 = min(min(y)); y2 = max(max(y));
ys = (y2-y1)/(max_i-2); yp = (y1-ys/2):ys:(y2+ys/2);
[xp, yp] = meshgrid(xp', yp);


% 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...
mxx = max(max(x)); mnx = min(min(x)); mxy = max(max(y)); mny = min(min(y));
ii = round((max_j-1)*(x - mnx)/(mxx - mnx))+1;
jj = round((max_i-1)*(y - mny)/(mxy - mny))+1;
ij = (ii-1)*max_i + jj;

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;

% ------------------------------------------------------
% Interpolate on a row by row basis....

%for i = 1:max_i,
%   ind = find(nn(i,:));
%   if (length(ind) > 1)
%      pffill(i,:) = interp1(xi(i,ind), pf(i, ind), xi(i,:));
%   elseif (length(ind) == 1)
%      pffill(i,:) = ones(1, max_j) * pf(i, ind);
%   else
%      pffill(i,:) = pffill(i-1,:);
%   end;
%end;
%pf = pffill;

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

for i = 1:500,
   % 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;

