function Eout = particleclaw_network_solver(par,V,E)
%PARTICLECLAW_NETWORK_SOLVER
%   Particle method for scalar hyperbolic conservation laws
%   laws on networks.
%
%   Description of method:
%   https://link.springer.com/chapter/10.1007/978-3-642-32979-1_13
%   https://arxiv.org/abs/1201.0057
%
% 12/2010 by Benjamin Seibold and Yossi Farjoun
%            http://www.math.temple.edu/~seibold/


%========================================================================
% Preprocessing
%========================================================================
% compute lengths of edges
for j = 1:size(E,1)
    E{j,2} = norm(V{E{j,1}(2),3}-V{E{j,1}(1),3});
end
% extend/truncate solution to edges
for j = 1:size(E,1)
    L = E{j,2}; x = E{j,4}; u = E{j,5};
    [x,u] = reduce_to_edge(x,u,L);
    E{j,4} = x; E{j,5} = u;
end
% make extremal values consistent
for j = 1:size(V,1)
    E = generalized_Riemann_problem(V,E,j);
end

%========================================================================
% Time Evolution
%========================================================================
par.t = 0;
plot_network(par,V,E)
while par.t<par.tfinal
    dt = min(par.dt,par.tfinal-par.t);
    par.t = par.t+dt;
    E = evolve_solution_step(par,V,E,dt);
    plot_network(par,V,E)
end
if nargout==1, Eout = E(:,1:end-1); end



%========================================================================
% Computation Functions
%========================================================================
function newE = evolve_solution_step(par,V,E,dt)
% Perform one update step of the full network, that is
% (1) Evolve edges for a time step dt independently
% (2) Synchronize edges by solving a generalized Riemann problem at
%     each node and re-distributing area near each node
nv = size(V,1); ne = size(E,1);
for j = 1:ne               % loop over edges (evolve edges independently)
    x = E{j,4}; u = E{j,5}; rhoM = E{j,3}{1}; vM = E{j,3}{2};
    [x,u,area_credit] = evolve_edge(par,x,u,rhoM,vM,dt);
    E{j,6} = area_credit;   % extend edge data structure locally for area
    E{j,4} = x; E{j,5} = u;
end
for j = 1:nv                           % loop over vertices (synchronize)
    E = synchronize(V,E,j,dt);
end
newE = E;

%------------------------------------------------------------------------

function [newx,newu,area_credit] = evolve_edge(par,x,u,rhoM,vM,t)
% Evolve solution on a single edge for a time increment t.
% Particles will not remain coinciding with domain of edge [0,L]
% First and last particle need special treatment when involved in a
% merge.
[x,u] = coarsen_and_refine(par,x,u);                % particle management
f1 = vM*(1-(2/rhoM)*u);                       % characteristic velocities
area_credit = [0 0];                      % area credit on left and right
while t>0 % repeat until time to advance is zero
    ts = -diff(x)./diff(f1); ts(ts<0) = inf;
    tmax = min(ts);                                   % maximal time step
    dt = min(tmax,t);
    x = x+dt*f1;                      % move particles characteristically
    t = t-dt;
    while 1 % repeat while particles remain to be merged
        dx = diff(x); du = diff(u);
        ind = find(dx<1e-8&dx./(abs(du)+eps)<1);
        if isempty(ind), break, end
        ind = ind(1); % resolve the first merge that is found
        % resolve shocks at ends by extrapolating solution
        if ind==length(x)-1 % last particle is involved
            x = [x x(end)+par.d]; u = u([1:end end]);
            area_credit(2) = area_credit(2)+par.d*u(end); % area on right
        end
        if ind==1 % first particle is involved
            x = [x(1)-par.d x]; u = u([1 1:end]); ind = 2;
            area_credit(1) = area_credit(1)+par.d*u(1);    % area on left
        end
        % potentially insert particles near shock
        if x(ind)-x(ind-1)>par.d*1.5            % if enough space, insert
            [x,u] = insert(x,u,x(ind)-par.d);    % particle left of shock
            ind = ind+1;
        end
        if x(ind+2)-x(ind+1)>par.d*1.5          % if enough space, insert
            [x,u] = insert(x,u,x(ind+1)+par.d); % particle right of shock
        end
        % perform the particle merge
        [xn,un] = merge(x(ind+(-1:2)),u(ind+(-1:2)));
        x = [x(1:ind-1) xn x(ind+2:end)];
        u = [u(1:ind-1) un u(ind+2:end)];
    end
    f1 = vM*(1-(2/rhoM)*u);                   % characteristic velocities
end
newx = x; newu = u;


