function [p] = plotPhaseUncertainty(W, freqs, args)
% plotPhaseUncertainty -
%
% Syntax: [p] = plotPhaseUncertainty(W, freqs, args)
%
% Inputs:
%    - W     - Multiplicative Uncertainty Weight
%    - freqs - Frequency Vector [Hz]
%    - args  - Optional Arguments:
%      - G
%      - color_i
%      - opacity
%
% Outputs:
%    - p - Plot Handle

    arguments
        W
        freqs double {mustBeNumeric, mustBeNonnegative}
        args.G = tf(1)
        args.unwrap    logical {mustBeNumericOrLogical} = false
        args.color_i (1,1) double {mustBeInteger, mustBeNonnegative} = 0
        args.opacity (1,1) double {mustBeNumeric, mustBePositive} = 0.3
        args.DisplayName char = ''
    end

    % Get defaults colors
    colors = get(groot, 'defaultAxesColorOrder');

    % Compute Phase Uncertainty
    Dphi = 180/pi*asin(abs(squeeze(freqresp(W, freqs, 'Hz'))));
    Dphi(abs(squeeze(freqresp(W, freqs, 'Hz'))) > 1) = 360;

    % Compute Plant Phase
    if args.unwrap
        G_ang = 180/pi*unwrap(angle(squeeze(freqresp(args.G, freqs, 'Hz'))));
    else
        G_ang = 180/pi*angle(squeeze(freqresp(args.G, freqs, 'Hz')));
    end

    p = patch([freqs flip(freqs)], [G_ang+Dphi; flip(G_ang-Dphi)], 'w', ...
              'DisplayName', args.DisplayName);

    if args.color_i == 0
        p.FaceColor = [0; 0; 0];
    else
        p.FaceColor = colors(args.color_i, :);
    end
    p.EdgeColor = 'none';
    p.FaceAlpha = args.opacity;

end