function m = splitSecondLayerEdges( m )
%m = splitSecondLayerEdges( m )
%   Split edges of the second layer that are longer than a certain amount.

    edgethreshold = m.secondlayer.splitThreshold;
    edgethresholdsq = edgethreshold*edgethreshold;
    numedges = size(m.secondlayer.edges,1);
    edgevecs = m.secondlayer.cell3dcoords( m.secondlayer.edges(:,1), :) ...
                - m.secondlayer.cell3dcoords( m.secondlayer.edges(:,2), :);
    edgelensqs = sum(edgevecs.*edgevecs,2)';
    splitmap = edgelensqs > edgethresholdsq;
    edgestosplit = find(splitmap);
    if isempty(edgestosplit)
        return;
    end
  % edgestosplit
    numvx = size(m.secondlayer.cell3dcoords,1);
  % perturbationSize = 1/(2*sqrt(3));
    perturbationSize = 0.7/(2*sqrt(3));
    for ei=edgestosplit
        v1 = m.secondlayer.edges(ei,1);
        v2 = m.secondlayer.edges(ei,2);
        p1 = m.secondlayer.cell3dcoords(v1,:);
        p2 = m.secondlayer.cell3dcoords(v2,:);
        p = (p1+p2)*0.5;
        numvx = numvx+1;
        numedges = numedges+1;
        v = numvx;
        fe1 = m.secondlayer.vxFEMcell(v1);
        fe2 = m.secondlayer.vxFEMcell(v2);
        edgevec = p2-p1;
        cn1 = m.unitcellnormals( fe1, : );
        cn2 = m.unitcellnormals( fe2, : );
        cn = (cn1+cn2)*0.5;
        perturbation = cross(edgevec,cn)*perturbationSize;
        if rand(1) < 0.5
            perturbation = -perturbation;
        end
        p = p + perturbation;
        m.secondlayer.cell3dcoords(v,:)= p;
        if fe1==fe2
            m.secondlayer.vxBaryCoords(v,:) = cellBaryCoords( m, fe1, p );
            m.secondlayer.vxFEMcell(v) = fe1;
        else
            [ m.secondlayer.vxFEMcell(v), ...
              m.secondlayer.vxBaryCoords(v,:), ...
              err ] = ...
                findFEMCell( m, p, [fe1 fe2] );
          % err
        end
        m.secondlayer.edges(ei,2) = v;
        m.secondlayer.edges(numedges,:) = ...
            [v,v2,m.secondlayer.edges(ei,[3 4])];
        if isfield( m.secondlayer, 'generation' )
            m.secondlayer.generation(numedges) = m.secondlayer.generation(ei);
        end
        c1 = m.secondlayer.edges(ei,3);
        c2 = m.secondlayer.edges(ei,4);
        if c1 ~= 0
            [m.secondlayer.cells(c1).vxs,m.secondlayer.cells(c1).edges] = ...
                insertvertex( ...
                    m.secondlayer.cells(c1).vxs, ...
                    m.secondlayer.cells(c1).edges, ...
                    numedges, v1, v2, v );
        end
        if c2 ~= 0
            [m.secondlayer.cells(c2).vxs,m.secondlayer.cells(c2).edges] = ...
                insertvertex( ...
                    m.secondlayer.cells(c2).vxs, ...
                    m.secondlayer.cells(c2).edges, ...
                    numedges, v1, v2, v );
        end
    end
    checkclonesvalid( m.secondlayer );
end

function [vxs,edges] = insertvertex( vxs, edges, ei, v1, v2, v )
    v1i = find( vxs==v1 );
    v2i = find( vxs==v2 );
    % Check v2i is next after v1i.
    if v2i ~= (mod(v1i,length(vxs))+1)
        complain('insertvertex: invalid v2i: %d %d', v1i, v2i );
        vxs
        edges
        v1
        v2
        v
    end
    vxs = [ vxs(1:v1i), v, vxs((v1i+1):end) ];
    edges = [ edges(1:v1i), ei, edges((v1i+1):end) ];
end