function [xnew,unew] = merge(x,u)
% Merge 4 points into 3, with the middle two replaced by a single one,
% preserving area. Outputs the one new particle only.
A = sum(diff(x).*(u(1:3)+u(2:4)));            % doubled area before merge
xnew = (x(2)+x(3))/2;              % merged point in middle of old points
unew = (A-(xnew-x(1))*u(1)-(x(4)-xnew)*u(4))/(x(4)-x(1));


function [xnew,unew] = insert(x,u,xn)
% For particles (x,u), insert a new particle at xn. Outputs all
% particles.
ind = sum(x<xn); % index of particle directly left to xn
if ind==0             % new particle left outside, constant extrapolation
    xnew = [xn n]; unew = u([1 1:end]);
elseif ind==length(x)% new particle right outside, constant extrapolation
    xnew = [x xn]; unew = u([1:end end]);
else       % new particle between two existing ones, linear interpolation
    un = u(ind)+(u(ind+1)-u(ind))*(xn-x(ind))/(x(ind+1)-x(ind));
    xnew = [x(1:ind) xn x(ind+1:end)];
    unew = [u(1:ind) un u(ind+1:end)];
end


function [xnew,unew] = coarsen_and_refine(par,x,u)
% Coarsen particles on an edge by removing redundant particles,
% and inserts new particles into too large gaps.
while 1
    d2x = x(3:end)-x(1:end-2); 
    u_mid = u(1:end-2)+(u(3:end)-u(1:end-2)).*...        % values between
        (x(2:end-1)-x(1:end-2))./d2x;           % particles of distance 2
    error = abs(u(2:end-1)-u_mid);
    error = error+(d2x>par.h*2);
    [mv,mi] = min(error);
    if mv>1e-9, break, end                             % error threshold
    x = x([1:mi mi+2:end]); u = u([1:mi mi+2:end]);
end
while 1
    dx = diff(x);
    [mv,mi] = max(dx);
    if mv<par.h*2, break, end       % minimum distance between particles
    [x,u] = insert(x,u,(x(mi)+x(mi+1))/2);
end
xnew = x; unew = u;


function [xnew,unew] = sharpen(x,u)
% Sharpen shocks into true discontinuities.




%------------------------------------------------------------------------

function newE = synchronize(V,E,j_node,t)
% For a network node, synchronize all adjancent edges:
% - define an interpolant on the intervals [0,L_i] again
% - solve generalized Riemann problems for extremal states
% - arrange excess and virtual areas
edge_in  = V{j_node,1}; n_in  = length(edge_in);
edge_out = V{j_node,2}; n_out = length(edge_out);
if n_in==1&&n_out==1                                         % bottleneck
    rhoM1 = E{edge_in {1},3}{1}; vM1 = E{edge_in {1},3}{2};
    rhoM2 = E{edge_out{1},3}{1}; vM2 = E{edge_out{1},3}{2};
    L = E{edge_in {1},2};                             % length of in-edge
    x_in  = E{edge_in {1},4}; u_in  = E{edge_in {1},5};   % states before
    x_out = E{edge_out{1},4}; u_out = E{edge_out{1},5};      % operations
    % truncate/extrapolate solution to edges
    [xn_in ,un_in ,I_in ] = fit_to_domain_right(x_in ,u_in ,L);
    [xn_out,un_out,I_out] = fit_to_domain_left(x_out,u_out);
    % solve Riemann problem
    [ub_in,ub_out] = Riemann_problem_11(...
        rhoM1,vM1,rhoM2,vM2,un_in(end),un_out(1));
    [xn_in,un_in] = add_jump(-xn_in(end:-1:1),un_in(end:-1:1),ub_in);
    xn_in = -xn_in(end:-1:1); un_in = un_in(end:-1:1);        % add jumps
    [xn_out,un_out] = add_jump(xn_out,un_out,ub_out);
    % manage areas
    C = u_out(1)*dflux(rhoM2,vM2,u_out(1))-...     % constant that equals
        u_in(end)*dflux(rhoM1,vM1,u_in(end));             % (d/dt)(I2-I1)
    I_add = I_out-E{edge_out{1},6}(1)+I_in-E{edge_in {1},6}(2)+C*t;
    if dflux(rhoM1,vM1,ub_in)>0           % information goes to the right
        [xn_out,un_out] = match_area_by_function(xn_out,un_out,I_add);
    else                                   % information goes to the left
        [xn_in,un_in] = match_area_by_function(...   % apply to rightmost
            -xn_in(end:-1:1),un_in(end:-1:1),I_add);      % particles, by
        xn_in = -xn_in(end:-1:1); un_in = un_in(end:-1:1);   % reflection
    end
    E{edge_in {1},4} = xn_in;  E{edge_in {1},5} = un_in;
    E{edge_out{1},4} = xn_out; E{edge_out{1},5} = un_out;
