% Convection - Diffusion - 2D; Dirichle BC on [0,1]*[0,1]
% upwind schemes, variable coeffisients
%
function A=cd_upwd2d(cdeps,flag,tita,hstep,nanr,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;
    xi=i*hstep;yj=j*hstep;
    [vx,vy,sx,sy]=cd_coeff(xi,yj,flag,tita);
    nextA = nextA + 1;
    AI(nextA) = ind;
    AJ(nextA) = ind;
    AV(nextA) = diag + sx*vx*hstep + sy*vy*hstep; 
    if i-1>=1,
        nextA = nextA + 1;
        AI(nextA) = ind;
        AJ(nextA) = ind-nanr;
        AV(nextA) = -cdeps;
       if sx==1,
          nextA = nextA + 1;
          AI(nextA) = ind;
          AJ(nextA) = ind-nanr;
          AV(nextA) = - vx*hstep;
       end
    end
    if i+1<=nanr,
        nextA = nextA + 1;
        AI(nextA) = ind;
        AJ(nextA) = ind+nanr;
        AV(nextA) = -cdeps;
       if sx==-1, 
          nextA = nextA + 1;
          AI(nextA) = ind;
          AJ(nextA) = ind+nanr;
          AV(nextA) = vx*hstep;
       end
    end
    if j-1>=1,
        nextA = nextA + 1;
        AI(nextA) = ind;
        AJ(nextA) = ind-1;
        AV(nextA) = -cdeps; 
       if sy==1,
          nextA = nextA + 1;
          AI(nextA) = ind;
          AJ(nextA) = ind-1;
          AV(nextA) = - vy*hstep;
       end
    end
    if j+1<=nanr,
        nextA = nextA + 1;
        AI(nextA) = ind;
        AJ(nextA) = ind+1;
        AV(nextA) = -cdeps; 
       if sy==-1,
          nextA = nextA + 1;
          AI(nextA) = ind;
          AJ(nextA) = ind+1;
          AV(nextA) = vy*hstep;
       end
    end
  end
end
AI = AI(1:nextA);
AJ = AJ(1:nextA);
AV = AV(1:nextA);

A = sparse(AI, AJ, AV); %
A = A/(hstep^2);
%
