function f = periodic_left_right(Nx, Ny, f, f_temp) 

%%%%    This function should be used in conjunction with the subroutine
%%%%    stream_2 or equivalent.  It streams the distribution function, 
%%%%    using f_temp to update f, for periodic boundary conditions.  
%%%%    Particles exiting on the right boundary re-enter on the left
%%%%    boundary, and vice-versa.  The function is particular to a D2Q9 
%%%%    lattice.

    in = 1;
    out = Nx;
    
    i = in;
        i_n = out;     % periodic boundary    
        for j = 1:Ny
            j_n = (j-1)*(j~=1) + (Ny)*(j==1);   % periodic boundary
            j_p = (j+1)*(j<Ny) + (1)*(j==Ny);
            f(1, i, j) = f_temp(1, i_n, j);
            f(5, i, j) = f_temp(5, i_n, j_n);
            f(8, i, j) = f_temp(8, i_n, j_p);     
        end
        
    i = out;
        i_p = in;        % periodic boundary
        for j = 1:Ny
            j_n = (j-1)*(j~=1) + (Ny)*(j==1);   % periodic boundary
            j_p = (j+1)*(j<Ny) + (1)*(j==Ny);
            f(3, i, j) = f_temp(3, i_p, j);
            f(6, i, j) = f_temp(6, i_p, j_n);
            f(7, i, j) = f_temp(7, i_p, j_p); 
        end    
        
end