function [I, A] = coverage(a, b, N, flag)

% Computes the coverage (area) of the polar pattern
% produced by a mic with magnitude response: 
% r = |a + b*cos(t)|, a + b = 1
% using a Monte Carlo integration with N points.
% The input flag (1 or 0) is used to toggle plotting.
% A is the estimated area.
% I is the exact area.

n = 0;  % Hits count.
hits = [];
misses = [];

% Compute integral:
for i = 1:N
    r = rand;  % Try different distributions on [0, 1] to weight larger r values.
    theta = rand*2*pi;
    if r <= abs(a + b*cos(theta))
        hits = [hits, [r; theta]];
        n = n + 1;
    else
        misses = [misses, [r; theta]];
    end
end

I = pi*(n/N); % Estimated area.
A = pi*(a^2 + (b^2)/2); % Exact area.

% Plot results if flag == 1:
if flag == 1
    t = 0:0.01:2*pi;
    r = abs(a + b*cos(t));
    h0 = polar(t, r, 'r');
    set(h0, 'LineWidth', 2)
    title(['Estimated Area: ', num2str(I)])
    xlabel(['Exact Area: ', num2str(A)])
    hold on
    h1 = polar(hits(2,:), hits(1,:), 'r.');
    set(h1, 'MarkerSize', 1)
    h2 = polar(misses(2,:), misses(1,:), 'b.');
    set(h2, 'MarkerSize', 1)
end