 function [mesh,ok] = addRandomSecondLayerCell( mesh, ...
                        diameter, axisratio, numsides, cloneindex, ...
                        ci, newci, allowoverlap, cellcolor )
%mesh = addRandomSecondLayerCell( mesh, ...
%                       diameter, axisratio, numsides, cloneindex, ...
%                       ci, newci, allowoverlap, cellcolor )
%   Add a second layer cell at a random position within the FEM cell ci. If
%   ci is omitted, a random FEM cell will be chosen.
%   The index of the new cell will be newci; by default, this will be 1
%   more than the number of existing second layer cells.
%   The new cell will have NUMSIDES sides.  DIAMETER is the diameter of the cell,
%   in absolute units.
%   If allowoverlap is false, it will not create the new cell unless it does
%   not overlap any existing cell.

    ok = true;
    if (nargin < 4) || isempty(cloneindex)
        cloneindex = 1;
    end
    if (nargin < 5) || isempty(ci)
        ci = randint( 1, length(mesh.celldata) );
    end
    if (nargin < 6) || isempty(newci)
        newci = numcellsBioA( mesh.secondlayer ) + 1;
    end
    if (nargin < 7) || isempty(allowoverlap)
        allowoverlap = true;
    end
    
    bccentre = randBaryCoords(1);  % [1 1 1]/3;
    femVxs = mesh.nodes( mesh.tricellvxs( ci, : ), : );
    cellcentre = bccentre * femVxs;
    cellradius = diameter/2; % sz * sqrt(mesh.globalDynamicProps.currentArea);
    % Test if cell center is within cellradius of any other vertex of the
    % bio layer.
    if ~allowoverlap
        cellradiussq = cellradius*cellradius;
        for i=1:size(mesh.secondlayer.cell3dcoords,1)
            if (mesh.secondlayer.cell3dcoords(i,:) - cellcentre).^2 < cellradiussq
                % Do not create new cell.
                ok = false;
                return;
            end
        end
    end
    havegradient = any( mesh.gradpolgrowth(ci,:) ~= 0 );
    if havegradient && (axisratio ~= 1)
        if ~isfield( mesh, 'cellFrames' )
            mesh = makeCellFrames( mesh );
        end
      % J = mesh.celldata(ci).cellFrame(:,[3,1,2]);
        J = mesh.cellFrames(:,[3,1,2],ci);
        adjustment = sqrt( axisratio );
        cellpts = circlepoints( cellradius, cellcentre, numsides, ...
            0, J(:,2)'*adjustment, J(:,3)'/adjustment );
    else
        n = mesh.unitcellnormals(ci,:);
        J = makebasis( n );
        cellpts = circlepoints( cellradius, cellcentre, numsides, ...
            rand(1), J(:,2)', J(:,3)' );
    end
    cis = zeros(size(cellpts,1),1);
    bcs = zeros(size(cellpts,1),3);
    bcerr = zeros(size(cellpts,1),1);
    abserr = zeros(size(cellpts,1),1);
    for i = 1:size(cellpts,1)
        [ cis(i), bcs(i,:), bcerr(i), abserr(i) ] = findFEMCell( mesh, cellpts(i,:), ci );
    end
  % cis
  % bcs
  % bcerr
  % abserr

    numvertexes = length( mesh.secondlayer.vxFEMcell );
    numedges = size( mesh.secondlayer.edges, 1 );

    newvi = numvertexes+1 : numvertexes+numsides;
    newei = numedges+1 : numedges+numsides;
    
    if cloneindex==1
        colorindex = 1;
    else
        colorindex = 2;
    end
    colorindex = trimnumber( 1, colorindex, size( mesh.globalProps.colorparams, 1 ) );
    if ~isempty(cellcolor)
        mesh.secondlayer.cellcolor(newci,:) = cellcolor(1,:);
    elseif isempty(mesh.secondlayer.cells) || ~isempty( mesh.secondlayer.cellcolor )
        mesh.secondlayer.cellcolor(newci,:) = ...
            secondlayercolor( 1, mesh.globalProps.colorparams(colorindex,:) );
    end
    mesh.secondlayer.cells(newci).vxs = newvi;
    mesh.secondlayer.cells(newci).edges = newei;
    mesh.secondlayer.cloneindex(newci,:) = cloneindex;
    mesh.secondlayer.side(newci,1) = rand(1) < 0.5;
    
    mesh.secondlayer.vxFEMcell(newvi) = cis;
    mesh.secondlayer.vxBaryCoords(newvi,:) = bcs;
    
    newedges = int32( [ newvi', ...
                        [ newvi(2:end)'; numvertexes+1], ...
                        ones(numsides,1)*newci, ...
                        zeros(numsides,1) ...
                      ] );
    mesh.secondlayer.edges( newei, : ) = newedges;
    mesh = calcCloneVxCoords( mesh, newvi );
    mesh = setSplitThreshold( mesh, 1.05, newci );
end

    
    