elseif n_in==1&&n_out==2                                    % bifurcation
    rhoM1    = E{edge_in {1},3}{1}; vM1    = E{edge_in {1},3}{2};
    rhoM2{1} = E{edge_out{1},3}{1}; vM2{1} = E{edge_out{1},3}{2};
    rhoM2{2} = E{edge_out{2},3}{1}; vM2{2} = E{edge_out{2},3}{2};
    L = E{edge_in {1},2};                             % length of in-edge
    A{1} = V{j_node,4}; A{2} = 1-A{1};                 % direction matrix
    x_in     = E{edge_in {1},4}; u_in     = E{edge_in {1},5};    % states
    x_out{1} = E{edge_out{1},4}; u_out{1} = E{edge_out{1},5};    % before
    x_out{2} = E{edge_out{2},4}; u_out{2} = E{edge_out{2},5};% operations
    % truncate/extrapolate solution to edges
    [xn_in ,un_in ,I_in ] = fit_to_domain_right(x_in ,u_in ,L);
    [xn_out{1},un_out{1},I_out{1}] = fit_to_domain_left(x_out{1},u_out{1});
    [xn_out{2},un_out{2},I_out{2}] = fit_to_domain_left(x_out{2},u_out{2});
    % solve Riemann problem
    [ub_in,ub_out{1},ub_out{2}] = Riemann_problem_12(...
        rhoM1,vM1,rhoM2{1},vM2{1},rhoM2{2},vM2{2},A{1},...
        un_in(end),un_out{1}(1),un_out{2}(1));
    [xn_in,un_in] = add_jump(-xn_in(end:-1:1),un_in(end:-1:1),ub_in);
    xn_in = -xn_in(end:-1:1); un_in = un_in(end:-1:1);        % add jumps
    [xn_out{1},un_out{1}] = add_jump(xn_out{1},un_out{1},ub_out{1});
    [xn_out{2},un_out{2}] = add_jump(xn_out{2},un_out{2},ub_out{2});
    % manage areas
    for j = 1:2
        C = u_out{j}(1)*dflux(rhoM2{j},vM2{j},u_out{j}(1))-...
            u_in(end)*dflux(rhoM1,vM1,u_in(end))*A{j};
        I_add{j} = I_out{1}-E{edge_out{j},6}(1)+...
            (I_in-E{edge_in {1},6}(2))*A{j}+C*t;
    end
    if dflux(rhoM2{1},vM2{1},ub_out{1})<=0    % out-edge 2 is influencing
        [xn_in,un_in] = match_area_by_function(...   % apply to rightmost
            -xn_in(end:-1:1),un_in(end:-1:1),I_add{1}/A{1});  % particles
        xn_in = -xn_in(end:-1:1); un_in = un_in(end:-1:1);% by reflection
        [xn_out{2},un_out{2}] = match_area_by_function(...
            xn_out{2},un_out{2},I_add{1}/A{1}*A{2});
    elseif dflux(rhoM2{2},vM2{2},ub_out{2})<=0 %out-edge 3 is influencing
        [xn_in,un_in] = match_area_by_function(...   % apply to rightmost
            -xn_in(end:-1:1),un_in(end:-1:1),I_add{2}/A{2});  % particles
        xn_in = -xn_in(end:-1:1); un_in = un_in(end:-1:1);% by reflection
        [xn_out{1},un_out{1}] = match_area_by_function(...
            xn_out{1},un_out{1},I_add{2}/A{2}*A{1});
    else                                % information enters from in-edge
        for j = 1:2
            [xn_out{j},un_out{j}] = match_area_by_function(...
                xn_out{j},un_out{j},I_add{j});
        end
    end
    % assign results to edge data structure
    E{edge_in {1},4} = xn_in;     E{edge_in {1},5} = un_in;
    E{edge_out{1},4} = xn_out{1}; E{edge_out{1},5} = un_out{1};
    E{edge_out{2},4} = xn_out{2}; E{edge_out{2},5} = un_out{2};
