% Convection - Diffusion - 2D; Dirichle BC on [0,1]*[0,1]
% central differences, variable coeffisients
% Lexicographical ordering, first along 'y'-direction
% The mesh   size is 'nanr=number of subintervals in one direction'
% The matrix size is 'nun=(nanr-1)^2'
%
function A=cd_cd2d(cdeps,flag,tita,hstep,nanr,nun)
%A=cpalloc(nun,nun,5*nun);
diag = 4*cdeps;
disp('Begin allocating memory...')
lengthA=nanr*nanr*5;
AI  = zeros(lengthA,1);
AJ  = zeros(lengthA,1);
AV  = zeros(lengthA,1);
nextA = 0;
disp('...end allocating memory.')%
for i=1:nanr,
for j=1:nanr,
    ind=(i-1)*nanr+j;
%    A(ind,ind) = diag; 
     nextA    = nextA + 1;
     AI(nextA) = ind;
     AJ(nextA) = ind;
     AV(nextA) = diag;
    xi=i*hstep;yj=j*hstep;
    [coefx,coefy,sgnx,sgny]=cd_coeff(xi,yj,flag,tita);
    if i-1>=1,
%       A(ind,ind-nanr)    = -(cdeps + coefx*hstep/2);
       nextA    = nextA + 1;
       AI(nextA) = ind;
       AJ(nextA) = ind-nanr;
       AV(nextA) = -(cdeps + coefx*hstep/2);
    end
    if i+1<=nanr,
%       A(ind,ind+nanr)    = -(cdeps - coefx*hstep/2);        
       nextA    = nextA + 1;
       AI(nextA) = ind;
       AJ(nextA) = ind+nanr;
       AV(nextA) = -(cdeps - coefx*hstep/2);
    end
    if j-1>=1,
%       A(ind,ind-1) = -(cdeps + coefy*hstep/2);
       nextA    = nextA + 1;
       AI(nextA) = ind;
       AJ(nextA) = ind-1;
       AV(nextA) = -(cdeps + coefy*hstep/2);
    end
    if j+1<=nanr,
%       A(ind,ind+1) = -(cdeps - coefy*hstep/2);
       nextA    = nextA + 1;
       AI(nextA) = ind;
       AJ(nextA) = ind+1;
       AV(nextA) = -(cdeps - coefy*hstep/2);
    end
end
end
AI = AI(1:nextA);
AJ = AJ(1:nextA);
AV = AV(1:nextA);

A  = sparse(AI, AJ, AV); %K0 = (K0+K0')*0.5;

  
A=A/(hstep^2);
%