function grains_constant_rho_left_right (grain_arrangement_file, Nx, Ny, nu, rho_left, rho_right, t_max)

% grains_constant_rho_left_right.m
%
% Created by Jeff Cunningham while on sabbatical from USF and visiting 
% CSIRO in Floreat, Western Australia, 2013-2014
%
% Uses lattice Boltzmann modeling to simulate 2D fluid flow through a
% porous medium.  The user specifies the geometry of the porous medium by
% providing the file grain_arrangement_file.
%
% Flow is left to right.
% This code differes from its predecessor, grains.m, in that this code does
%   not assume periodic boundary conditions on the left and right boundaries.
%   Instead, it uses specified density (pressure)on the left and right
%   boundaries.
%
% Some parts of this lattice Boltzmann code are modeled after codes posted
% on-line by Jonas Latt and co-workers.  I also have received assistance
% from Shadab Anwar and Mike Sukop as I have developed this code.  
%
% grain_arrangement_file is the name of a file containing data that specify
%   the geometry of the porous domain.  A 1 indicates that the node is
%   occupied by a grain; a 0 indicates that the node is fluid.
% Nx and Ny are the size of the domain in lattice units.  Nx and Ny should
%   agree with the data in grain_arrangement_file.
% nu is the fluid viscosity in LBM units (lu^2/ts).  
%   Recommend nu = 1/6.
% rho_left and rho_right are the density on the left and right boundaries,
%   respectively.  Recommend both close to 1.0.  For instance, rho_left =
%   1.005 and rho_right = 0.995.
% t_max is the number of time steps for the simulation to run.
%   t_max must be high enough that the system reaches steady state.
%   As a rough guideline, use t_max = Nx*Ny/10 .




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Geometry of the problem    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
fid = fopen(grain_arrangement_file, 'r');
pore_wall = fscanf(fid, '%i', [Nx, Ny]);
fclose(fid); 
is_fluid = 1 - pore_wall;       % if it's not a pore wall, then it is fluid
[y, x] = meshgrid(1:Ny, 1:Nx);      

