function [nodes,abserr, relerr, stretchchange,numiters,ok] = BSequilibrate( ...
            nodes, edgeends, restlengths, springconst, t, iters, maxerr, finald, ...
            fig, progressCallback )
%[nodes,abserr, relerr, stretchchange,numiters,ok] = BSequilibrate( ...
%            nodes, edgeends, restlengths, springconst, t, iters, maxerr, finald, ...
%            fig, progressCallback )
%   Use the ball and spring method to equilibrate a mesh.

    if ~exist( 'progressCallback', 'var' )
        progressCallback = [];
    end
    dt = t/iters;
    numiters = iters;
    minlen = min(abs(restlengths));
    oldstretches = zeros( size(edgeends,1), 1 );
    ok = true;
    for i=1:iters
        [d,abserr,relerr,stretches] = springdisplacementR0( nodes, edgeends, restlengths, springconst, dt );
        scale = 1; % max( 1, 0.01 * abserr/max(abs(d(:))) );
        stretchchange = max(abs((stretches-oldstretches)./restlengths));
        nodes = nodes + d*scale;
        maxd = max(abs(d(:)));
        relmaxd = maxd/minlen;
        status = [];
        if (~isempty(fig)) && (mod(i,10)==0)
            handles = guidata( fig );
          % plotsimplemesh( handles.plotAxes, nodes, edgeends );
            plotsimplemesh( handles.plotAxes, nodes, edgeends, stretches );
            axis equal
            drawnow;
            if ~isempty(progressCallback)
                status = progressCallback( fig, ...
                        struct( 'iter', i, ...
                                'iters', iters, ...
                                'relstretchchange', stretchchange, ...
                                'target', finald ) );
              % fprintf( 1, 'BSequilibrate progress = %d\n', ok );
            end
        end
        if (~isempty(status)) || (stretchchange <= finald)
            numiters = i;
            break;
        end
        oldstretches = stretches;
    end
end