elseif n_in==2&&n_out==1                                     % confluence
    rhoM1{1} = E{edge_in {1},3}{1}; vM1{1} = E{edge_in {1},3}{2};
    rhoM1{2} = E{edge_in {2},3}{1}; vM1{2} = E{edge_in {2},3}{2};
    rhoM2    = E{edge_out{1},3}{1}; vM2    = E{edge_out{1},3}{2};
    L{1} = E{edge_in {1},2}; L{2} = E{edge_in {2},2};%lengths of in-edges
    x_in{1} = E{edge_in {1},4}; u_in{1} = E{edge_in {1},5};      % states
    x_in{2} = E{edge_in {2},4}; u_in{2} = E{edge_in {2},5};      % before
    x_out   = E{edge_out{1},4}; u_out   = E{edge_out{1},5};  % operations
    % truncate/extrapolate solution to edges
    [xn_in{1},un_in{1},I_in{1}] = fit_to_domain_right(x_in{1},u_in{1},L{1});
    [xn_in{2},un_in{2},I_in{2}] = fit_to_domain_right(x_in{2},u_in{2},L{2});
    [xn_out  ,un_out  ,I_out  ] = fit_to_domain_left(x_out,u_out);
    % solve Riemann problem
    [ub_in{1},ub_in{2},ub_out] = Riemann_problem_21(...
        rhoM1{1},vM1{1},rhoM1{2},vM1{2},rhoM2,vM2,...
        un_in{1}(end),un_in{2}(end),un_out(1));
    for j = 1:2                                               % add jumps
        [xn_in{j},un_in{j}] = add_jump(...
            -xn_in{j}(end:-1:1),un_in{j}(end:-1:1),ub_in{j});
        xn_in{j} = -xn_in{j}(end:-1:1); un_in{j} = un_in{j}(end:-1:1);
    end
    [xn_out,un_out] = add_jump(xn_out,un_out,ub_out);
    % manage areas
    C = u_out(1)*dflux(rhoM2,vM2,u_out(1))-...
        u_in{1}(end)*dflux(rhoM1{1},vM1{1},u_in{1}(end))-...
        u_in{2}(end)*dflux(rhoM1{2},vM1{2},u_in{2}(end));
    I_add = I_out-E{edge_out{1},6}(1)+C*t+...
        I_in{1}-E{edge_in {1},6}(2)+I_in{2}-E{edge_in {2},6}(2);
    if dflux(rhoM2,vM2,ub_out)>=0                  % out-edge is affected
        [xn_out,un_out] = match_area_by_function(xn_out,un_out,I_add);
    elseif dflux(rhoM1{1},vM1{1},ub_in{1})<=0     % in-edge 1 is affected
        [xn_in{1},un_in{1}] = match_area_by_function(...
            -xn_in{1}(end:-1:1),un_in{1}(end:-1:1),I_add);
        xn_in{1} = -xn_in{1}(end:-1:1); un_in{1} = un_in{1}(end:-1:1);
    elseif dflux(rhoM1{2},vM1{2},ub_in{2})<=0     % in-edge 2 is affected
        [xn_in{2},un_in{2}] = match_area_by_function(...
            -xn_in{2}(end:-1:1),un_in{2}(end:-1:1),I_add);
        xn_in{2} = -xn_in{2}(end:-1:1); un_in{2} = un_in{2}(end:-1:1);
    end
    % assign results to edge data structure
    E{edge_in {1},4} = xn_in{1}; E{edge_in {1},5} = un_in{1};
    E{edge_in {2},4} = xn_in{2}; E{edge_in {2},5} = un_in{2};
    E{edge_out{1},4} = xn_out;   E{edge_out{1},5} = un_out;
elseif n_in==0                                  % ingoing boundary vertex
    x_out = E{edge_out{1},4}; u_out = E{edge_out{1},5};
    [xn_out,un_out,~] = fit_to_domain_left(x_out,u_out);
    E{edge_out{1},4} = xn_out; E{edge_out{1},5} = un_out;
elseif n_out==0                                % outgoing boundary vertex
    L = E{edge_in {1},2};
    x_in  = E{edge_in {1},4}; u_in  = E{edge_in {1},5};
    [xn_in ,un_in ,~] = fit_to_domain_right(x_in,u_in,L);
    E{edge_in {1},4} = xn_in;  E{edge_in {1},5} = un_in;
end
newE = E;


function [xnew,unew] = add_jump(x,u,u0)
% Add a jump to the value u0 before the first particle, but do so only
% if u0 is different from u(1).
if abs(u0-u(1))>1e-10
    xnew = [x(1) x(1)+1e-10 x(2:end)];
    unew = [u0,u];
else
    xnew = x; unew = u;
end


function [xnew,unew] = match_area_by_function(x,u,area_add)
% Modify function represented by particles to add area_add to the area.
% successive attempts are undertaken with boxes [x(1),x(j)]x[u(1),u(j)].
if abs(area_add)<1e-12     % if no significant area change has to be done
    xnew = x; unew = u; return               % just leave input unchanged
