function Nr=CreateReachNetwork(N,P)
% Create reach network from pixel network

maxdegreein=max(N.degreein);
nodes_r=find(N.degreein==0);  %leaves of network
for d=2:maxdegreein
    nodes_r=[nodes_r; find(N.degreein==d)];
end
nodes_r(nodes_r==N.outlet)=[];
nodes_r=[nodes_r; N.outlet]; %added outlet at the end of the array

%Nr structure with data of the reach network
Nr.X=N.X(nodes_r);
Nr.Y=N.Y(nodes_r);
Nr.n_nodes=length(nodes_r);
Nr.AD=(sparse(Nr.n_nodes,Nr.n_nodes,0));
Nr.L=zeros(Nr.n_nodes,1);
Nr.pixel2reach=zeros(N.n_nodes,1); %for each pixel, it indicates the ID of the reach it belongs to
Nr.A=zeros(Nr.n_nodes,1); 
for int=1:Nr.n_nodes-1             %It done not search the connection of the outlet (the last element of vector nodescond)
    int;
    LL=0;
    nodo=nodes_r(int);
    Nr.pixel2reach(nodo)=int;
    drainage_A_last=N.A(nodo);
    
    for passo=1:N.n_nodes
        nodosucc=N.down_node(nodo);
        LL=LL+N.L(nodo);
        if ismember(nodosucc,nodes_r)
            Nr.L(int)=LL;
            Nr.A(int)=drainage_A_last;
            Nr.AD(int,(nodes_r==nodosucc))=1;
            break
        end
        nodo=nodosucc;
        Nr.pixel2reach(nodo)=int;
        drainage_A_last=N.A(nodo);
    end 
end
Nr.A(end)=N.A(N.outlet);
Nr.pixel2reach(N.outlet)=Nr.n_nodes;
Nr.DDA=zeros(Nr.n_nodes,1);
for  int=1:Nr.n_nodes
    Nr.DDA(int)=Nr.A(int)-sum(Nr.A(full(Nr.AD(:,int))==1));
end
Nr.outlet=find(sum(Nr.AD,2)==0);

%In the reach network the outlet is a fake node with zero length, thus it
%removed
Nr.AD=Nr.AD(1:end-1,1:end-1);  %Because the outlet is the last node for construction
DDA_outlet=Nr.DDA(Nr.outlet);
Nr.DDA(Nr.outlet)=[];
Nr.A(Nr.outlet)=[];
Nr.X(Nr.outlet)=[];
Nr.Y(Nr.outlet)=[];
Nr.L(Nr.outlet)=[];

Nr.n_nodes=length(Nr.A);
old_outlet_r=Nr.outlet;
Nr.outlet=find(sum(Nr.AD,2)==0);
Nr.DDA(Nr.outlet)=Nr.DDA(Nr.outlet)+DDA_outlet; 
Nr.A(Nr.outlet)=Nr.A(Nr.outlet)+DDA_outlet;
Nr.pixel2reach(Nr.pixel2reach==old_outlet_r)=Nr.outlet;
clear DDA_outlet old_outlet_r

% Create Nr.down_node (down_node of each reach)
Nr.down_node=zeros(Nr.n_nodes,1);
for i=1:Nr.n_nodes
    if not(i==Nr.outlet)
        Nr.down_node(i)=find(Nr.AD(i,:));
    end
end

% Compute slope for each reach
Nr.slope=P.slope_outlet*(Nr.A/Nr.A(Nr.outlet)).^P.alpha_slope;

% Create CATCH matrix, whose element ij is 1 if j is in the subcatchent closed at i
Nr.CATCH=(sparse(Nr.n_nodes,Nr.n_nodes,0));
for j=1:Nr.n_nodes
    Nr.CATCH(j,j)=1;
    node=j;
    for k=1:Nr.n_nodes
        if Nr.down_node(node)>0
            node=Nr.down_node(node);
            Nr.CATCH(node,j)=1;
        else
            break
        end
    end
end

