function m = computeResiduals( m, retainFraction )
%m = computeResiduals( m, retainFraction )
%   Given m.displacements, this computes the following components of
%   m.celldata:
%       displacementStrain
%       residualStrain
%   If required (i.e. if celldata.actualGrowthTensor exists), it also
%   computes:
%       actualGrowthTensor

    if nargin < 2
        retainFraction = 0;
    end
    numCells = size(m.tricellvxs,1);
    vorticities = [];
    if isempty(m.displacements) || (m.globalProps.timestep==0)
        for ci=1:numCells
            m.celldata(ci).residualStrain = zeros( size(m.celldata(ci).residualStrain) );
        end
    else
        vorticities = zeros(3,3,6,numCells);
        for ci=1:numCells
            trivxs = m.tricellvxs(ci,:);
            prismvxs = [ trivxs*2-1, trivxs*2 ];
            if retainFraction==0
                 m.celldata(ci) = computeDisplacementStrains( ...
                    m.celldata(ci), m.displacements(prismvxs,:) );
            else
                % We only need to compute the vorticities when there is
                % non-zero retained strain, because we need to rotate the
                % residual strain according to the growh in the last step.
                [m.celldata(ci),vort] = computeDisplacementStrains( ...
                    m.celldata(ci), m.displacements(prismvxs,:) );
                vorticities(:,:,:,ci) = vort;
            end
        end
        m = computeResidualStrains( m, retainFraction, vorticities );
    end
end