end
area = 0;
for j = 1:6
    dx = x(j+1)-x(j);                              % size of one interval
    dxcum = x(j+1)-x(1);                    % cumulated size of intervals
    area = area+dx*(u(j)+u(j+1))/2;          % actual area under function
    area_min_2 = min(u(1:j+1))*dxcum;           % minimum area in big box
    area_max_2 = max(u(1:j+1))*dxcum;           % maximum area in big box
    area_min_1 = min(u([1,j+1]))*dxcum;       % minimum area in small box
    area_max_1 = max(u([1,j+1]))*dxcum;       % maximum area in small box
    area_des = area+area_add;                              % desired area
    if area_min_1<=area_des&&area_des<=area_max_1       % small box works
        x_mid = (area_des-x(j+1)*u(j+1)+x(1)*u(1))/(u(1)-u(j+1));
        xnew = [x(1),x_mid,x_mid,x(j+1:end)];       % represent area by a
        unew = [u(1),u(1),u(j+1),u(j+1:end)];      % single discontinuity
        break
    end
    if (area_min_2<=area_des&&area_des<=area_max_2)...    % big box works
            ||j==6            % or accept over/under-shoot as last resort
        u_mid = area_des/dxcum;  % constant value to achieve desired area
        xnew = [x(1),x(1),x(j+1),x(j+1:end)];       % represent area by a
        unew = [u(1),u_mid,u_mid,u(j+1:end)];            % constant piece
        break
    end
end


function [xnew,unew,area] = fit_to_domain_left(x,u)
% Truncates ot extrapolates particle solution at left side of edge x=0.
% Yields also the exact excess/virtual area beyond x=0
area = 0;
while x(2)<0           % remove intervals beyond a single one left of x=0
    area = area+(x(2)-x(1))*(u(1)+u(2))/2;
    x = x(2:end); u = u(2:end);
end
if x(1)<0                                      % one point is left of x=0
    ub = u(1)+(u(2)-u(1))*(-x(1))/(x(2)-x(1));
    area = area+(0-x(1))*(u(1)+ub)/2;       % contributes a positive area
    x = [0,x(2:end)]; u = [ub,u(2:end)];
else                                  % leftmost point is not left of x=0
    ub = u(1);                   % define value by constant extrapolation
    area = (0-x(1))*u(1);                   % contributes a negative area
    x = [0,x]; u = [ub,u];
end
if abs(x(2)-x(1))+abs(u(2)-u(1))<1e-8             %first two points close
    x = x([1 3:end]); u = u([1 3:end]);               % remove second one
end
xnew = x; unew = u;


function [xnew,unew,area] = fit_to_domain_right(x,u,L)
% Truncates ot extrapolates particle solution at right side of edge x=L.
% Yields also the exact excess/virtual area beyond x=L.
area = 0;
while x(end-1)>L      % remove intervals beyond a single one right of x=L
    area = area+(x(end)-x(end-1))*(u(end-1)+u(end))/2;
    x = x(1:end-1); u = u(1:end-1);
end
if x(end)>L                                   % one point is right of x=L
    ub = u(end-1)+(u(end)-u(end-1))*(L-x(end-1))/(x(end)-x(end-1));
    area = area+(x(end)-L)*(ub+u(end))/2;   % contributes a positive area
    x = [x(1:end-1),L]; u = [u(1:end-1),ub];
else                                % rightmost point is not right of x=L
    ub = u(end);                 % define value by constant extrapolation
    area = (x(end)-L)*u(end);               % contributes a negative area
    x = [x,L]; u = [u,ub];
end
if abs(x(end)-x(end-1))+abs(u(end)-u(end-1))<1e-8 % last two points close
    x = x([1:end-2 end]); u = u([1:end-2 end]);   % remove second to last
end
xnew = x; unew = u;
% For future: This could also be calling left sides routine after y = L-x

%------------------------------------------------------------------------

function newE = generalized_Riemann_problem(V,E,j_node)
% Call generalized Riemann problems at j-th network node and add jumps
% at extremal points on adjacent edges, if required. This routine
% assumes that the first particle is at x=0, and the last is at x=L.
edge_in  = V{j_node,1}; n_in  = length(edge_in);
edge_out = V{j_node,2}; n_out = length(edge_out);
if n_in==0||n_out==0 % boundary vertex
    newE = E; return
end
if n_in==1&&n_out==1                                         % bottleneck
    rhoM1 = E{edge_in {1},3}{1}; vM1 = E{edge_in {1},3}{2};
    rhoM2 = E{edge_out{1},3}{1}; vM2 = E{edge_out{1},3}{2};
    u_in {1} = E{edge_in {1},5}(end);        % rightmost value on in-edge
    u_out{1} = E{edge_out{1},5}(1);          % leftmost value on out-edge
    [un_in{1},un_out{1}] = Riemann_problem_11(...              % solve RP
        rhoM1,vM1,rhoM2,vM2,u_in{1},u_out{1});
    % in case of shock right on network node, place it on adjacent edge
    if dflux(rhoM1,vM1,un_in{1})>0&&dflux(rhoM2,vM2,un_out{1})<0
        if abs(un_in{1}-u_in{1})<=1e-10        % in-state has not changed
            f = flux(rhoM1,vM1,un_in{1});
            un_in{1} = new_state(rhoM1,vM1,-1,f,inf);
        elseif abs(un_out{1}-u_out{1})<=1e-10 % out-state has not changed
            f = flux(rhoM2,vM2,un_out{1});
            un_out{1} = new_state(rhoM2,vM2,1,f,inf);
        end
    end
