function f_eq = calc_f_eq ( Nx, Ny, is_fluid, u_x, u_y, rho, acceleration, tau)

%%%%    This function calculates the equilibrium distribution function f_eq
%%%%    that is used for collision or relaxation in the lattice Boltzmann
%%%%    method.  This function is particular for 2-D problems with a D2Q9
%%%%    lattice.  Also, it accounts for an external acceleration or force,
%%%%    but that force can only act in the x-direction.  The function will
%%%%    need to be modified if a force in the y-direction is also desired.

    f_eq = zeros(9, Nx, Ny);

    w  = [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9];
    e_x = [  1,   0,  -1,   0,    1,   -1,   -1,    1,   0];
    e_y = [  0,   1,   0,  -1,    1,    1,   -1,   -1,   0];

    for i = 1:Nx
        for j = 1:Ny
            if is_fluid(i,j)
                                
                u_x_eq = u_x(i,j) + tau*acceleration;
                u_y_eq = u_y(i,j);
                u_squared = u_x_eq^2 + u_y_eq^2;
                
                for k = 1:9
                    ex = e_x(k);    ey = e_y(k);
                    f_eq(k,i,j) = w(k) * rho(i,j) * (1 + ...
                        3*(ex*u_x_eq + ey*u_y_eq) + ...
                        (9/2)*(ex*u_x_eq + ey*u_y_eq)^2 - ...
                        (3/2)*u_squared);
                end
                
            end
        end
    end
    

end