function m = calcPolGrad( m, cis )
%m = calcPolGrad( m )
%   Calculate the polarising gradients for growth and bending in the
%   specified cells of the mesh, by default all of them.

    if size( m.morphogens, 2 ) < 3
        return;
    end
    
    if m.globalProps.userpolarisation
        return;
    end
    
    if nargin < 2
        cis = 1:size(m.tricellvxs,1);
    end
    
    pol_mgen = polariserIndex( m );
    
    rawgrad = mgenCellGradient( m, pol_mgen, cis );
    rawgradsq = sum( rawgrad.^2, 2 );
    if m.globalProps.relativepolgrad
        grads = perVertextoperFE( m, m.morphogens(:,pol_mgen), cis );
        biggrads = rawgradsq >= (grads .* m.globalProps.mingradient).^2;
    else
        biggrads = rawgradsq >= m.globalProps.mingradient^2;
    end
    bigcis = cis(biggrads);
    m.gradpolgrowth( bigcis, : ) = rawgrad( biggrads, : );
    m.polfreeze( bigcis, : ) = reshape( ...
                                    m.morphogens( m.tricellvxs(bigcis,:), pol_mgen ), ...
                                    length(bigcis), ...
                                    size(m.tricellvxs,2) );
    m.polfreeze = procrustesHeight( m.polfreeze, size(m.tricellvxs,1) );
    m.polfreezebc = procrustesHeight( m.polfreezebc, size(m.tricellvxs,1) );
    bigcisindexes = find(bigcis);
    for i=1:length(bigcisindexes)
        ci = bigcisindexes(i);
        m.polfreezebc(ci,:) = vec2bc( m.gradpolgrowth( ci, : ), ...
                                       m.nodes( m.tricellvxs(ci,:), : ), ...
                                       m.unitcellnormals(ci,:) );
    end
    smallcis = cis(~biggrads);
    if m.globalProps.usefrozengradient
        for ci=smallcis
            vxs = m.tricellvxs(ci,:);
            if m.globalProps.usepolfreezebc
                m.gradpolgrowth( ci, : ) = bc2vec( m.polfreezebc(ci,:), m.nodes( vxs, : ) );
            else
                m.gradpolgrowth( ci, : ) = -trianglegradient( m.nodes( vxs, : ), m.polfreeze(ci,:) );
            end
        end
    else
        m.gradpolgrowth( smallcis, : ) = 0;
    end
    m.polfrozen(cis) = ~biggrads(:);
end
