function [rho, u_x, u_y] = calc_macroscopic_variables_5 ( Nx, Ny, f )
%
% The purpose of this function is to compute the fluid density (rho) and
% the fluid velocity (u_x, u_y) at each location on the lattice.  The
% lattice Boltzmann method computes these macroscopic variables from the
% distribution function f.
%
% Input variables Nx and Ny indicate the size of the domain.
% Input variable f is the LBM distribution function.
%
% This function is specific to a D2Q9 lattice.

    e_x  = [ 1,   0,  -1,   0,    1,   -1,   -1,    1,   0];
    e_y  = [ 0,   1,   0,  -1,    1,    1,   -1,   -1,   0];

    rho = zeros(Nx, Ny);
    u_x = zeros(Nx, Ny);
    u_y = zeros(Nx, Ny);
    
    for i = 1:Nx
       for j = 1:Ny
            for k = 1:9
                rho(i, j) = rho(i,j) + f(k, i, j);
                u_x(i, j) = u_x(i,j) + e_x(k) * f(k, i, j);
                u_y(i, j) = u_y(i,j) + e_y(k) * f(k, i, j);
            end
            
            u_x(i,j) = u_x(i,j) / rho(i,j);
            u_y(i,j) = u_y(i,j) / rho(i,j);
            
        end
    end

end