%% 2D Convection-diffusion problem
% Parameters:
%            n        -> number of intervals in each direction
%        nun=(n-1)^2) -> total number of unknowns   
%  1 grad = 0.017453 rad

clear A b residual x sol

rpi = 0.017453;
n   = input('Refinement parameter: ');
nanr= n-1;
nun = nanr^2;
h   = 1/n;
cdeps=input('Epsilon: ');
disp('Choose a vector field:')
disp('flag = 0 - no convection	     ')
disp('flag = 1 - rotating wind	     ')
disp('flag = 2 - nonlinear convection')
disp('flag = 3 - horizontal wind     ')
disp('flag = 4 - vertical wind	     ')

flag =input('Vector field flag: ');
if flag==1,
   tita=30*rpi; 
else
  tita=1;
end
%
A_cd = cd_cd2d(cdeps,flag,tita,h,nanr,nun);

A_up = cd_upwd2d(cdeps,flag,tita,h,nanr,nun);
sol  = cd_sol(h,nanr);
b_cd = A_cd*sol;
b_up = A_up*sol;
 
disp('Matrices and rhs are ready.')

%% Solution central differences ___
% % disp('Central differences')
%  [x_cdfo,r_cdfo,k_cdfo] = GaussSeidel(A_cd,b_cd,'fo');
% %wait
% [x_cdba,r_cdba,k_cdba] = GaussSeidel(A_cd,b_cd,'ba');
% 
% figure(1),clf
% plot(sol),hold
% plot(x_cdfo,'r'),plot(x_cdba,'g')
 %% Solution with upwind discretization ___
 [x_upfo,r_upfo,k_upfo] = GaussSeidel(A_up,b_up,'fo');
 %wait
 [x_upba,r_upba,k_upba] = GaussSeidel(A_up,b_up,'ba');
 
 figure(2),clf
 plot(sol),hold
 plot(x_upfo,'r'),plot(x_upba,'g')
