% assign boundary conditions
outlet_ic = 1;
inlet_ic = 0;

count = 0;
tolerence = 1e-4;
alpha = 1.825;
TotErrSq = 0;
converge = 0;

% set initial conditions
for i=1:15
    for j=1:15
        phi(i,j) = 0;
    end
end
for i=1:9
    phi(i,15) = outlet_ic;
end

phi2 = phi;



while(~converge)
    
    TotErrSq = 0;

    for j=14:-1:1
        for i=14:-1:1
    
            % make sure the point is in the domain of interest
            if( ~((j>5) && (i>9)) )
    
                % top left corner
                if((i==1) && (j==1))
                    phi2(1,j) = 0.25*(2*phi(i+1,j)+2*phi(i,j+1));
                               
                % top wall
                elseif(j==1)
                    phi2(i,j) = 0.25*( phi(i+1,j)+phi(i-1,j)+2*phi(i,j+1) );
                
                % left wall
                elseif(i==1)
                    phi2(i,j) = 0.25*(phi(i,j+1)+phi(i,j-1)+2*phi(i+1,j));

                % right wall
                elseif((i==9) && (j>5))
                    phi2(i,j)= 0.25*(phi(i,j+1)+phi(i,j-1)+2*phi(i-1,j));

                % bottom wall (right)
                elseif((j==5) && (i>9))
                    phi2(i,j) = 0.25*(phi(i+1,j)+phi(i-1,j)+2*phi(i,j-1) );

                % all other points
                else
                    phi2(i,j) = 0.25*(phi(i-1,j) + phi(i+1,j) + phi(i,j-1) + phi(i,j+1));

                end
                
                % calculate new phi
                phi2(i,j)=(1-alpha)*phi(i,j) + alpha*phi2(i,j);
               
                % sum the errors
                err = (phi2(i,j) - phi(i,j))^2;
                TotErrSq = TotErrSq + err;
                
                phi(i,j) = phi2(i,j);
            end
            
        end
    end
        
                
    count = count + 1;
    
    % check for convergence
    if(TotErrSq < tolerence)
        converge = 1;
    end

end

count
TotErrSq
contourf(phi);
colorbar;


% Calculate the volume averaged velcoity at inlet

i=14;
for j=1:5
    ui(j) = phi(i,j) - phi(i+1,j);
end
U_inlet = (0.5*ui(1)+ui(2)+ui(3)+ui(4)+0.5*ui(5))/4;

% Calculate velocity (divided by U_inlet) at point B
uB = (phi(14,1) - phi(15,1)) / U_inlet;

% calculate the u (x-direction) component of velocity between nodes (center
% of cell) along the walls

% top wall
j=1;
for i=1:14
    u(i,j) = (phi(i,j) - phi(i+1,j)) / U_inlet;
    cp(i,j) = uB^2 - u(i,j)^2;
end


% bottom, right wall
j=5;
for i=10:14
    u(i,j) = (phi(i,j) - phi(i+1,j)) / U_inlet;
    cp(i,j) = uB^2 - u(i,j)^2;
end

% left wall
i=1;
for j=1:14
    v(i,j) = (phi(i,j+1) - phi(i,j)) / U_inlet;
    cp(i,j) = uB^2 - v(i,j)^2;
end

% right wall
i=9;
for j=6:14
    v(i,j) = (phi(i,j+1) - phi(i,j)) / U_inlet;
    cp(i,j) = uB^2 - v(i,j)^2;
end



% print all phi values
for j=1:15
    phi(:,j)
end

% print top wall 
u(:,1)
cp(:,1)

%print bottom wall 
u(10:14,5)
cp(10:14,5)

%print left wall
v(1,:)
cp(1,:)

%print right wall
v(9,6:14)
cp(9,6:14)