function sqdist_contour(theta, samples, reality)
if ~exist('reality','var')
    reality = [nan nan];
end
t = reshape(theta, [1 1 length(theta)]);
[r,f] = meshgrid(linspace(0,0.99,samples), linspace(0,2*pi,samples*pi));
x = r.*cos(f);
y = r.*sin(f);
S = sum(log(sqdist(t, x, y)), 3);

mean = [sum(cos(t)) sum(sin(t))] / length(t);
[M I] = max(S(:));
[X Y] = ind2sub(size(x),I);
max_likelihood = [x(X,Y) y(X,Y)];

hold on
axis equal
scatter([mean(1) max_likelihood(1)], [mean(2) max_likelihood(2)]);
text([mean(1) max_likelihood(1) reality(1)], [mean(2) max_likelihood(2) reality(2)], ["MM", "MMV", "\mu"], 'VerticalAlignment', 'top', 'HorizontalAlignment', 'center');
contour(x,y,S);
scatter(cos(t),sin(t),'X')
scatter(reality(1), reality(2), '*', 'black')
hold off