function [Xi,Yi,Zi] = polarplot3d(Zp,varargin)
% POLARPLOT3D  Plot a 3D surface from polar coordinate data
%   [Xi,Yi,Zi] = polarplot3d(Zp,varargin)
% Input
%   Zp         two dimensional matrix of input magnitudes or intensities
%              each column of Zp contains information along a single meridian
%              each row    of Zp contains information at    a single radius
%              the direction sense of the rows and columns is determined by
%              the relative order of the angular and radial vectors
%              by default Zp is increasing in radius down each column and
%              increasing in angle (counter-clockwise) along each row
%              the default plot is a full disc surface plot with unit radius
%
%   varargin   'property','value' pairs that modify plot characteristics
%              property names must be specified as quoted strings
%              property names and their values are not case sensitive
%              for each property the default value is given below
%
%   'AngularRange'  scalar or 2 element vector, radians     (default = [0 2*pi])
%   'RadialRange'   scalar or 2 element vector              (default = [0 1])
%                   if a scalar is given for either range it is used as
%                   the maximum value and zero is used for the minimum
%   'ColorData'     a matrix of color values, size must be equal to size(Zp)
%                   the default coloring is according to Zp magnitude
%   'MeshScale'     mesh element scale factor               (default = 1)
%                   values > 1 increase mesh element size   (coarser mesh)
%                   values < 1 decrease mesh element size   (finer   mesh)
%   'TickSpacing'   spacing of polar tick marks in degrees  (default = 10)
%                   every other tick mark is labeled
%   'PolarGrid'     polar grid density vector (2 elements)  (default = [10 36])
%                   number of grid lines in radial and azimuthal directions
%   'ContourLines'  scalar (number of contours) or a vector (contour locations)
%                   the default is auto selection by the contour function
%   'PlotType'      'surfn'   surface plot,               no mesh (default)
%                   'surfcn'  surface plot with contours, no mesh
%                   'surfc'   surface plot with contours
%                   'surf'    surface plot
%                   'mesh'    mesh    plot
%                   'meshc'   mesh    plot with contours
%                   'meshl'   mesh    plot with no interpolation
%                             full data matrix plotted, 'MeshScale' is ignored
%                   'contour' 2D contour plot
%                             use 'ContourLines' to set the number of contours
%                   'wire'    wireframe polar grid plot only, no surface plot
%                             a wireframe and polar axis overlay can be added to
%                             an existing plot (the current plot is not cleared)
%                   'off'     no plot is produced, interpolation only
%   'AxisLocation'  'surf'    polar axis along edge of surface     (default)
%                   'min'     polar axis at minimum Zp for largest radius
%                   'max'     polar axis at maximum Zp for largest radius
%                   'mean'    polar axis at mean    Zp for largest radius
%                   'top'     polar axis at top     of plot box
%                   'bottom'  polar axis at bottom  of plot box
%                   value     polar axis is drawn at specified height
%                   'off'     no polar axis
%   'InterpMethod'  'cubic'   bicubic          interpolation on Zp (default)
%                   'linear'  bilinear         interpolation on Zp
%                   'spline'  spline           interpolation on Zp
%                   'nearest' nearest neighbor interpolation on Zp
%   'PlotProps'     a cell array of 'property','value' pairs
%                   plotprops are passed to mesh, surf and contour
%                   the default value is an empty cell array
%
% Output
%   Xi,Yi,Zi   Cartesian locations corresponding to polar coordinates (T,R,Zp)
%              T and R are created from Trange and Rrange using meshgrid and
%              converted to Cartesian coordinates with pol2cart
%              Xi,Yi,Zi are square matrices with size equal to the maximum
%              dimension of Zp after interpolation
%              matrix sizes are reduced or enlarged by mscale
%              if 'meshl' is specified, size(Xi)=size(Yi)=size(Zi)=size(Zp)
%
% Notes        Zp is the only required input argument
%              if no input arguments are given an example plot is produced
%              and this help text is displayed in the command window
%
% Example
%   1) a transformation of data generated by the peaks function is plotted
%      in 3d polar coordinates between 30 and 270 degrees and from radius
%      80 to 400
%   2) a surface colored according to height above the xy plane is shown
%      without Cartesian grid lines
%   3) a 4 by 24 polar grid is drawn on top of the surface
%   4) a polar axis is drawn along the perimeter of the surface
%   5) tick marks are drawn at 8 degree intervals with labels every 16 degrees
%
%   P = peaks(101);                           % generate interesting data
%   P = (P+rot90(P)-rot90(3*fliplr(P)));
%   P(62:end,:) = [];
%
%   figure;                                   % draw 3d polar plot
%   polarplot3d(P,'PlotType','surfn','PolarGrid',[4 24],'TickSpacing',8,...
%                 'AngularRange',[30 270]*pi/180,'RadialRange',[80 400]);
%
%                                             % set plot attributes
%   set(gca,'DataAspectRatio',[10 10 1],'View',[-64,26],...
%           'Xlim',[-450 450],'Xtick',[-400 -200 0 200 400],...
%           'Ylim',[-450 450],'Ytick',[-400 -200 0 200 400]);
%   title('polarplot3d example');
%
% Ken Garrard, North Carolina State University, 2006-2007
% version 1    Original function based on POLAR3D by J De Freitas
%         2    Changed argument method to 'Property',value pairs
%         2.1  Added 'ColorData' property

