function m = splitclonecell( m, ci, v, splitpoint )
%m = splitclonecell( m, ci, v )
%   Split cell ci perpendicular to direction v, through splitpoint.
%   splitpoint defaults to the centroid of the vertexes.

%   One new cell, two new vertexes, and three new edges are created.

%   dumpsecondlayer( m.secondlayer );

    if all(v==0)
        return;
    end

    % Should pull-in be implemented for edges that do not have a cell on
    % both sides?
    pullRim = false;
    
    % Force edges to be split at least this distance from their endpoints.
    min_d = 0.05 * sqrt(m.secondlayer.averagetargetarea/1.4);

VERBOSE = 0;
if VERBOSE
fprintf( 1, 'splitclonecell: splitting cell %d.\n', ci );
end
    % Find out how many of everything there is, and get the cell vertex
    % coordinates.
    numcells = length( m.secondlayer.cells );
    numvertexes = length( m.secondlayer.vxFEMcell );
    numedges = size( m.secondlayer.edges, 1 );
    cell3dcoords = m.secondlayer.cell3dcoords( m.secondlayer.cells(ci).vxs, : );
    numcellvxs = size( cell3dcoords, 1 );

    [split1cv1,split1cv2,newvx3d1,split2cv1,split2cv2,newvx3d2] = splitPoly( v, splitpoint, cell3dcoords );
    if split1cv1 == 0
        fprintf( 1, ...
            '%s: Error: cell %d is not split by a plane through the splitpoint: c (%.3f %.3f %.3f) v (%.3f %.3f %.3f).\n', ...
            mfilename(), ci, splitpoint, v );
        return;
    end
    splitvxs = m.secondlayer.cells(ci).vxs( [ split1cv1, split1cv2, split2cv1, split2cv2 ] );
    splitedgeindexes = m.secondlayer.cells(ci).edges( [ split1cv1, split2cv1 ] );
    splitedgedata = m.secondlayer.edges( splitedgeindexes, : );
% if VERBOSE
%     fprintf( 1, 's1cv1 %d s2cv1 %d s1cv2 %d s2cv2 %d\n', ...
%         split1cv1, split2cv1, split1cv2, split2cv2 );
%     fprintf( 1, 'splitvxs [ %d %d %d %d ]\n', splitvxs );
%     fprintf( 1, 'splitedgeindexes [ %d %d ]\n', splitedgeindexes );
%     fprintf( 1, 'cellvxs [' );
%     fprintf( 1, ' %d', m.secondlayer.cells(ci).vxs );
%     fprintf( 1, ' ]\n' );
%     fprintf( 1, 'celledges [' );
%     fprintf( 1, ' %d', m.secondlayer.cells(ci).edges );
%     fprintf( 1, ' ]\n' );
% end
    % Create indexes for the new cell, edges, and vertexes.
    newci = numcells+1;
    newvi1 = numvertexes+1;
    newvi2 = numvertexes+2;
    newei1 = numedges+1;
    newei2 = numedges+2;
    newei3 = numedges+3;
    lastGeneration = m.secondlayer.generation( length(m.secondlayer.generation) );
    m.secondlayer.generation([newei1,newei2,newei3]) = ...
         [ m.secondlayer.generation(splitedgeindexes(1));
           m.secondlayer.generation(splitedgeindexes(2));
           lastGeneration + 1 ];
