function x=ogrid(z,T,alpha,tol,omega,nrange)
%
% DGRID  x=dgrid(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);
   [m,n]=size(z);
%
% flag cells to be filled by gridding algorithm
%
%  rand('seed',0);
   ifill=ceil(rand(sp*sp*.8,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));
%   x=spgrid(z,T,alpha,tol,omega,nrange);        % do gridding
   x=ogrid(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));

   figure;
   subplot(1,2,1)                               % show data
   pcolor(z)
   set(gca,'aspectratio',[alpha alpha])
   hold on
   plot([im(1) im(in)+1 im(in)+1 im(1) im(1)],...
     [jm(1) jm(1) jm(jn)+1 jm(jn)+1 jm(1)])
   title('Irregular Data')
   hold off

   subplot(1,2,2)                               % show gridded results
   pcolor(x)
   set(gca,'aspectratio',[alpha alpha])
   hold on
   plot([im(1) im(in)+1 im(in)+1 im(1) im(1)],...
      [jm(1) jm(1) jm(jn)+1 jm(jn)+1 jm(1)])
   title('Gridded Data')
   hold off
   figure(1);
   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 the index of points where gridding shouldn't be done ("land").

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
nfind = length(ifind);                         % number of points to grid

ntot=m4*n4;                                    % dimensions of y

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

%..........................................................................
% Create masks for each grid point where interpolation will occur.  Vectors 
% are 1 if point will be interpolated, 0 if masked as "land" (=-Inf).
%
% NB! this is not the most mathematically rigorous way of handling differ-
% encing due to missing "influence" cells; however, it is a quick and dirty
% method that seems to work so long as we keep omega small (may not make
% SOR most efficient) so that problem is stable.
%

msk=zeros(nfind,13);

msk(:, 1) = (y(ifind - 2 ) ~= -Inf);         % c1
msk(:, 2) = (y(ifind-m4-1) ~= -Inf);         % c2
msk(:, 3) = (y(ifind - 1 ) ~= -Inf);         % c3
msk(:, 4) = (y(ifind+m4-1) ~= -Inf);         % c4         
msk(:, 5) = (y(ifind-2*m4) ~= -Inf);         % c5
msk(:, 6) = (y(ifind - m4) ~= -Inf);         % c6
msk(:, 8) = (y(ifind + m4) ~= -Inf);         % c8
msk(:, 9) = (y(ifind+2*m4) ~= -Inf);         % c9
msk(:,10) = (y(ifind-m4+1) ~= -Inf);         % c10
msk(:,11) = (y(ifind + 1 ) ~= -Inf);         % c11
msk(:,12) = (y(ifind+m4+1) ~= -Inf);         % c12
msk(:,13) = (y(ifind + 2 ) ~= -Inf);         % c13

%..........................................................................
% Multiply mask vectors by coeffients. 

d=msk.*((ones(nfind,1)*c));

d(:,1) = d(:,1).*msk(:,3);     % if c3 is masked, mask c1
d(:,5) = d(:,5).*msk(:,6);     % if c6 is masked, mask c5
d(:,9) = d(:,9).*msk(:,8);     % if c8 is masked, mask c9
d(:,13) = d(:,13).*msk(:,11);  % if c11 is masked, mask c13


% If there is no mask, then d(:,7) is all set.  But if a mask was used then
% the sum of weights is wrong and needs to be adjusted to reflect the
% actual weights used.

d(:,7) = sum(d')';

%..........................................................................
% Make data values finite in masked regions instead of -Inf so that when
% they are multiplied by zero, you get zero!

y(imask)=ones(size(imask));



% have plot window open for updates at each step
%loop_fig = figure;
%subplot(1,2,1);
%pcolor(z2);
%set(gca,'aspectratio',[alpha alpha])
%subplot(1,2,2);
%pcolor(y);
%set(gca,'aspectratio',[alpha alpha])

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

max_residual=Inf;
res1=Inf;
res2=Inf;
res3=Inf;
while (max_residual > tol),
  for i=1:nfind;
    ynew = (d(i,1).*y(ifind(i) -2) + ...
            d(i,2).*y(ifind(i)-m4-1) + ...
            d(i,3).*y(ifind(i)-1) + ...
	    d(i,4).*y(ifind(i)+m4-1) +...
            d(i,5).*y(ifind(i)-2*m4) + ...
	    d(i,6).*y(ifind(i)-m4) + ...
            d(i,8).*y(ifind(i)+m4) + ...
            d(i,9).*y(ifind(i)+2*m4) + ...
   	    d(i,10).*y(ifind(i)-m4+1) + ...
            d(i,11).*y(ifind(i)+1) + ...
	    d(i,12).*y(ifind(i)+m4+1) + ...
            d(i,13).*y(ifind(i)+2))./d(i,7);
    residual(i)=y(ifind(i))-ynew;
    y(ifind(i))=(1-omega)*y(ifind(i)) + omega*ynew;
  end
%  res3=res2;
%  res2=res1;
%  res1=max_residual;
  max_residual=max(abs(residual));
%  if (max_residual>res1 & max_residual>res2 & max_residual>res3)
%     max_residual=0;
%  end;

  % at each step, let's see how we're doing
%  subplot(1,2,2);
%  pcolor(y);
%  set(gca,'aspectratio',[alpha alpha])
  disp(max_residual);
%  pause;


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!');

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

