function zout=linADI(K, alpha, beta, gammafun, deltafun,...
   InitCondfun, Bdryfuns, Rect, N, M, T, L, disp)
% linADI  Solves a linear PDE
%          u_t = K(u_xx + u_yy) + alpha u_x + beta u_y 
%                               + gamma u + delta
%       on a rectangular region.
%    The constant coefficients K, alpha, beta may be complex.  
%    VERSION I simply always uses gamma=delta=0 and ignors 
%       gammafun and deltafun.
%    Domain is a  rectangle described by "Rect".  
%    Evaluates name (string), or inline object, in InitCondfun
%       to determine initial condition in grid coordinates.
%    Uses the names (strings), or inline objects, in Bdryfuns
%       to evaluate Dirichlet boundary data.  Boundary functions
%       must be time independent.  Boundary functions
%       are specified in the following order:
%           lower, right, upper, left.
%    Grid is N subdivisions in x, M in y. 
%    Computes L steps from t=0 to t=T.  Displays every 
%       "disp"th step.
%    See "exampleADI.m".
% (Ed Bueler; 4/1/02)

clear global NN MM ax ay bx by dx dy BF
global NN MM ax ay bx by dx dy BF

gamma=0.0;delta=0.0; % KLUDGE FOR NOw
ax=Rect(1); bx=Rect(2); ay=Rect(3); by=Rect(4);
dx=(bx-ax)/N; dy=(by-ay)/M; dxdx=dx*dx; dydy=dy*dy; dt=T/L;
Q=(N-1)*(M-1); % size of matrices and vectors
NN=N; MM=M; BF=Bdryfuns; % global variables seem to need this

% global <--> local index matrices: 
%   r(j,k), resp. s(j,k), are global names of j,k grid points; 
%   if z is Q by 1 vector in r order then z(r) is in grid coords;
%   if z is Q by 1 vector in s order then z(s) is in grid coords;
%   if w is in grid coords then w(:) is in r order.
r=reshape(1:Q,N-1,M-1);  s=reshape(1:Q,M-1,N-1)';
% build permutation vector: 
%    if z is Q by 1 vector in r order then z(prs) is in s order
%    if z is Q by 1 vector in s order then z(psr) is in r order
rt=r';  prs=rt(:)';  psr=s(:)';
% to see perm matrices: I=eye(Q); P=I(p,:), P*P'

%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct equations in r order:  Ar new = Br old + cr.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
Ar=spdiags((2/dt+2*K/dxdx)*ones(Q,1),0,Q,Q);
Br=spdiags((2/dt-2*K/dydy+gamma)*ones(Q,1),0,Q,Q);
cr=zeros(Q,1);
aLr=-K/dxdx+alpha/(2*dx);  aRr=-K/dxdx-alpha/(2*dx);
bDr=K/dydy-beta/(2*dy);  bUr=K/dydy+beta/(2*dy);
%%%% k=1 case
% j=1 case
rr=r(1,1);
cr(rr,1)=bDr*bdryval(1,0)+delta-aLr*bdryval(0,1);
Ar(rr,r(2,1))=aRr;
Br(rr,r(1,2))=bUr;
% general j case
for j=2:N-2
   rr=r(j,1);
   cr(rr,1)=bDr*bdryval(j,0)+delta;
   Ar(rr,r(j-1,1))=aLr;  Ar(rr,r(j+1,1))=aRr;
   Br(rr,r(j,2))=bUr;
end;
% j=N-1 case
rr=r(N-1,1);
cr(rr,1)=bDr*bdryval(N-1,0)+delta-aRr*bdryval(N,1);
Ar(rr,r(N-2,1))=aLr;
Br(rr,r(N-1,2))=bUr;
%%%% general k case 
for k=2:M-2
   % j=1 case
   rr=r(1,k);
   cr(rr,1)=delta-aLr*bdryval(0,k);
   Ar(rr,r(2,k))=aRr;
   Br(rr,r(1,k-1))=bDr;   Br(rr,r(1,k+1))=bUr;
   % general j case
   for j=2:N-2
      rr=r(j,k);
      cr(rr,1)=delta;
      Ar(rr,r(j-1,k))=aLr;  Ar(rr,r(j+1,k))=aRr;
      Br(rr,r(j,k-1))=bDr;  Br(rr,r(j,k+1))=bUr;
   end;
   % j=N-1 case
   rr=r(N-1,k);
   cr(rr,1)=delta-aRr*bdryval(N,k);
   Ar(rr,r(N-2,k))=aLr;
   Br(rr,r(N-1,k-1))=bDr;  Br(rr,r(N-1,k+1))=bUr;
end;
%%%% k=M-1 case
% j=1 case
rr=r(1,M-1);
cr(rr,1)=bUr*bdryval(1,M)+delta-aLr*bdryval(0,M-1);
Ar(rr,r(2,M-1))=aRr;
Br(rr,r(1,M-2))=bDr;
% general j case
for j=2:N-2
   rr=r(j,M-1);
   cr(rr,1)=bUr*bdryval(j,M)+delta;
   Ar(rr,r(j-1,M-1))=aLr;  Ar(rr,r(j+1,M-1))=aRr;
   Br(rr,r(j,M-2))=bDr;
