function [u_x, f] = constant_rho_left_right(Nx, Ny, f, u_x, rho_left, rho_right)
% This function should be used in conjuction with stream_2 or similar.
% It sents the boundary conditions on the left and right boundaries.
% Both boundaries are constant-pressure (constant-density) boundaries.
% rho_left is the density on the left boundary.
% rho_right is the density on the right boundary.
% To avoid compressibility (mass balance) problems, these should not be
% very far apart from each other. I have used rho_left = 1.01 and
% rho_right = 0.99 and there was no problem with those values.
i = 1; % left boundary
for j = 1:Ny
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);
end
i = Nx; % right boundary
for j = 1:Ny
u_x(i,j) = -1 + (f(9,i,j) + f(2,i,j) + f(4,i,j) + 2*(f(1,i,j) + f(5,i,j) + f(8,i,j))) / rho_right;
f(3,i,j) = f(1,i,j) - (2/3)*rho_right*u_x(i,j);
f(6,i,j) = f(8,i,j) - (1/6)*rho_right*u_x(i,j) - (1/2)*(f(2,i,j) - f(4,i,j));
f(7,i,j) = f(5,i,j) - (1/6)*rho_right*u_x(i,j) + (1/2)*(f(2,i,j) - f(4,i,j));
end
end