% if VERBOSE
% fprintf( 1, 'newci %d newvi1 %d newvi2 %d newei1 %d newei2 %d newei3 %d\n', ...
%     newci, newvi1, newvi2, newei1, newei2, newei3 );
% end
% if VERBOSE
% fprintf( 1, 'splitedgeindexes %d %d\n', ...
%     splitedgeindexes );
% fprintf( 1, 'splitedgedata [ %d %d %d %d ] [ %d %d %d %d ]\n', ...
%     splitedgedata' );
% end
    % Find the cells that lie on the other side of each of the edges to be
    % split.
    otherci1 = splitedgedata( 1, 3 );
    if otherci1==ci
        otherci1 = splitedgedata( 1, 4 );
    end
    
    otherci2 = splitedgedata( 2, 3 );
    if otherci2==ci
        otherci2 = splitedgedata( 2, 4 );
    end
    
    % Find the intersections of the split edges with the cutting plane.
    % These are the new vertexes.
    edgeendCoords = m.secondlayer.cell3dcoords( splitvxs, : );
%     newvx3d1 = v1; % lineplaneIntersection( edgeendCoords(1,:), edgeendCoords(2,:), v, splitpoint )
%     newvx3d2 = v2; % lineplaneIntersection( edgeendCoords(3,:), edgeendCoords(4,:), v, splitpoint )

    newvx3d1 = avoidEnds( newvx3d1, edgeendCoords(1,:), edgeendCoords(2,:), min_d );
    newvx3d2 = avoidEnds( newvx3d2, edgeendCoords(3,:), edgeendCoords(4,:), min_d );

    b = m.globalProps.bioApullin;
    a = 1 - b;
    if (~pullRim) && (otherci1 == 0)
        adjnewvx3d1 = newvx3d1;
    else
        adjnewvx3d1 = a*newvx3d1 + b*newvx3d2;
    end
    if (~pullRim) && (otherci2 == 0)
        adjnewvx3d2 = newvx3d2;
    else
        adjnewvx3d2 = a*newvx3d2 + b*newvx3d1;
    end
    newvx3d1 = adjnewvx3d1;
    newvx3d2 = adjnewvx3d2;
    
    % Find the parent FEM cell and barycentric coordinates for the first
    % new vertex.
    % Find which FEM cells the edge endpoints are in.
    femCells = m.secondlayer.vxFEMcell( splitvxs );
if VERBOSE
    fprintf( 1, 'femCells %d %d\n', femCells );
end
    if 0
        [ testvi1, testbc1 ] = splitcloneedge( m, splitedgeindexes(1) );
        [ testvi2, testbc2 ] = splitcloneedge( m, splitedgeindexes(2) );
    end
    if 0 && (testvi1 ~= 0) && (testvi2 ~= 0)
        fprintf( 1, 'splitclonecell: new method used for cell %d.\n', ci );
        newcell1 = testvi1;
        bc1 = testbc1;
    else
        if femCells(1)==femCells(2)
            newcell1 = femCells(1);
            bc1 = cellBaryCoords( m, newcell1, newvx3d1 );
        else
            [ newcell1, bc1, err ] = findFEMCell( m, newvx3d1, femCells(1:2) );
        end
if VERBOSE
        fprintf( 1, 'newcell1 %d\n', newcell1 );
end
    end
    
    if isempty(bc1)
        % Error: cannot find new point.
        fprintf( 1, 'splitclonecell: cannot find edge midpoint 1.  err %.3f vx3d1 (%.3f %.3f %.3f), cells %d %d\n', ...
            err, newvx3d1, femCells(1:2));
        return;
    end
    
    % Find the parent FEM cell and barycentric coordinates for the second
    % new vertex.
    if 0 && (testvi1 ~= 0) && (testvi2 ~= 0)
        newcell2 = testvi2;
        bc2 = testbc2;
    else
        if femCells(3)==femCells(4)
            newcell2 = femCells(3);
            bc2 = cellBaryCoords( m, newcell2, newvx3d2 );
        else
            [ newcell2, bc2, err ] = findFEMCell( m, newvx3d2, femCells(3:4) );
        end
if VERBOSE
        fprintf( 1, 'newcell2 %d\n', newcell2 );
end
    end
    
    if isempty(bc2)
        % Error: cannot find new point.
        fprintf( 1, 'splitclonecell: cannot find edge midpoint 2.  err %.3f vx3d2 (%.3f %.3f %.3f), cells %d %d\n', ...
            err, newvx3d2, femCells(3:4));
        return;
    end

    % Create the new vertexes.
    m.secondlayer.vxFEMcell(newvi1) = newcell1;
    m.secondlayer.vxBaryCoords(newvi1,:) = normaliseBaryCoords( bc1 );
    m.secondlayer.vxFEMcell(newvi2) = newcell2;
    m.secondlayer.vxBaryCoords(newvi2,:) = normaliseBaryCoords( bc2 );
    
    % Create the new edges, and modify the existing ones.
    m.secondlayer.edges(newei1,:) = [ splitvxs(2), newvi1, newci, otherci1 ];
    m.secondlayer.edges(newei2,:) = [ splitvxs(4), newvi2, ci, otherci2 ];
    m.secondlayer.edges(newei3,:) = [ newvi1, newvi2, ci, newci ];
    m.secondlayer.edges(splitedgeindexes(1),1:2) = ...
        replaceEdgeEnd( m.secondlayer.edges(splitedgeindexes(1),1:2), ...
            splitvxs(2), newvi1 );
    m.secondlayer.edges(splitedgeindexes(2),1:2) = ...
        replaceEdgeEnd( m.secondlayer.edges(splitedgeindexes(2),1:2), ...
            splitvxs(4), newvi2 );
    if split2cv2 <= split1cv1
        oldrange = split2cv2:split1cv1;
    else 
        oldrange = [split2cv2:numcellvxs, 1:split1cv1];
    end
    if split1cv2 <= split2cv1
        newrange = split1cv2:split2cv1;
    else 
        newrange = [split1cv2:numcellvxs, 1:split2cv1];
    end
if VERBOSE
    fprintf( 1, 'newei1 %d data [ %d %d %d %d ]\n', newei1, m.secondlayer.edges(newei1,:) );
    fprintf( 1, 'newei2 %d data [ %d %d %d %d ]\n', newei2, m.secondlayer.edges(newei2,:) );
    fprintf( 1, 'newei3 %d data [ %d %d %d %d ]\n', newei3, m.secondlayer.edges(newei3,:) );
end
    % Create the new cell's vertex and edge lists.
    m.secondlayer.cells(newci).vxs = ...
        [ m.secondlayer.cells(ci).vxs( newrange ), newvi2, newvi1 ];
    m.secondlayer.cells(newci).edges = ...
        [ m.secondlayer.cells(ci).edges( newrange ), newei3, newei1 ];
if VERBOSE
fprintf( 1, 'newci %d vxs [', newci );
fprintf( 1, ' %d', m.secondlayer.cells(newci).vxs );
fprintf( 1, ' ]\n' );
fprintf( 1, 'newci %d edges [', newci );
fprintf( 1, ' %d', m.secondlayer.cells(newci).edges );
fprintf( 1, ' ]\n' );
end

    % Make the new cell's edges refer to the new cell.
    for ei=m.secondlayer.cells(ci).edges( newrange )
    %    oldedgesei = m.secondlayer.edges(ei,:)
        m.secondlayer.edges(ei,3:4) = ...
            replaceEdgeEnd( m.secondlayer.edges(ei,3:4), ci, newci );
    end
if VERBOSE
fprintf( 1, 'newci %d edge data [\n', newci );
fprintf( 1, '    v %d e %d: %d %d %d %d\n', ...
    [ m.secondlayer.cells(newci).vxs; ...
      m.secondlayer.cells(newci).edges; ...
      m.secondlayer.edges(m.secondlayer.cells(newci).edges,:)' ] );
fprintf( 1, ' ]\n' );
end
    % Modify the old cell's vertex and edge lists.
    m.secondlayer.cells(ci).vxs = ...
        [ m.secondlayer.cells(ci).vxs( oldrange ), newvi1, newvi2 ];
    m.secondlayer.cells(ci).edges = ...
        [ m.secondlayer.cells(ci).edges( oldrange ), newei3, newei2 ];
if VERBOSE
fprintf( 1, 'ci %d edge data [\n', ci );
fprintf( 1, '    v %d e %d: %d %d %d %d\n', ...
    [ m.secondlayer.cells(ci).vxs; ...
      m.secondlayer.cells(ci).edges; ...
      m.secondlayer.edges(m.secondlayer.cells(ci).edges,:)' ] );
fprintf( 1, ' ]\n' );
end

    % Insert the new vertexes and edges into otherci1 and otherci2.
    % Find splitvxs(1:2) in otherci1.vxs.  Insert newvi1 between.
    % Insert newei1 one place before.
    if otherci1 ~= 0
        numouter1vxs = length(m.secondlayer.cells(otherci1).vxs);
        i = find( m.secondlayer.cells(otherci1).vxs==splitvxs(2), 1 );
        i1 = mod(i,numouter1vxs) + 1;
        if m.secondlayer.cells(otherci1).vxs(i1) == splitvxs(1)
if VERBOSE
            fprintf( 1, 'Same orientation for cell otherci1 %d.\n', otherci1 );
end
            m.secondlayer.cells(otherci1).vxs = ...
                [ m.secondlayer.cells(otherci1).vxs(1:i), ...
                  newvi1, ...
                  m.secondlayer.cells(otherci1).vxs(i+1:numouter1vxs) ];
            m.secondlayer.cells(otherci1).edges = ...
                [ m.secondlayer.cells(otherci1).edges(1:i-1), ...
                  newei1, ...
                  m.secondlayer.cells(otherci1).edges(i:numouter1vxs) ];
        else
if VERBOSE
            fprintf( 1, 'Opposite orientation for cell otherci1 %d.\n', otherci1 );
end
            m.secondlayer.cells(otherci1).vxs = ...
                [ m.secondlayer.cells(otherci1).vxs(1:i-1), ...
                  newvi1, ...
                  m.secondlayer.cells(otherci1).vxs(i:numouter1vxs) ];
            m.secondlayer.cells(otherci1).edges = ...
                [ m.secondlayer.cells(otherci1).edges(1:i-1), ...
                  newei1, ...
                  m.secondlayer.cells(otherci1).edges(i:numouter1vxs) ];
        end
    end

      
    % Repeat for the other cell.
    if otherci2 ~= 0
        numouter2vxs = length(m.secondlayer.cells(otherci2).vxs);
        i = find( m.secondlayer.cells(otherci2).vxs==splitvxs(4), 1 );
        i1 = mod(i,numouter2vxs) + 1;
        if m.secondlayer.cells(otherci2).vxs(i1) == splitvxs(3)
if VERBOSE
            fprintf( 1, 'Same orientation for cell otherci2 %d.\n', otherci2 );
end
            m.secondlayer.cells(otherci2).vxs = ...
                [ m.secondlayer.cells(otherci2).vxs(1:i), ...
                  newvi2, ...
                  m.secondlayer.cells(otherci2).vxs(i+1:numouter2vxs) ];
            m.secondlayer.cells(otherci2).edges = ...
                [ m.secondlayer.cells(otherci2).edges(1:i-1), ...
                  newei2, ...
                  m.secondlayer.cells(otherci2).edges(i:numouter2vxs) ];
        else
if VERBOSE
            fprintf( 1, 'Opposite orientation for cell otherci2 %d.\n', otherci2 );
fprintf( 1, 'otherci2 %d cellvxs [', otherci2 );
fprintf( 1, ' %d', m.secondlayer.cells(otherci2).vxs );
fprintf( 1, ' ]\n' );
fprintf( 1, 'otherci2 %d celledges [', otherci2 );
fprintf( 1, ' %d', m.secondlayer.cells(otherci2).edges );
fprintf( 1, ' ]\n' );
fprintf( 1, 'otherci2 %d edge data [\n', otherci2 );
fprintf( 1, '    %d %d %d %d\n', m.secondlayer.edges(m.secondlayer.cells(otherci2).edges,:)' );
fprintf( 1, ' ]\n' );
end
            m.secondlayer.cells(otherci2).vxs = ...
                [ m.secondlayer.cells(otherci2).vxs(1:i-1), ...
                  newvi2, ...
                  m.secondlayer.cells(otherci2).vxs(i:numouter2vxs) ];
            m.secondlayer.cells(otherci2).edges = ...
                [ m.secondlayer.cells(otherci2).edges(1:i-1), ...
                  newei2, ...
                  m.secondlayer.cells(otherci2).edges(i:numouter2vxs) ];
if VERBOSE
fprintf( 1, 'otherci2 %d cellvxs [', otherci2 );
fprintf( 1, ' %d', m.secondlayer.cells(otherci2).vxs );
fprintf( 1, ' ]\n' );
fprintf( 1, 'otherci2 %d celledges [', otherci2 );
fprintf( 1, ' %d', m.secondlayer.cells(otherci2).edges );
fprintf( 1, ' ]\n' );
fprintf( 1, 'otherci2 %d edge data [\n', otherci2 );
fprintf( 1, '    %d %d %d %d\n', m.secondlayer.edges(m.secondlayer.cells(otherci2).edges,:)' );
fprintf( 1, ' ]\n' );
end
        end
        if 0
            if m.secondlayer.cells(otherci2).vxs(i1) ~= splitvxs(3)
                i1 = i-1;  if i1==0, i1 = numouter2vxs; end
                if m.secondlayer.cells(otherci2).vxs(i1) == splitvxs(3)
                    fprintf( 1, 'Warning: orientation violation for cells %d and %d.\n', ci, otherci2 );
                    i = i1;
                else
                    fprintf( 1, 'Warning: inconsistency for cells %d and %d.\n', ci, otherci2 );
                end
            end
            m.secondlayer.cells(otherci2).vxs = ...
                [ m.secondlayer.cells(otherci2).vxs(1:i), ...
                  newvi2, ...
                  m.secondlayer.cells(otherci2).vxs(i+1:numouter2vxs) ];
            m.secondlayer.cells(otherci2).edges = ...
                [ m.secondlayer.cells(otherci2).edges(1:i-1), ...
                  newei2, ...
                  m.secondlayer.cells(otherci2).edges(i:numouter2vxs) ];
        end
    end

    m = calcCloneVxCoords( m, [newvi1,newvi2] );
    if ~isempty( m.secondlayer.cellcolor )
        numcolorparams = size( m.globalProps.colorparams, 1 );
        if (m.secondlayer.cloneindex( ci ) == 0) ...
                || (numcolorparams < m.secondlayer.cloneindex( ci ))
            newcolor = m.secondlayer.cellcolor( ci, : );
        else
            if numcolorparams==1
                colorindex = 1;
            else
                colorindex = m.secondlayer.cloneindex( ci );
                if colorindex > numcolorparams
                    colorindex = numcolorparams;
                end
            end
            if true
                newcolor = randcolornear( 1, ...
                                m.secondlayer.cellcolor( ci, : ), ...
                                m.globalProps.colorvariation );
                if VERBOSE
                    fprintf( 1, 'old [ %.3f %.3f %.3f ]\nnew [ %.3f %.3f %.3f ]  %f\n\n', ....
                        m.secondlayer.cellcolor( ci, : ), ...
                        newcolor, m.globalProps.colorvariation );
                end
            else
                newcolor = secondlayercolor( 1, ...
                	m.globalProps.colorparams( colorindex, : ) );
            end
        end
        m.secondlayer.cellcolor( newci, : ) = newcolor;
    end

    m.secondlayer.celltargetarea( newci ) = m.secondlayer.celltargetarea( ci );
    m.secondlayer.areamultiple( newci ) = m.secondlayer.areamultiple( ci );
    m.secondlayer.cloneindex( newci ) = m.secondlayer.cloneindex( ci );
    m.secondlayer.side( newci, 1 ) = m.secondlayer.side( ci, 1 );
    m.secondlayer.cellarea( newci ) = polyarea3( ...
            m.secondlayer.cell3dcoords( m.secondlayer.cells(newci).vxs, : ) );

    if isa( m.globalProps.bioApostsplitproc, 'function_handle' )
        m = m.globalProps.bioApostsplitproc( m, ...
                ci, length(m.secondlayer.cells(ci).edges) - 1, ...
                newci, length(m.secondlayer.cells(newci).edges) - 1, ...
                splitedgeindexes(1), splitedgeindexes(2), newei1, newei2, newei3 );
    end
end

function p0 = avoidEnds( p0, p1, p2, min_d )
    d11 = norm( p0 - p1 );
    d12 = norm( p0 - p2 );
    d = norm( p2 - p1 );
    minalpha = min_d/d;
    
    alpha = 0;
    beta = 0;
    if d11==0
        alpha = 1 - minalpha;
        beta = minalpha;
    elseif d12==0
        alpha = minalpha;
        beta = 1 - minalpha;
    elseif d11/d12 < minalpha
        alpha = 1 - minalpha;
        beta = minalpha;
    elseif d12/d11 < minalpha
        alpha = minalpha;
        beta = 1 - minalpha;
    end
    if alpha > 0
        p0 = alpha*p1 + beta*p2;
    end
end