% argument string values for PlotType, AxisLocation and InterpMethod properties
plst = {'mesh','meshc','meshl','wire',...
        'surf','surfc','surfn','surfcn','contour','off'};
alst = {'off','const','min','max','mean','surf','top','bottom'};
mlst = {'cubic','linear','spline','nearest'};

% if function was called without arguments
% draw example plot in a new figure and display help text
if nargin < 1
   P = peaks(101);                             % generate interesting data
   P = (P+rot90(P)-rot90(3*fliplr(P)));
   P(62:end,:) = [];

   figure;                                     % draw 3d polar plot
   polarplot3d(P,'PlotType','surfn','PolarGrid',[4 24],'TickSpacing',8,...
                 'AngularRange',[30 270]*pi/180,'RadialRange',[80 400]);

                                               % set plot attributes
   set(gca,'DataAspectRatio',[10 10 1],'View',[-64,26],...
           'Xlim',[-450 450],'Xtick',[-400 -200 0 200 400],...
           'Ylim',[-450 450],'Ytick',[-400 -200 0 200 400]);
   title('polarplot3d example');
                                               % display help text
   error(['No input arguments given\n'...
          'Please consult the help text and the example plot\n'...
          '--------\n%s'],help(mfilename));
end

%-- parse and validate input arguments

% setup property structure with default values
p.angularrange = [0 2*pi];    % angular range
p.radialrange  = [0 1];       % radial  range
p.plottype     = 'surfn';     % surface plot, no rectangular grid
p.meshscale    = 1;           % no mesh scaling
p.polargrid    = [10 36];     % number of azimuthal and radial grid lines
p.tickspacing  = 10;          % ten degree tick mark spacing
p.axislocation = '';          % handle as special case after parsing varargs
p.interpmethod = 'cubic';     % bicubic interpolation
p.colordata    = [];          % default coloring according to Zp values
p.contourlines = '';          % default contour specification
p.plotprops    = {};          % no additional plot properties

% parse property value pairs, replacing defaults with those specified by caller
try
   p = parse_pv_pairs(p,varargin);
catch
   error('Error parsing varargin list\n??? %s',lasterr);
end

% check plot type specification
if ~ischar(p.plottype) || isempty(strmatch(lower(p.plottype),plst,'exact'))
   error('Invalid ''PlotType'' property value');
end
p.plottype = lower(p.plottype);

% choose default polar axis location
if isequal(p.axislocation,'')            
   if ~isempty(strmatch(p.plottype,{'meshc','surfc','surfcn'},'exact'))
        p.axislocation = 'bottom';       % plot box bottom for contour plots
   else p.axislocation = 'surf';         % along perimeter of surface otherwise
   end
end

% user specified polar axis location as a numeric value
if isnumeric(p.axislocation)
   [r,c] = size(p.axislocation);
   if (((r ~= 1) || (c ~= 1)) || ~isreal(p.axislocation))
      error('''AxisLocation'' property value must be scalar, positive and real');
   end
   polax          = p.axislocation;
   p.axislocation = 'const';
end

% check polar axis location specification
if ~ischar(p.axislocation) || isempty(strmatch(lower(p.axislocation),alst,'exact'))
   error('Invalid ''AxisLocation'' property value');
end
p.axislocation = lower(p.axislocation);

% check angular and radial range vectors
p.angularrange = p.angularrange(:);
if isscalar(p.angularrange), p.angularrange = [0; p.angularrange]; end
if length(p.angularrange) ~= 2 || ~isreal(p.angularrange)
   error(['''AngularRange'' property value must be a scalar'...
            ' or a two element numeric vector']);
end