elseif n_in==1&&n_out==2                                    % bifurcation
    rhoM1 = E{edge_in {1},3}{1}; vM1 = E{edge_in {1},3}{2};
    rhoM2 = E{edge_out{1},3}{1}; vM2 = E{edge_out{1},3}{2};
    rhoM3 = E{edge_out{2},3}{1}; vM3 = E{edge_out{2},3}{2};
    a12 = V{j_node,4};                                 % direction matrix
    u_in {1} = E{edge_in {1},5}(end);        % rightmost value on in-edge
    u_out{1} = E{edge_out{1},5}(1);        % leftmost value on out-edge 2
    u_out{2} = E{edge_out{2},5}(1);        % leftmost value on out-edge 3
    [un_in{1},un_out{1},un_out{2}] = Riemann_problem_12(...    % solve RP
        rhoM1,vM1,rhoM2,vM2,rhoM3,vM3,a12,u_in{1},u_out{1},u_out{2});
elseif n_in==2&&n_out==1                                     % confluence
    rhoM1 = E{edge_in {1},3}{1}; vM1 = E{edge_in {1},3}{2};
    rhoM2 = E{edge_in {2},3}{1}; vM2 = E{edge_in {2},3}{2};
    rhoM3 = E{edge_out{1},3}{1}; vM3 = E{edge_out{1},3}{2};
    u_in {1} = E{edge_in {1},5}(end);      % rightmost value on in-edge 1
    u_in {2} = E{edge_in {2},5}(end);      % rightmost value on in-edge 2
    u_out{1} = E{edge_out{1},5}(1);          % leftmost value on out-edge
    [un_in{1},un_in{2},un_out{1}] = Riemann_problem_21(...    % solve RP
        rhoM1,vM1,rhoM2,vM2,rhoM3,vM3,u_in{1},u_in{2},u_out{1});
end
for j = 1:n_in                                         % for all in-edges
    if abs(un_in{j}-u_in{j})>1e-10                 % if state has changed
        E{edge_in {j},4}(end+1) = E{edge_in {j},4}(end);       % add jump
        E{edge_in {j},4}(end-1) = E{edge_in {j},4}(end-1)-1e-10; % at x=L
        E{edge_in {j},5}(end+1) = un_in{j};
    end
end
for j = 1:n_out                                       % for all out-edges
    if abs(un_out{j}-u_out{j})>1e-10               % if state has changed
        E{edge_out{j},4} = E{edge_out{j},4}([1 1:end]);        % add jump
        E{edge_out{j},4}(2) = E{edge_out{j},4}(2)+1e-10;         % at x=0
        E{edge_out{j},5} = E{edge_out{j},5}([1 1:end]);
        E{edge_out{j},5}(1) = un_out{j};
    end
end
newE = E;


function [u1n,u2n] = Riemann_problem_11(rhoM1,vM1,rhoM2,vM2,u1,u2)
% Generalized Riemann problem for one in-edge and one out-edge.
f1d = demand(rhoM1,vM1,u1);                 % demand function for in-edge
f2s = supply(rhoM2,vM2,u2);                % supply function for out-edge
f = min([f1d,f2s]);                   % new flux through the intersection
f1n = f; f2n = f;                                   % new fluxes on edges
u1n = new_state(rhoM1,vM1,-1,f1n,u1);              % new state on in-edge
u2n = new_state(rhoM2,vM2, 1,f2n,u2);             % new state on out-edge


function [u1n,u2n,u3n] = Riemann_problem_12(rhoM1,vM1,...
    rhoM2,vM2,rhoM3,vM3,a12,u1,u2,u3)
% Generalized Riemann problem for one in-edge and two out-edges.
a12 = min(max(a12,1e-14),1-1e-14); a13 = 1-a12;    % cap direction matrix
f1d = demand(rhoM1,vM1,u1);                 % demand function for in-edge
f2s = supply(rhoM2,vM2,u2);              % supply function for out-edge 2
f3s = supply(rhoM3,vM3,u3);              % supply function for out-edge 3
f = min([f1d,f2s/a12,f3s/a13]);       % new flux through the intersection
f1n = f; f2n = a12*f; f3n = a13*f;                  % new fluxes on edges
u1n = new_state(rhoM1,vM1,-1,f1n,u1);              % new state on in-edge
u2n = new_state(rhoM2,vM2, 1,f2n,u2);           % new state on out-edge 2
u3n = new_state(rhoM3,vM3, 1,f3n,u3);           % new state on out-edge 3


