function varargout = checkerplot(f, varargin)
%CHECKERPLOT   Show the complex transform of some grid pattern.
%
%   CHECKERPLOT(f, __), with f a function handle, draws a checkerboard
%      pattern's tranform by f(z) in the complex plane.
%
%   CHECKERPLOT(f, [rmin rmax imin imax], __) uses the provided range. If
%      none is provided, [-1 1 -1 1] is used.
%
%   Further options:
%    - 'angle'  or  'a'
%         Show only contours of the angle of f(z).
%    - 'magnitude'  or  'm'
%         Show only contours of the logarithm of the magnitude of f(z).
%    - 'carthesisan'  or  'c'
%         Use carthesian instead of polar grid.
%    - 'real'  or  'r'
%         Show only contours of the real part of f(z).
%    - 'imaginary'  or  'i'
%         Show only contours of the imaginary part of f(z).
%    - 'thin'  or  't'
%         Use thin black lines insead of equal width.
%    - 'pixels', [N M]  or  'p', [N M]
%         Sample N x M pixels, default [720 720].

%  Copyright (c) 2023, Evert Provoost,
%     Department of Computer Science, KU Leuven.

%% Parse input
% set defaults
ax = [-1 1 -1 1]; % default axis
pxs = [720 720];      % default number of pixels
preal = 0;            % don't show real part by default
pimag = 0;            % don't show imaginary part by default
pangl = 1;            % show angle by default
pmagn = 1;            % show magnitude by default
thin = 0;             % don't use thin lines by default

setax = 0;
if nargin < 1
    error('CHECKERPLOT:inputs', 'No function f(z) provided.')
end
j = 1;
while j < nargin
    j = j+1;
    v = varargin{j-1};
    % user specified axis
    if ~ischar(v) && ~isstring(v)
        if setax
            warning('CHECKERPLOT:inputs', 'Axes already provided.')
        end
        ax = v;
        if length(ax) == 1
            ax = ax * [-1 1 -1 1];
        end
        setax = 1;
    % user wants carthesisan grid
    elseif isequal(v, 'carthesisan') || isequal(v, 'c')
        preal = 1; pimag = 1; pmagn = 0; pangl = 0;
    % user wants only real part
    elseif isequal(v, 'real') || isequal(v, 'r')
        preal = 1; pimag = 0; pmagn = 0; pangl = 0;
    % user wants only real part
    elseif isequal(v, 'imaginary') || isequal(v, 'i')
        preal = 0; pimag = 1; pmagn = 0; pangl = 0;
    % user wants only maginitude part
    elseif isequal(v, 'magnitude') || isequal(v, 'm')
        preal = 0; pimag = 0; pmagn = 1; pangl = 0;
    % user wants only angle part
    elseif isequal(v, 'angle') || isequal(v, 'a')
        preal = 0; pimag = 0; pmagn = 0; pangl = 1;
    % User wants thin lines
    elseif isequal(v, 'thin') || isequal(v, 't')
        thin = 1;
    % user specified number of pixels
    elseif isequal(v, 'pixels') || isequal(v, 'p')
        if j >= nargin || ~isnumeric(varargin{j})
            error('CHECKERPLOT:inputs', 'Number of pixels not provided.')
        end
        pxs = varargin{j};
        if length(pxs) == 1
            pxs = pxs * [1 1];
        end
        j = j+1;
    else
        error('CHECKERPLOT:inputs', ['Unrecognized input "',v,'".'])
    end
end

%% Checker plot
% compute f(z)
x = linspace(ax(1), ax(2), pxs(1));
y = linspace(ax(3), ax(4), pxs(2));
[xx,yy] = meshgrid(x, y);
zz = xx + 1i*yy;
try
    fzz = f(zz);
catch
    warning('CHECKERPLOT:inputs', ['f fails on array inputs, ' ...
            'plotting could be slow.'])
    fzz = arrayfun(f, zz);
end

C = ones(size(fzz));
% add requested parts
if preal
    C = C.*sin(5*pi*real(fzz));
end
if pimag
    C = C.*sin(5*pi*imag(fzz));
end
if pangl
    C = C.*sin(15*angle(fzz));
end
if pmagn
    C = C.*sin(15*log(abs(fzz)));
end

% what kind of grid?
if thin
    C = abs(C).^.08;
else
    C = .9*min(1,sign(C)+1)+.08;
end

% now show plot
C = repmat(C,[1 1 3]);
h = image([ax(1), ax(2)], [ax(3), ax(4)], C);
% image defaults to y inverted, fix that
set(gca, 'YDir', 'normal') 
% smooth image on newer versions
if isprop(h, 'Interpolation')
    set(h, 'Interpolation', 'bilinear')
end
% keep axes proportional
axis equal

% return image if asked
if nargout > 0
    varargout = {h};
end
