function [grid, L, Phi_0]  = hewitt_grid_v2(n_rows,n_cols,L_x,L_y,visualize)
%hewitt_grid(n_rows,n_cols,L_x,L_y)
%   Generate trapezoidal grid oriented, with opposite corners of trapezes
%   oriented down / across-slope, and alternating across and down-slope
%   diagonals included as network edges similar to Hewitt (2013, EPSL)
%   Input arguments:
%       n_rows:     number of rows in grid
%       n_columns:  number of columns in grid
%       L_x:        spacing between rows
%       L_y:        spacing between columns
%       visualize:  flag set to true if the code should plot the network;
%                   this results in a legible plot only if n_rows and
%                   n_columns is small
%   Output:
%       grid: structure with basic network information
%               n_nodes: number of nodes in network
%               n_edges: number of network edges
%               up_node: list of upstream nodes for each edge
%               down_node: list of downstream nodes for each edge
%                   (up_node and down_node together give the data expected
%                   in a connectivity array)
%               x,y: coordinates of grid nodes
%       L:      network edge lengths, used in the computation of gradients
%               along network edges
%       Phi_0:  A default background hydraulic gradient, defined as Phi_0 =
%               -x at each node
%To visualize the network, 


%check input
if n_rows/2 == round(n_rows/2)
    error('n_rows must be odd')
end

%initialize output
grid.n_nodes = n_rows*(n_cols+1/2) - 1/2;
grid.n_edges = (n_rows-1)*(3*n_cols-1/2)-n_cols+1;
%(n_rows-1)*(3*n_cols-1/2)-n_cols+2
grid.x = zeros(grid.n_nodes,1);
grid.y = zeros(grid.n_nodes,1);
grid.up_node = zeros(grid.n_edges,1);
grid.down_node = zeros(grid.n_edges,1);
grid.bdy_nodes = zeros(n_cols,1);

%work by rows of nodes, except for last two rows
for ii_row=1:n_rows
    %odd row number
    if ii_row/2 ~= round(ii_row/2)
        jj_col = 1:n_cols;
        ii_node = (ii_row-1)*(n_cols+1/2)+jj_col;
        grid.x(ii_node) = (ii_row-1)*L_x;
        grid.y(ii_node) = (2*jj_col-1)*L_y;
        %last row has no edges attached
        if ii_row < n_rows
            for kk_edge = 1:3
                ii_edge = (ii_row-1)*(3*n_cols-1/2)+3*(jj_col-1)+kk_edge;
                grid.up_node(ii_edge) = ii_node;
                switch kk_edge
                    case 1
                        grid.down_node(ii_edge) = ii_node+2*n_cols+1;
                    case 2
                        grid.down_node(ii_edge) = ii_node+n_cols;
                    case 3
                        grid.down_node(ii_edge) = ii_node+n_cols+1;
                end
            end
        else
           grid.bdy_nodes = ii_node.'; 
        end
    %even row numbers    
    else
        jj_col = 1:n_cols+1;
        ii_node = ii_row*(n_cols+1/2)-n_cols-1+jj_col;
        grid.x(ii_node) = (ii_row-1)*L_x;
        grid.y(ii_node) = 2*(jj_col-1)*L_y;
        if ii_row < n_rows-1
            for kk_edge = 1:3
                switch kk_edge
                    case 1
                        %asymmetry of left-most node
                        jj_col1 = jj_col(2:end);
                        ii_node1 = ii_node(2:end);
                        ii_edge = ii_row*(3*n_cols-1/2)-3*n_cols+1+3*(jj_col1-1)+kk_edge-2;
                        grid.up_node(ii_edge) = ii_node1;
                        grid.down_node(ii_edge) = ii_node1+n_cols;
                    case 2
                        jj_col2 = jj_col(2:end-1);
                        ii_node2 = ii_node(2:end-1);
                        ii_edge = ii_row*(3*n_cols-1/2)-3*n_cols+1+3*(jj_col2-1)+kk_edge-2;
                        grid.up_node(ii_edge) = ii_node2+n_cols;
                        grid.down_node(ii_edge) = ii_node2+n_cols+1;
                    case 3 
                        %asymmetry of right-most node
                        jj_col3 = jj_col(1:end-1);
                        ii_node3 = ii_node(1:end-1);
                        ii_edge = ii_row*(3*n_cols-1/2)-3*n_cols+1+3*(jj_col3-1)+kk_edge-2;
                        grid.up_node(ii_edge) = ii_node3;
                        grid.down_node(ii_edge) = ii_node3+n_cols+1;
                end
            end
        else
            for kk_edge = 1:2
                switch kk_edge
                    case 1
                        %asymmetry of left-most node
                        jj_col1 = jj_col(2:end);
                        ii_node1 = ii_node(2:end);
                        ii_edge = ii_row*(3*n_cols-1/2)-3*n_cols+1+2*(jj_col1-1)+kk_edge-1;
                        grid.up_node(ii_edge) = ii_node1;
                        grid.down_node(ii_edge) = ii_node1+n_cols;
                    case 2 
                        %asymmetry of right-most node
                        jj_col2 = jj_col(1:end-1);
                        ii_node2 = ii_node(1:end-1);
                        ii_edge = ii_row*(3*n_cols-1/2)-3*n_cols+1+2*(jj_col2-1)+kk_edge-1;
                        grid.up_node(ii_edge) = ii_node2;
                        grid.down_node(ii_edge) = ii_node2+n_cols+1;    
                end
            end
                        
        end
    end
end

%edge lengths and potentials
L = sqrt((grid.x(grid.up_node)-grid.x(grid.down_node)).^2 + (grid.y(grid.up_node)-grid.y(grid.down_node)).^2);
Phi_0 = -grid.x;

%visualization
if visualize
    %load geometry
    X = grid.x;
    Y = grid.y;
    up_node = grid.up_node;
    down_node = grid.down_node;

    %generate plots
    figure, hold on
    h=plotpanel(1,plotparams);
    hold on
    plot([X(up_node), X(down_node)]',[Y(up_node), Y(down_node)]','k');
    plot(X,Y,'ro','MarkerSize',5)
    for ii=1:length(X)
        text(X(ii),Y(ii),num2str(ii),'FontSize',16)
    end
    xlim([min(X) max(X)])
    ylim([min(Y) max(Y)])
    box on
end

end