```function PP = pltphase(z,f,cs,t,pres,height)
% phase plot of complex function f(z), internal use

% Usage: PP = pltphase(z,f,cs,t,pres,height)

% Part of the Complex Function Explorer (former Phase Plot) package
% Version 1.1, February 1, 2014
% Copyright(c) 2012-2014 by Elias Wegert (elias@wegert.com, www.wegert.com)

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if nargin<6, height=0; end

if nargin<5 || isempty(pres), pres=20; end

if nargin<4 || isempty(t), t=exp(1i*linspace(-pi,pi,pres+1)); t=t(1:pres); end

if nargin<3 || isempty(cs), cs = 'p'; end

if nargin<2, z=zdomain; f=sin(5*z); end

% convert scalar to a constant function
if length(f)==1, f=f*ones(size(z)); end

if nargin<6
% if height is not given, choose shape of surface on which the plot lives

%rep='a';  standard analytic landscape
%rep='l';  logarithmic analytic landscape
%rep='c';  compressed analytic landscape

rep='0';  height=zeros(size(z)); hmin=0; hmax=1; %#ok<NASGU>

if length(cs)>=2

% modify if color scheme requests special surface
incsa = strfind(cs(2:end),'a');
if ~isempty(incsa)>0, rep=cs(incsa(1)+1); cs(incsa+1) = ''; end

incsl = strfind(cs(2:end),'l');
if ~isempty(incsl)>0, rep=cs(incsl(1)+1); cs(incsl+1) = ''; end

incsc = strfind(cs(2:end),'c');
if ~isempty(incsc)>0, rep=cs(incsc(1)+1); cs(incsc+1) = ''; end

if rep=='a'
% analytic landscape, cutoff at hmax
hmax=5;
height=abs(f);
height(height>hmax)=hmax;
% normalize to [0,1]
hmin = min(min(height)); hmax = max(max(height));
if hmax>hmin
height = (1/(hmax-hmin))*(height-hmin);
end

elseif rep=='l'
% logarithmic analytic landscape
hmin =-3; hmax=3;
absf = abs(f);
height = log(absf);
small = (height<hmin);
large = (height>hmax);
height(small) = hmin;
height(large) = hmax;
% normalize height to [0,1]
hmin = min(min(height)); hmax = max(max(height));
if hmax>hmin
height = (1/(hmax-hmin))*(height-hmin);
end

elseif rep=='c'
% compressed analytic landscape, zeros and poles at levels 0 and 1
height=(2/pi)*atan(abs(f));

end

end

end

% flip arrays for correct orientation of surface normal
z=flipud(z);
f=flipud(f);
height=flipud(height);

hmin = min(min(height));
hmax = max(max(height));

% scale height for nice shape of surface
xmin=min(min(real(z)));
xmax=max(max(real(z)));
ymin=min(min(imag(z)));
ymax=max(max(imag(z)));

scale=0.4*(xmax-xmin+ymax-ymin);
height=scale*height;

% get RGB colors according to chosen color scheme
RGB  = colscheme(f,cs,t,pres);

% plot surface
PP = surf(real(z),imag(z),height,RGB);
set(PP,'EdgeColor','none');

% set appropriate axis
hmin=scale*hmin;
hmax=scale*hmax;

axis equal
axis tight
ax = axis;
if hmax>hmin, axis([ax,hmin,hmax]); end
axis off

```