function [u1n,u2n,u3n] = Riemann_problem_21(rhoM1,vM1,...
    rhoM2,vM2,rhoM3,vM3,u1,u2,u3)
% Generalized Riemann problem for two in-edges and one out-edge.
f1d = demand(rhoM1,vM1,u1);               % demand function for in-edge 1
f2d = demand(rhoM2,vM2,u2);               % demand function for in-edge 2
f3s = supply(rhoM3,vM3,u3);                % supply function for out-edge
if f1d+f2d<=f3s                 % case (a): sufficient capacity available
    u1n = u1; u2n = u2;                          % new states on in-edges
    f = f1d+f2d;                      % new flux through the intersection
    f3n = f;                                       % new flux on out-edge
else                      % case (b): demand exceeds supply, a jam occurs
    q1 = .5; q2 = 1-q1;                                   % merging ratio
    f = min([f1d/q1,f2d/q2,f3s]);     % new flux through the intersection
    f1n = q1*f; f2n = q2*f;                         % new fluxes on edges
        % Nash-shift (optional)
        f1e = f3s-f2n; f2e = f3s-f1n;                     % excess fluxes
        f1n = min(f1d,f1e); f2n = min(f2d,f2e);  % new fluxes on in-edges
    f3n = f1n+f2n;                                 % new flux on out-edge
    u1n = new_state(rhoM1,vM1,-1,f1n,u1);        % new state on in-edge 1
    u2n = new_state(rhoM2,vM2,-1,f2n,u2);        % new state on in-edge 2
end
u3n = new_state(rhoM3,vM3, 1,f3n,u3);             % new state on out-edge


function un = new_state(rhoM,vM,dir,fn,u)
% New state for in-road (dir=-1) or out-roads (dir=1).
% This amounts to inverting the flux function.
f = flux(rhoM,vM,u);
if abs(fn-f)<1e-12                              % if flux has not changed
    un = u;                          % the new state equals the old state
else                                                          % otherwise
    un = rhoM/2*(1-dir*sqrt(1-4/(vM*rhoM)*fn));  % find new state to flux
end                                   % that moves away from intersection


function f = demand(rhoM,vM,u)
% Demand function for in-road.
f = flux(rhoM,vM,min(u,rhoM/2));


function f = supply(rhoM,vM,u)
% Supply function for out-road.
f = flux(rhoM,vM,max(u,rhoM/2));


function f = flux(rhoM,vM,u)
% Flux function.
f = vM*(1-(1/rhoM)*u)*u;


function f = dflux(rhoM,vM,u)
% Characteristic velocity.
f = vM*(1-(2/rhoM)*u);



%========================================================================
% Plotting Functions
%========================================================================
function plot_network(par,V,E)
if ~par.flag_plot, return, end
t = linspace(0,2*pi,64); t = t(2:end); vrad = 0.05*2;
pcx = cos(t)*vrad; pcy = sin(t)*vrad;
nv = size(V,1); ne = size(E,1);
clf, hold on
rfac = 1.56;
box = zeros(nv+ne,4);
if par.flag_color, col = [1 0 0]; else col = [.3 .4 .4]; end
for j = 1:nv % loop over vertices
    pos = V{j,3};
    patch(pos(1)+pcx*rfac,pos(2)+pcy*rfac,col,'edgecolor','none')
    box(j,:) = [pos(1)+rfac*vrad*[-1 1],pos(2)+rfac*vrad*[-1 1]];
end
for j = 1:ne % loop over edges
    edgeverts = E{j,1};
    ps = V{edgeverts(1),3}; pe = V{edgeverts(2),3};
    L = E{j,2}; x = E{j,4}; u = E{j,5};
    rhoM = E{j,3}{1};
    box(nv+j,:) = plot_cedge(par,ps,pe,L,vrad,x,u/rhoM,rhoM*2);
    vM = E{j,3}{2};
    mid = (ps+pe)/2; % coordinates of midpoint
    text(mid(1),mid(2),sprintf('u_{m}=%g',vM),'fontsize',10,...
        'horizontalalign','center','verticalalign','middle')
end
for j = 1:nv % loop over vertices
    pos = V{j,3};
    patch(pos(1)+pcx,pos(2)+pcy,[0 0 0],'edgecolor','none')