p.rrange = p.radialrange(:);
if isscalar(p.radialrange), p.radialrange = [0; p.radialrange]; end
if length(p.radialrange) ~= 2 || ~isreal(p.radialrange)
   error(['''RadialRange'' property value must be a scalar'...
            ' or a two element numeric vector']);
end

% angular and radial range vectors define data order
% reflect Zp left-right and/or up-down depending on angular and radial ranges
if p.angularrange(1) > p.angularrange(2), Zp = fliplr(Zp); end
if p.radialrange (1) > p.radialrange (2), Zp = flipud(Zp); end

[Tmin,Tmax] = deal(min(p.angularrange),max(p.angularrange));   % Tmin < Tmax
[Rmin,Rmax] = deal(min(p.radialrange ),max(p.radialrange ));   % Rmin < Rmax

% angular range cannot be more than one full circumference
if abs(Tmax - Tmin) > 2*pi
   error('Angular range cannot be greater than 2*pi');
end

% extract radial and azimuthal grid density from parameter structure
p.polargrid = p.polargrid(:);
if length(p.polargrid) ~= 2 || ~isreal(p.polargrid)
   error('''PolarGrid'' property value must be a two element numeric vector');
end
[radgrid,azmgrid] = deal(p.polargrid(1),p.polargrid(2));

% check azimuthal and radial grid density
if ~isnumeric(radgrid) || radgrid < 0
   error('Non-numeric or non-positive radial grid parameter');
end
if ~isnumeric(azmgrid) || azmgrid < 0
   error('Non-numeric or non-positive azimuthal grid parameter');
end

% check mesh scale factor property value
[r,c] = size(p.meshscale);
if (((r ~= 1) || (c ~= 1)) || ~isreal(p.meshscale)) || p.meshscale <= 0
   error('''MeshScale'' property value must be scalar, positive and real');
end

% check tick spacing property value
[r,c] = size(p.tickspacing);
if (((r ~= 1) || (c ~= 1)) || ~isreal(p.tickspacing)) || p.tickspacing <= 0
   error('''TickSpacing'' property value must be scalar, positive and real');
end

% check contour lines property value
if ~isnumeric(p.contourlines) && ~isequal(p.contourlines,'');
   error('''ContourLines'' property value must be numeric');
end
p.contourlines = p.contourlines(:);

% check interpolation method
p.interpmethod = lower(p.interpmethod);
if ~ischar(p.interpmethod) || isempty(strmatch(p.interpmethod,mlst,'exact'))
   error('Invalid ''InterpMethod'' property value');
end

% check input data size
[r,c] = size(Zp);
if (r < 5) || (c < 5)
   error('Zp matrix size must be greater than 4 x 4');
end

% check if mesh scale factor is compatible with input data dimension
if round(max(r,c)/p.meshscale) < 4
   error('Mesh scale factor is too large, not enough data remaining to plot');
end

% choose default color matrix
if isempty(p.colordata), p.colordata = Zp; end       % surface height coloring

% check color matrix
if ~isnumeric(p.colordata) || ~isequal(size(p.colordata),[r,c])
   error('Color matrix must be numeric and the same size as Zp');
end

% no polar axis or grid is generated if no plot was requested
% or if the plot type is a 2D contour plot
if ~isempty(strmatch(p.plottype,{'off','contour'},'exact'))
   p.axislocation  = 'off';
   azmgrid = 0;
   radgrid = 0;
end

%-- create polar grid and interpolate data

% create radius and angle vectors and polar grid for input data matrix
rho  = linspace(Rmin,Rmax,r);                        % radius vector
angl = linspace(Tmin,Tmax,c);                        % angle  vector
[xx,yy] = meshgrid(angl,rho);                        % mesh   matrices

% check for 'meshl' plot
if isequal(p.plottype,'meshl')                       % no interpolation
   Xi = rho' * cos(angl);                            % matrix of x's
   Yi = rho' * sin(angl);                            % matrix of y's
   Zi = Zp;                                          % z's == input data

% interpolate Zp to a scaled polar grid
else
   q    = fix(max(r,c)/p.meshscale);                 % new mesh size
   rho  = linspace(Rmin,Rmax,q);                     %     radius vector
   angl = linspace(Tmin,Tmax,q);                     %     angle  vector
   [theta,rad] = meshgrid(angl,rho);                 % create polar grid
   T = interp2(xx,yy,Zp,theta,rad,p.interpmethod);   % interpolate to grid
   [Xi,Yi,Zi] = pol2cart(theta,rad,T);               % convert     to Cartesian
end

% interpolate color data to same grid as height data
if isequal(p.plottype,'meshl')                       % no interpolation needed
     Ci = p.colordata;                               %    for 'meshl' plots
else Ci = interp2(xx,yy,p.colordata,theta,rad,p.interpmethod);
end

%-- plot the Cartesian surface

switch p.plottype
   case 'wire',    grid on;
   case 'meshc',   meshc(Xi,Yi,Zi,Ci);
   case 'mesh',    mesh (Xi,Yi,Zi,Ci,p.plotprops{:});
   case 'meshl',   mesh (Xi,Yi,Zi,Ci,p.plotprops{:});
   case 'surf',    surf (Xi,Yi,Zi,Ci,p.plotprops{:});
   case 'surfc',   surfc(Xi,Yi,Zi,Ci,p.plotprops{:});
   case 'surfn',   surf (Xi,Yi,Zi,Ci,'LineStyle','none',p.plotprops{:});
   case 'surfcn',  surfc(Xi,Yi,Zi,Ci,'LineStyle','none',p.plotprops{:});
   case 'contour', axis equal; polar([Tmin Tmax],[Rmin Rmax]); hold on;
                   contour(Xi,Yi,Zi,p.contourlines); colorbar; hold off;
end

%-- plot the polar axis

% axis and tick font attributes
fontargs = {'FontName','Arial','FontSize',10,'FontWeight','bold'};

if ~isequal(p.axislocation,'off')

   % create polar axis data just outside the largest radius
   xa = Rmax * 1.005 .* cos(angl);
   ya = Rmax * 1.005 .* sin(angl);

   % create polar axis tic marks at p.tickspacing intervals
   ts = 180/p.tickspacing;
   ta = pi/ts * (round(Tmin*ts/pi):1:round(Tmax*ts/pi));
   tr = Rmax  * [1.005; 1.03; 1.1];
   xt = tr * cos(ta);
   yt = tr * sin(ta);

   % set polax to z location of polar axis
   switch p.axislocation
      case 'min',    polax = min (Zi(end,:));
      case 'max',    polax = max (Zi(end,:));
      case 'mean',   polax = mean(Zi(end,:));
      case 'top',    zlim  = get(gca,'zlim'); polax = zlim(2);
      case 'bottom', zlim  = get(gca,'zlim'); polax = zlim(1);
   end

   % z values for polar axis and tick marks
   if isequal(p.axislocation,'surf'),        % vary along edge of surface
      za = Zi(end,:);
      zt = interp1(angl,za,ta,'linear');
   else
      za = zeros(size(xa))     + polax;      % constant location
      zt = zeros(1,length(ta)) + polax;
   end

   % label every other tick mark
   nl = round(length(ta)/2);

   % beginning and end of a full polar axis are identical, label once
   if abs(Tmin-Tmax) == 2*pi, nl = nl-1; end

   % plot the polar axis
   hold on;
   line(xa,ya,za,'Color','k','LineWidth',1);
   line(xt(1:2,:),yt(1:2,:),[zt; zt],'Color','k');   % add tick marks

   % draw tick labels
   for k = 2 * (1:nl) - 1
      text(xt(3,k),yt(3,k),zt(k),num2str(ta(k)*180/pi),...
           'HorizontalAlignment','Center',fontargs{:});
   end
   hold off;
end

%-- draw polar grid lines

if azmgrid > 0 || radgrid > 0

   % interpolate data to a fine mesh to produce smooth grid lines
   N = 360;                                          % new grid size
   m = abs(fix(N/max(1,azmgrid)));                   % angular grid spacing
   r = abs(fix(N/max(1,radgrid+1)));                 % radial  grid spacing
   rho  = linspace(Rmin,Rmax,N+1);                   % radius  vector
   angl = linspace(Tmin,Tmax,N+1);                   % angle   vector
   [theta,rad] = meshgrid(angl,rho);                 % create polar grid
   T = interp2(xx,yy,Zp,theta,rad,p.interpmethod);   % interpolate to grid
   [xi,yi,zi] = pol2cart(theta,rad,T);               % convert     to Cartesian
   hold on;                                          % plot in same figure

   % don't draw polar grid lines at minimum radius unless the plot
   % region is an annulus or only meridian grid lines are requested
   k = r+1;
   if Rmin > 0 || radgrid == 0, k = 1; end
   
   % draw radial grid lines in azimuthal direction (meridians)
   if azmgrid > 0
      xr = xi(k:end,1:m:end);
      yr = yi(k:end,1:m:end);
      zr = zi(k:end,1:m:end);
      plot3(xr,yr,zr,'Color','k','LineWidth',1);     % plot grid lines
   end

   % draw azimuthal grid lines in radial direction (concentric arcs)
   if radgrid > 0
      xm = xi(k:r:end-1,:);
      ym = yi(k:r:end-1,:);
      zm = zi(k:r:end-1,:);
      plot3(xm',ym',zm','Color','k','LineWidth',1);  % plot grid lines
   end
   hold off;
end

%-- set axis font
if ~isequal(p.plottype,'off'), set(gca,fontargs{:}); end
