function f = periodic_top_bottom_2(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 top boundary re-enter on the bottom
%%%%    boundary, and vice-versa.  The function is particular to a D2Q9 
%%%%    lattice.

    j = 1;              % bottom boundary
        j_n = Ny;       % periodic boundary
        
        i = 1;          % bottom-left corner
            i_p = i+1;
            f(2,i,j) = f_temp(2, i,   j_n);
            f(6,i,j) = f_temp(6, i_p, j_n);
     
        for i = 2:Nx-1  
            i_n = i-1;   
            i_p = i+1;
            f(2,i,j) = f_temp(2, i,   j_n);
            f(5,i,j) = f_temp(5, i_n, j_n);
            f(6,i,j) = f_temp(6, i_p, j_n);
        end
        
        i = Nx;         % bottom-right corner
            i_n = i-1;
            f(2,i,j) = f_temp(2, i,   j_n);
            f(5,i,j) = f_temp(5, i_n, j_n);
        
            
            
    j = Ny;             % top boundary
        j_p = 1;        % periodic boundary
               
        i = 1;          % top-left corner
            i_p = i+1;
            f(4,i,j) = f_temp(4, i,   j_p);
            f(7,i,j) = f_temp(7, i_p, j_p);
        
        for i = 2:Nx-1
            i_n = i-1;   
            i_p = i+1;
            f(4,i,j) = f_temp(4, i,   j_p);
            f(7,i,j) = f_temp(7, i_p, j_p);
            f(8,i,j) = f_temp(8, i_n, j_p);
        end
       
        i = Nx;         % top-right corner
            i_n = i-1;
            f(4,i,j) = f_temp(4, i,   j_p);
            f(8,i,j) = f_temp(8, i_n, j_p);        
        
        
end