function x=ngrid(z,T,alpha,tol,omega,nrange)
%
% NGRID  x=ngrid(z,T,alpha,tol,omega,nrange)
%
%    Input a grid (z) populated by data cells, cells to be filled (flagged
%    by Inf), and cells to be masked (flagged by -Inf).
%
%    Parameters:
%	omega		SOR parameter (experience shows that to be safe,
%                                      keep this thing near 1. (~1.05)
%	T		Tension:  T=0 is pure Bi-Harmonic. 
%				  T=1 is pure Laplacian.  Try T=0.35
%	alpha		anisotropy of grid: alpha = dy/dx
%	tol		tolerance for convergence test.  Stop if maximum
%				change at any grid element is less than tol
%	nrange		number of neighboring points influenced by data
%

%    Implements the GMT gridding algorithm 
%    See "Gridding with Continuous Splines in Tension" by
%    Smith and Wessel, 1990 in  Geophysics, 55(3), P. 293-305
%    The data is gridded by solving 
%
%         (1 - T) * L (L (z)) - T * L (z) = 0    at the cells to be filled,
%
%    where T is a tension factor between 0 and 1, and L indicates the Lapla-
%    cian operator.  T = 0 gives the "minimum curvature" solution which is
%    equivalent to SuperMISP and the ISM packages.  Minimum curvature can 
%    cause undesired oscillations and false local maxima or minima (See Smith
%    and Wessel manuscript), and you may wish to use T > 0 to suppress these
%    effects.  Experience suggests T ~ 0.25 usually looks good for potential
%    field data and T should be larger (T ~ 0.35) for steep topography data.
%    T = 1 gives a harmonic surface (no maxima or minima are possible except
%    at control data points). 
%

% Derek Fong and Rich Signell received (zgrid.m) 3/24/95
% Modified ClH 5/10/95 for T<> 1
% Modified and bug fixes by DAF 9/27/95

%..........................................................................
% if Nargin==0, demo routine

if(nargin==0), 
   omega=1.4;                          	% SOR parameter 
   T=.5;      				% tension                   
   alpha=1.0; 				% anisotropy of grid:  alpha = dy/dx
   tol=.1;    				% tolerance for convergence test
   nrange=5;                            % neighboring influence range
%
% create data from "peaks"
%
   sp=30;     				% size of peaks data set for gridding
   z=peaks(sp);
   z2=z;
   [m,n]=size(z);
%
% flag cells to be filled by gridding algorithm
%
%   rand('seed', 0);
   ifill=ceil(rand(sp*sp*.90,1)*(m*n-1));
   z(ifill)=Inf*ones(size(ifill));
%
% flag cells to be masked, i.e. not filled by gridding algorithm
%
   im=[ceil(m/3):1:floor(m/2)];
   jm=[ceil(n/3):1:floor(n/2)];
   in=length(im);
   jn=length(jm);
   z(im,jm)=-Inf*ones(length(im),length(jm));
   z2(im,jm)=NaN*ones(length(im),length(jm));

   x=ngrid(z,T,alpha,tol,omega,nrange);        % do gridding


   aa=find(~finite(x));                   %black out non-working points
   x(aa)=nan*ones(size(aa));
   bb=find(~finite(z));
   z(bb)=nan*ones(size(bb));
   set(gcf,'pos',[10 10 800 300]);

   subplot(1,3,2)                               % show data
   pcolor(z);shading('flat')
   set(gca,'dataaspectratio',[1 alpha 1])
%  show "land region"
   line([im(1) im(in)+1 im(in)+1 im(1) im(1)],...
     [jm(1) jm(1) jm(jn)+1 jm(jn)+1 jm(1)],'color','red')
   text(.5*(im(1)+im(in)),.5*(jm(1)+jm(jn)),'land','horiz','center');
   title('Sampled Data')

   subplot(1,3,3)                               % show gridded results
   pcolor(x);shading('flat')
   set(gca,'dataaspectratio',[1 alpha 1])
%  show "land region"
   line([im(1) im(in)+1 im(in)+1 im(1) im(1)],...
      [jm(1) jm(1) jm(jn)+1 jm(jn)+1 jm(1)],'color','red')
   text(.5*(im(1)+im(in)),.5*(jm(1)+jm(jn)),'land','horiz','center');
   title('Gridded Data')
   
   subplot(1,3,1)                               % show "true" field
   pcolor(z2);shading('flat')
   set(gca,'dataaspectratio',[1 alpha 1])
%  show "land region"
   line([im(1) im(in)+1 im(in)+1 im(1) im(1)],...
      [jm(1) jm(1) jm(jn)+1 jm(jn)+1 jm(1)],'color','red')
   text(.5*(im(1)+im(in)),.5*(jm(1)+jm(jn)),'land','horiz','center');
   title('True Field')
   return
end

% end demo
%..........................................................................
%..........................................................................
% Initialize constants.

a2=alpha^2;
a4=alpha^4;

%..........................................................................
% Add two layers of land cells on all sides for boundary conditions.

[m,n]=size(z);
m4=m+4;
n4=n+4;
z2=-Inf*ones(m4,n4);
z2(3:m4-2,3:n4-2)=z;

%..........................................................................
% Create a matrix y where the neighboring influence has been set.
% Set idata as the index of the original data points (before spreading).
%     ifind as the index of points which need gridding (interpolation).
%     imask as land points to mask out.
y = z2;
y = spread(y,nrange);

idata = find(~isinf(z2));                      % original data
ifind = find(z2==Inf & y~=-Inf);               % points to grid
imask = find(y==-Inf);                         % "land" points
ntot=m4*n4;                                    % dimensions of y



% .........................................................................
% Determine on which points the full biharmonic equation is to be applied,
% and upon which we will just use the laplacian.
% When we finish, "lfind" will be the index of points to use the laplacian
% on, and "bfind" will be the index of points to use the full biharmonic.
bfind = [];
lfind = [];
for i = 1:length(ifind)

  biharmonic_influence = [ ...
      ifind(i)-2; ...
      ifind(i)-m4-1; ...
      ifind(i)-1; ...
      ifind(i)+m4-1; ...
      ifind(i)-2*m4; ...
      ifind(i)-m4; ...
      ifind(i)+m4; ...
      ifind(i)+2*m4; ...
      ifind(i)-m4+1; ...
      ifind(i)+1; ...
      ifind(i)+m4+1; ...
      ifind(i)+2; ...
      ];
  
  % if any grid points within the biharmonic range of influence are -inf,
  %then we must list them to use the laplacian only.
  if ( ~isempty((find(z2(biharmonic_influence) == -inf))) | ...
       ~isempty((find(y(biharmonic_influence) == -inf))) )
    lfind = [lfind; ifind(i)];
  else
    bfind = [bfind; ifind(i)];
  end
  
end


bn = length(bfind);
ln = length(lfind);
			

%..........................................................................
% Compute coefficient array.
%
% Biharmonic
%                      --                                      --
%         c1          |                   a4                     |
%     c2  c3  c4      |      2*a2     -4*(1+a2)*a2     2*a2      |
% c5  c6  c7  c8  c9  | 1  -4*(1+a2) -(6+8*a2+6*a4) -4*(1+a2)  1 |*(1-T)  -
%    c10  c11 c12     |      2*a2     -4*(1+a2)*a2     2*a2      |
%         c13         |                   a4                     |
%                      --                                      --
% Laplacian            --                                      --
%         c3          |                   a2                     |
%      c6 c7 c8       |        1        2*(1+a2)      1          |*T = 0
%         c11         |                   a2                     |
%                      --                                      -- 
%  where a2=alpha^2,   a4=alpha^4, and alpha = dy/dx
%
%  from Smith & Wessel Eqn (A-4)
%

cB = zeros(1,13);
cB([1 13])        = a4*([1 1]);
cB([2  4 10 12])  = 2*a2*([1 1 1 1]);
cB([6  8])        = -4*(1+a2)*([1 1]);    
cB([3 11])        = -4*(1+a2)*a2*([1 1]);
cB([5 9])         =[1 1];
cB(7)             = -(6 + 8*a2 + 6*a4);          % sum of weights

cL = zeros(1, 13);
cL([3 11]) = a2*[1 1];
cL([6 8]) = [1 1];
cL(7) = 2*(1 + a2);                              % sum of weights

c = (1-T)*cB - T*cL;

biharmonic = (1-T)*cB - T*cL;
biharmonic(7) = [];

laplacian = [a2 1 1 a2]';


%..........................................................................
% Successive overrelaxation to find gridded matrix.

max_residual=Inf;

while (max_residual > tol)

  % do biharmonic part
  for i = 1:bn
  
    ynew = dot ( biharmonic, ... 
                 [ ...
			 y(bfind(i) -2) ...
			 y(bfind(i)-m4-1) ...
			 y(bfind(i)-1) ...
			 y(bfind(i)+m4-1) ...
			 y(bfind(i)-2*m4) ...
			 y(bfind(i)-m4) ...
			 y(bfind(i)+m4) ...
			 y(bfind(i)+2*m4) ...
			 y(bfind(i)-m4+1) ...
			 y(bfind(i)+1) ...
			 y(bfind(i)+m4+1) ...
			 y(bfind(i)+2) ...
	         ] ) / c(7);
    residual(i)=y(bfind(i))-ynew;
    y(bfind(i))=(1-omega)*y(bfind(i)) + omega*ynew;

  end
  
  
 
  % do laplacian part
  for i=1:ln

  
    % grab 3x3 square on which to apply laplacian
    lsq = [ ...
	        y(lfind(i) - 1) ...
		y(lfind(i)-m4) ...
		y(lfind(i)+m4) ...
		y(lfind(i) + 1) ...
	  ]';
  
    % If any border cells are -Inf, set them equal to the center point.
    % This has the effect of making the normal derivative zero there.
    landfind = find (lsq == -inf);
    lsq(landfind) = y(lfind(i)) * ones(size(landfind));
    
    ynew = dot ( laplacian, lsq ) / cL(7);

    residual(bn+i)=y(lfind(i))-ynew;
    y(lfind(i))=(1-omega)*y(lfind(i)) + omega*ynew;
  end

  max_residual=max(abs(residual));

  
  disp(sprintf('Residual error: %e\n',max_residual)) 
  

end


%..........................................................................
% Here's the final gridded result.

y(imask) = -Inf*ones(size(imask));            % set zeros back to "land"
x=y(3:m4-2,3:n4-2);
disp('done gridding!');

%..........................................................................
%..........................................................................

