function f = constant_rho_left_open_right_6 (Nx, Ny, f, u_x, rho_left)

    % This function should be used in conjuction with stream_2 or similar.
    % It sets the boundary conditions on the left and right boundaries.
    % The left boundary is a constant-pressure (constant-density) boundary.
    % The right boundary is open.  For the open boundary condition, it uses
    % the CBC-local algorithm of Lou et al. (2013, Physical Review E).



    w   = [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9];


    
    i = 1;                  % left boundary
        for j = 1:Ny
            if rho_left > 0.1
                u_x(i,j) = 1 - (f(9,i,j) + f(2,i,j) + f(4,i,j) + 2*(f(3,i,j) + f(6,i,j) + f(7,i,j)))/rho_left;
                f(1,i,j) = f(3,i,j) + (2/3)*rho_left*u_x(i,j);
                f(5,i,j) = f(7,i,j) - (1/2)*(f(2,i,j) - f(4,i,j)) + (1/6)*rho_left*u_x(i,j);
                f(8,i,j) = f(6,i,j) + (1/2)*(f(2,i,j) - f(4,i,j)) + (1/6)*rho_left*u_x(i,j);
            else
                f(1,i,j) = rho_left * w(1);
                f(5,i,j) = rho_left * w(5);
                f(8,i,j) = rho_left * w(8);         
            end

        end
    
      
        
    i = Nx;                 % right boundary
        for j = 1:Ny
            lambda = u_x(i-1, j);           
            f(3,i,j) = (f(3,i,j) + lambda*f(3,i-1,j)) / (1+lambda);
            f(6,i,j) = (f(6,i,j) + lambda*f(6,i-1,j)) / (1+lambda);
            f(7,i,j) = (f(7,i,j) + lambda*f(7,i-1,j)) / (1+lambda);
        end    
    

        
end