function varargout = gtGetMaxDisorientation(gvdm, Symm, mode)
    if (nargout == 0)
        silent = false;
    else
        silent = true;
    end
    if (~exist('Symm', 'var') || isempty(Symm))
        Symm = gtCrystGetSymmetryOperators('cubic');
    end
    if (~exist('mode', 'var'))
        mode = 'minimal';
    end
    if (~silent)
        fprintf(' * Computing maximum internal disorientation angle (mode: %s)..', mode);
        t = tic();
    end

    try
        [vec_1, vec_2] = getExtremes(gvdm, mode, silent);

        max_dis_angle = gtDisorientation(vec_1, vec_2, Symm, 'mode', 'active');
    catch mexc
        gtPrintException(mexc, 'Probably no spread in the distribution')

        max_dis_angle = 0;
    end

    if (silent)
        varargout{1} = max_dis_angle;
    else
        fprintf('\b\b: %f deg (in: %f s)\n', max_dis_angle, toc(t));
    end
end

function [vec_1, vec_2] = getExtremes(gvdm, mode, silent)

    num_vecs = size(gvdm, 2);

    switch (mode)
        case 'diameter_zero'
            lengths = sum(gvdm .^ 2, 1);
            [~, i_max] = max(lengths);

            vec_1 = gvdm(:, i_max);
            vec_2 = -gvdm(:, i_max);
        case 'diameter_average'
            center = mean(gvdm, 2);

            lengths = sum((gvdm - center(:, ones(num_vecs, 1))) .^ 2, 1);
            [~, i_max] = max(lengths);

            vec_1 = gvdm(:, i_max);
            vec_2 = center - (gvdm(:, i_max) - center);
        case 'minimal'
            best_match = struct('length', 0.0, 'vec1', 1, 'vec2', 1);
            if (~silent)
                fprintf('\b\b');
                num_chars = fprintf(': ');
            end

            k = convhull(double(gvdm'));
            gvdm = gvdm(:, k);
            num_vecs = size(gvdm, 2);

            for ii = 1:num_vecs
                if (~silent && ~mod(ii, 100))
                    num_chars = fprintf('%d/%d', ii, num_vecs);
                end
                rod_vecs_1 = gvdm(:, ii * ones(num_vecs, 1));

                rod_vecs_1 = rod_vecs_1 - gvdm;
                lengths = sum(rod_vecs_1 .^ 2, 1);
                [max_length, i_max] = max(lengths);

                if (best_match.length < max_length)
                    best_match.length = max_length;
                    best_match.vec1 = ii;
                    best_match.vec2 = i_max;
                end
                if (~silent && ~mod(ii, 100))
                    fprintf(repmat('\b', [1 num_chars]));
                end
            end
            if (~silent)
                fprintf([repmat(' ', [1, num_chars]) repmat('\b', [1, num_chars]) '\b\b..'])
            end

            vec_1 = gvdm(:, best_match.vec1);
            vec_2 = gvdm(:, best_match.vec2);
        otherwise
            error('gtGetMaxDisorientation:wrong_argument', 'No such mode: "%s"', mode)
    end
end