figure(1);
%imagesc(pore_wall');
surf(pore_wall');   shading flat;
axis('equal');
axis([1 Nx 1 Ny]);
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Map of the porous medium', 'FontSize', 16);
drawnow;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Set some LBM variables    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
average_density = 0.5*(rho_left + rho_right);

omega  = 1. / (3*nu + 1./2.);   % relaxation parameter
tau = 1/omega;

green_light = true;                % parameter to tell if we should proceed


    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    D2Q9 lattice constants    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% Use this numbering scheme for the directions:
% 6 2 5
% 3 9 1
% 7 4 8
%
w  =  [1/9, 1/9, 1/9, 1/9, 1/36, 1/36, 1/36, 1/36, 4/9];
opp = [  3,   4,   1,   2,    7,    8,    5,    6,   9];
%e_x = [ 1,   0,  -1,   0,    1,   -1,   -1,    1,   0];
%e_y = [ 0,   1,   0,  -1,    1,    1,   -1,   -1,   0];




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Pre-allocate the f matrices                       %%%%
%%%%    This makes Matlab happy                           %%%%
%%%%    Some of the matrices don't need pre-allocation    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
f = zeros(9, Nx, Ny);
f_temp = zeros(9, Nx, Ny);
%f_eq = zeros(9, Nx, Ny);   
%u_x = zeros(Nx, Ny);
%u_y = zeros(Nx, Ny);
%rho = zeros(Nx, Ny);

 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Initialize the distribution function f    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for k = 1:9
    f(k, 1:Nx, 1:Ny) = w(k) * average_density;
end
initial_fluid_mass = sum(sum(   is_fluid .* reshape(sum(f),Nx,Ny)   ));    % Track the fluid mass over time.
disp(initial_fluid_mass);                                                   % Verifies code is running properly.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Begin main loop    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t = 0;

while (green_light) && (t < t_max)

    t = t+1;
       
   
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Calculate macroscopic variables    %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [rho, u_x, u_y] = calc_macroscopic_variables_5 ( Nx, Ny, f );

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Every 50 time steps, check the mass balance    %%%%
    %%%%    Also display the velocity field                %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    if (mod(t, 50) == 0)
        green_light = check_mass_balance(t, initial_fluid_mass, rho, is_fluid);
        u = reshape(sqrt(u_x.^2 + u_y.^2), Nx, Ny);  
        figure(2);  
        %imagesc(u');    
        %surf(u');   colorbar;    
        pcolor(u');
        axis('equal');    axis([1 Nx 1 Ny]);    shading interp;
        drawnow;        
    end  
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%    
    %%%%    Calculate equilibrium distribution function     %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Use accel = 0, i.e., no external force applied in this simulation.
    f_eq = calc_f_eq ( Nx, Ny, is_fluid, u_x, u_y, rho, 0, tau);
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Collision or relaxation    %%%%
    %%%%    Use standard BGK method    %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    for i = 1:Nx
        for j = 1:Ny
            if is_fluid(i,j)    % fluid nodes relax to equilibrium
                for k = 1:9
                    f_temp(k,i,j) = f(k,i,j) - (f(k,i,j) - f_eq(k,i,j))/tau;
                end
            else                % solid nodes use bounce-back
                for k = 1:9
                    f_temp(k,i,j) = f(opp(k),i,j); 
                end
            end
        end
    end
    

    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%%%    Streaming step                                               %%%%
    %%%%    Periodic boundaries on top and bottom boundaries.            %%%%
    %%%%    Constant density (pressure) on left and right boundaries.    %%%%
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    f = stream_2(Nx, Ny, f, f_temp);
    f = periodic_top_bottom_2(Nx, Ny, f, f_temp);
    [u_x, f] = constant_rho_left_right(Nx, Ny, f, u_x, rho_left, rho_right);
    % I think it is important to put periodic_top_bottom_2 before
    % constant_rho_left_right in this case.  There are 32 values of f that
    % are not determined by stream_2.  Fourteen are covered by
    % periodic_top_bottom_2, and the other eighteen are covered by
    % constant_rho_left_right.  However, some of the calculations in
    % constant_rho_left_right depend upon knowing certain values of f.
    % Therefore I think we should update as much of f as possible before
    % running constant_rho_left_right.

end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Compute streamlines    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
u_x = u_x';                         % Need to transpose these to get the 
u_y = u_y';                         % figure to turn out properly

% I have hard-wired the starting locations for the streamlines.  That might
% not be a great way to do this, but I have found that it seems to work
% pretty well for most possible grain configurations.
y_start = 3:6:Ny;
vec_length = length(y_start);
x_start = ones(vec_length, 1);
x_start = [x_start', 6:6:Nx, 6:6:Nx];
new_length = length(6:6:Nx);
y_start = [y_start, ones(1,new_length), (Ny-1)*ones(1,new_length)];
XY = stream2(x', y', u_x, u_y, x_start, y_start);




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Graph streamlines    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(3);
orient(3, 'portrait');
set(3, 'PaperPosition', [1 2.5 6.5 6]);
h = streamline(XY);
set(h, 'Color', 'blue');
set(h, 'LineWidth', 1.);
axis('equal');
axis([1 Nx 1 Ny]);
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Streamlines of flow through a porous medium', 'FontSize', 16);
hold on;
contour(x, y, pore_wall);
hold off;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    Graph velocity vectors    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(4);
orient(4, 'portrait');
contour(x, y, pore_wall);
hold on;
h = quiver(x(4:4:Nx,4:4:Ny), y(4:4:Nx,4:4:Ny), (u_x(4:4:Nx,4:4:Ny))', (u_y(4:4:Nx,4:4:Ny))' );
set(h, 'Color', 'blue');
xlabel('x coordinate (lattice units)', 'FontSize', 16);
ylabel('y coordinate (lattice units)', 'FontSize', 16);
title('Visualization of velocity vectors', 'FontSize', 16);
axis('equal');
axis([1 Nx 1 Ny]);
hold off;




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%    End the program    %%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end