function [P,S,PE,PH] = findroutes(N,L,SD,H,coord)P = zeros(size(H,1),1);                 % link-path indent matrixS = zeros(size(SD,1),1);                % path-flow indent matrixlp = 0;                                 % number of pathsPE = zeros(size(H,1),size(SD,1));       % min energy (next-hop) routingPH = zeros(size(H,1),size(SD,1));       % min hop (direct) routingfor s=1:size(SD,1)    % mstack - local stack. format: i - observed link, ip - path length    mstack = [];    path = [];    is = 1;    ip = 1;    visited = zeros(N,1);    src = SD(s,1);    dst = SD(s,2);    visited(src) = 1;     minEpathcost = 1e6;    minHpathcost = 1e6;    i = 1;    while i <= size(L,1)        if dst == src            lp = lp+1;            pathEcost = 0;            for j=1:ip-1                P(path(j),lp) = 1;                pathEcost = pathEcost + 1/H(path(j),path(j));            end                        if pathEcost < minEpathcost                minEpath = path(1:ip-1);                minEpathcost = pathEcost;            end            if ip-1 < minHpathcost                minHpath = path(1:ip-1);                minHpathcost = ip-1;            end                        S(s,lp) = 1;            %mstack pop            is = is-1;            ms = mstack(is,:);            i = ms(1);            src = L(i,1);            visited(L(i,2)) = 0;            ip = ms(2);                    elseif L(i,1) == src & visited(L(i,2)) == 0 & dist(src,dst,coord) > dist(L(i,2),dst,coord) & dist(src,dst,coord) >= dist(src,L(i,2),coord)         % don't go through an already visited node, or a more distant node                    % mstack push            mstack(is,:) = [i,ip];            is = is+1;                        src = L(i,2);            path(ip) = i;            ip = ip+1;            visited(L(i,2)) = 1;            i = 0;        end        i = i+1;        while i == size(L,1) & is > 1            %mstack pop            is = is-1;            ms = mstack(is,:);            i = ms(1);            src = L(i,1);            visited(L(i,2)) = 0;            ip = ms(2);             i = i+1;        end    end        for j=1:size(minEpath,2)        PE(minEpath(j),s) = 1;    end    for j=1:size(minHpath,2)        PH(minHpath(j),s) = 1;    end        end