end;
% j=N-1 case
rr=r(N-1,M-1);
cr(rr,1)=bUr*bdryval(N-1,M)+delta-aRr*bdryval(N,M-1);
Ar(rr,r(N-2,M-1))=aLr;
Br(rr,r(N-1,M-2))=bDr;


%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Construct equations in s order:  As wnew = Bs wold + cs.
%%%%%%%%%%%%%%%%%%%%%%%%%%%
As=spdiags((2/dt+2*K/dydy-gamma)*ones(Q,1),0,Q,Q);
Bs=spdiags((2/dt-2*K/dxdx)*ones(Q,1),0,Q,Q);
cs=zeros(Q,1);
aDs=-K/dydy+beta/(2*dy);  aUs=-K/dydy-beta/(2*dy);
bLs=K/dxdx-alpha/(2*dx);  bRs=K/dxdx+alpha/(2*dx);
%%%% k=1 case
% j=1 case
ss=s(1,1);
cs(ss,1)=-aDs*bdryval(1,0)+delta+bLs*bdryval(0,1);
As(ss,s(1,2))=aUs;
Bs(ss,s(2,1))=bRs;
% general j case
for j=2:N-2
   ss=s(j,1);
   cs(ss,1)=-aDs*bdryval(j,0)+delta;
   As(ss,s(j,2))=aUs;
   Bs(ss,s(j-1,1))=bLs;  Bs(ss,s(j+1,1))=bRs;
end;
% j=N-1 case
ss=s(N-1,1);
cs(ss,1)=-aDs*bdryval(N-1,0)+delta+bRs*bdryval(N,1);
As(ss,s(N-1,2))=aUs;
Bs(ss,s(N-2,1))=bLs;
%%%% general k case 
for k=2:M-2
   % j=1 case
   ss=s(1,k);
   cs(ss,1)=delta+bLs*bdryval(0,k);
   As(ss,s(1,k-1))=aDs;  As(ss,s(1,k+1))=aUs;
   Bs(ss,s(2,k))=bRs;
   % general j case
   for j=2:N-2
      ss=s(j,k);
      cs(ss,1)=delta;
      As(ss,s(j,k-1))=aDs;  As(ss,s(j,k+1))=aUs;
      Bs(ss,s(j-1,k))=bLs;  Bs(ss,s(j+1,k))=bRs;
   end;
   % j=N-1 case
   ss=s(N-1,k);
   cs(ss,1)=delta+bRs*bdryval(N,k);
   As(ss,s(N-1,k-1))=aDs;  As(ss,s(N-1,k+1))=aUs;
   Bs(ss,s(N-2,k))=bLs;
end;
%%%% k=M-1 case
% j=1 case
ss=s(1,M-1);
cs(ss,1)=bLs*bdryval(0,M-1)+delta-aUs*bdryval(1,M);
As(ss,s(1,M-2))=aDs;
Bs(ss,s(2,M-1))=bRs;
% general j case
for j=2:N-2
   ss=s(j,M-1);
   cs(ss,1)=-aUs*bdryval(j,M)+delta;
   As(ss,s(j,M-2))=aDs; 
   Bs(ss,s(j-1,M-1))=bLs;  Bs(ss,s(j+1,M-1))=bRs;
end;
% j=N-1 case
ss=s(N-1,M-1);
cs(ss,1)=-aUs*bdryval(N-1,M)+delta+bRs*bdryval(N,M-1);
As(ss,s(N-1,M-2))=aDs;
Bs(ss,s(N-2,M-1))=bLs;

% to see the matrices: full(Ar), full(Br), full(As), full(Bs)
% to see cr, cs: cr, cs
% note   wic=ic(:)  is in r order, wic(prs) is in s order

x=linspace(ax+dx,bx-dx,N-1); % interior points
y=linspace(ay+dy,by-dy,M-1);
[X,Y]=meshgrid(x,y); X=X'; Y=Y'; % this way X(j,k)=x_jk etc.

zin=feval(InitCondfun,X,Y);
figure(1);
clf reset;
axis([ax bx ay by 0 1]);
hold on;
mp = mesh(x,y,zin');
drawnow 
xlabel('x axis'); ylabel('y axis'); zlabel('u axis');
title('Initial condition.  Hit any key to continue.');
axis auto
pause

%time loop: solve two tridiagonal systems at each time step
zr=zin(:); %r order
for ll=1:L
   tzr=Ar\(Br*zr+cr);
   tzs=tzr(prs);
   zs=As\(Bs*tzs+cs);
   zr=zs(psr);
   zout=zr(r);
   if mod(ll,disp)==0
      set(mp,'ZData',zout');
      set(mp,'CData',zout');
      drawnow
      title(['Time t=' num2str(ll*dt,'%5.5f')]);
      axis auto
   end;
end;
hold off;

% returns boundary vals in grid coords
function u=bdryval(j,k);

global NN MM ax ay bx by dx dy BF

if j==0  % left boundary
   temp=feval(BF,ax,ay+k*dy); u=temp(4);
elseif j==NN  % right boundary
   temp=feval(BF,bx,ay+k*dy); u=temp(2);
elseif k==0  % lower boundary
   temp=feval(BF,ax+j*dx,ay); u=temp(1);
elseif k==MM  % upper boundary
   temp=feval(BF,ax+j*dx,by); u=temp(3);
else error('cannot evaluate bdry at nonbdry point'); end;