end
hold off, axis equal; caxis([0 1])
ax = [min(box(:,1)),max(box(:,2)),min(box(:,3)),max(box(:,4))];
ax = ax+[-1 1 -1 1]*.01; axis(ax)
title(sprintf('%s at t =%5.2f',par.name,par.t))
if par.flag_save_figure&&abs(par.t-par.tfinal)<1e-8
    res = 200;
    set(gca,'visible','off')
    set(gca,'position',[0 0 1 1])
    set(gcf,'paperpositionmode','auto')
    set(gcf,'position',[10 400 res*(ax(2)-ax(1)) res*(ax(4)-ax(3))])
    framenr = round(par.t/par.dt);
    filename = sprintf('fig_network_%s_%04d',par.name_figure,framenr);
    print(filename,'-dpng','-r320','-opengl');
else
    set(gca,'position',[.04 .04 .92 .92])
    drawnow
%    pause
end


function box = plot_cedge(par,ps,pe,L,r,x,u,plot_height)
% Create ploygonal road from p1 to p2
% with thickness r
% and scalar function colored, given by
% x (positions 0->L) and u (values)
if nargin<1, plot_height = 3; end       % height of bounding box for plot
if par.flag_plot_edge_truncate
    [x,u] = reduce_to_edge(x,u,L);%truncate/extend edge solution to [0,L]
end
v = pe-ps; t = v/norm(v); n = [t(2) -t(1)];          % tangent and normal
ps = ps+t*r*1.2; pe = pe-t*r*1.2;              % move start and end point
nx = length(x);
C = ps'*ones(1,nx)+(pe-ps)'*x/L;                     % center line points
e = ones(1,nx); eo = (-1).^(1:nx);
S1 = C+(n'*eo)*r; S2 = C-(n'*eo)*r;              % oscillating side lines
% Draw colormaps
verts = [reshape([S1(1,:);S2(1,:)],1,[])',...
         reshape([S1(2,:);S2(2,:)],1,[])'];
cdata = reshape([1;1]*u,1,[])';
faces = [1:2:nx*2-3;2:2:nx*2-2;3:2:nx*2-1;4:2:nx*2]';
patch('Faces',faces,'Vertices',verts,...
    'FaceColor','interp','FaceVertexCData',cdata,...
    'EdgeColor','none','CDataMapping','scaled')
if par.flag_color, colormap(jet(256)), else colormap(gray(256)), end
% Draw arrows
la = [.05;.20]; ar = (1-la([1 1 2]))*ps+la([1 1 2])*pe+[n;-n;n*0]*r*.3;
patch(ar(:,1),ar(:,2),[1 1 1]*.2)
la = 1-la([2 1]); ar = (1-la([1 1 2]))*ps+la([1 1 2])*pe+[n;-n;n*0]*r*.3;
patch(ar(:,1),ar(:,2),[1 1 1]*.2)
% Draw function plot
switch par.plot_style
case 1
    SX = C+(n'*e)*r;                               % left sides side line
    plot(SX(1,:),SX(2,:),'k.','markersize',10)           % plot particles
    SX = C-(n'*e)*r;                              % right sides side line
    plot(SX(1,:),SX(2,:),'k.','markersize',10)           % plot particles
    box = [min(SX,[],2)' max(SX,[],2)']; box = box([1 3 2 4]);
case 2
    B = [C(:,[1 end])-(n'*[1 1])*r,...                     % bounding box
        C(:,[1 end])-(n'*[1 1])*(1+plot_height)*r];            % for plot
    B = B(:,[1 2 4 3]);
    SU = C-(n'*(e+u*plot_height))*r;              % added function values
    patch(B(1,:),B(2,:),[.9 .9 1])                    % plot bounding box
    plot(SU(1,:),SU(2,:),'k.-','linewidth',1,'markersize',10)  % and fct.
    box = [min(B,[],2)' max(B,[],2)']; box = box([1 3 2 4]);
end


function [xnew,unew] = reduce_to_edge(x,u,L)
% For the puropse of plotting, bring data (x,u) to the
% interval [0,L].
x = [min(-1,x(1)-1),x,max(L+1,x(end)+1)];           % extend beyond [0,L]
u = u([1 1:end end]);                            % constant extrapolation
ind1 = find(x<0); ind2 = find(x>L);         % make exactly first and last
ind = ind1(end):ind2(1);                    % points lie outside of [0,L]
x = x(ind); u = u(ind);
u(1) = u(1)+(u(2)-u(1))*(-x(1))/(x(2)-x(1));       % interpolate to [0,L]
u(end) = u(end-1)+(u(end)-u(end-1))*(L-x(end-1))/(x(end)-x(end-1));
x([1,end]) = [0,L];
ind = [true abs(diff(x))+abs(diff(u))>1e-10];      % remove double points
xnew = x(ind); unew = u(ind);
