% --I-- build example for final project
%==========================================
% This function creates a grid for the problem
% du/dx + dv/dy = f
% du/dy - dv/dx = g
%===========================================
% Remark the matrix represents the unit squre in the XY plane where
% increase in rows corresponds to Y-change and columns for X-change
%   |
%   |
%  -|----------->
%=====================================================
% Y - first dimension  | X - second dimension
%=====================================================
full_weighting_mask = .25* [ .25 .5 .25; .5 1 .5; .25 .5 .25];
%==============
res = input('Case or Default? (1,..,5 / def - 0)  ');
if (res > 0)
    n = 2^7;
else    
    n = input('n = ');
end
h = 1 / n;
X = [0:h:1];
Y = [0:h:1];
[X,Y] = meshgrid(X,Y);
%-=============================
switch res
    case 1
        U = Y + sin(X +Y);
        V = X + X.*sin(Y);
       'U = Y + sin(X +Y)'
       'V = X + X.*sin(Y)'
       U_name= 'U = Y + sin(X +Y)';
       V_name='V = X + X.*sin(Y)';
    case 2
        U = exp(X).*cos(Y);
        V = exp(X).*sin(Y);
       'U = exp(X).*cos(Y)'
       'V = exp(X).*sin(Y)'
       U_name='U = exp(X).*cos(Y)';
       V_name='V = exp(X).*sin(Y)';
   case 3
        U = sin(3*X);
        V = 4*Y + X;
       'U = sin(3*X)'
       'V = 4*Y + X'
       U_name='U = sin(3*X)';
       V_name='V = 4*Y + X';
    case 4
        U = sin(5*X + Y);
        V = exp(X + 1);
       'U = sin(5*X + Y)'
       'V = exp(X + 1)' 
       U_name='U = sin(5*X + Y)';
       V_name='V = exp(X + 1)' ;
    case 5
        U = exp(sin(X + Y));
        V = cos( X.*Y + 2);
       'U = exp(sin(X + Y))'
       'V = cos( X.*Y + 2)'
       U_name='U = exp(sin(X + Y))';
       V_name='V = cos( X.*Y + 2)';
   otherwise
        U = input('U(X,Y) = ');
        V = input('V(X,Y) = ');
        U_name= 'user-defined';
        V_name= 'user-defined';

end

%U_grid = imfilter(U,full_weighting_mask,'same','replicate');
%V_grid = imfilter(V,full_weighting_mask,'same','replicate');
%==============================
%subplot(2,2,1);mesh(X,Y,U);title('U')
%subplot(2,2,2);mesh(X,Y,V);title('V')

s = size(U); rows = s(1) ; cols = s(2);

F = zeros(s);
G = zeros(s);


%==================================
% Build F & G
%==================================
for i= 2:2:rows - 1
    for j = 2:2:cols - 1
        F(i,j) = (U(i,j+1) - U(i,j-1))/(2*h) + (V(i+1,j) - V(i-1,j))/(2*h);
    end
end

for i= 3:2:rows - 2
    for j = 3:2:cols - 2
        G(i,j) = (U(i+1,j) - U(i-1,j))/(2*h) - (V(i,j+1) - V(i,j-1))/(2*h) ;
    end
end

%==================================
% Build Guess Matrix
%==================================
%Guess = rand(s);
flag = 1;
% real Guess
%===========
if (flag == 1) 
    Guess = rand(s);%zeros(s);
    
    Guess = F + G;
    
    Guess(1,:) = V(1,:);
    Guess(rows,:) = V(rows,:);
    
    Guess(:,1) = U(:,1);
    Guess(:,cols) = U(:,cols);
    
    Guess(1,1) = 0; Guess(1,cols) = 0; Guess(rows,1) = 0;Guess(rows,cols) = 0;
else
    % trivial
    %========
    Guess = zeros(s);
    Guess = F+G;
    Guess(2:2:rows-1,1:2:cols) = U(2:2:rows-1,1:2:cols);
    Guess(1:2:rows,2:2:cols-1) = V(1:2:rows,2:2:cols-1);
end