classdef polariCommon
    % Package of helper-functions for polari plot.
     
    % Listener support across many classes (not just polari)
    %
    methods (Static)
        function Lvec = deleteListenerStruct(Lvec)
            % L is a structure with listener handles retained in each
            % field.  Identifies each field, deletes each listener, and
            % replaces the handle with an empty in the struct.
            %
            % If any field contains a struct, recursive descent into that
            % field is performed.
            %
            % L may be an empty matrix or struct.
            % Returns the updated struct.
            if ~isempty(Lvec)
                validateattributes(Lvec,{'struct'}, {}, 'deleteListenerStruct');
            end
            NLvec = numel(Lvec);
            for j = 1:NLvec
                L = Lvec(j);
                f = fieldnames(L);
                for i = 1:numel(f)
                    f_i = f{i};   % get next field name
                    v = L.(f_i);  % get value from field
                    if isstruct(v)
                        % Recurse on child struct:
                        L.(f_i) = internal.polariCommon.deleteListenerStruct(v);
                    else
                        % Delete listener handle
                        delete(v);
                    end
                    L.(f_i) = []; % replace value in field
                end
                Lvec(j) = L;
            end
        end
    end
    
    % Used in user-clickable links shown in the command window.
    %
    % Scope of invocation is base workspace: can't put these in private.
    %
    methods (Static)
        function methodsForDisplay(varname, varargin)
            if nargin == 1
                s = getString(message('MATLAB:graphicsDisplayText:FooterLinkFailureMissingVariable',varname));
                disp(s)
            else
                obj = varargin{1};
                classname = varargin{2};
                if ~isa(obj,classname)
                    dots = strfind(classname,'.');
                    if ~isempty(dots)
                        classname = classname(dots(end)+1:end);
                    end
                    s = getString(message('MATLAB:graphicsDisplayText:FooterLinkFailureClassMismatch',varname, classname));
                    disp(s)
                else
                    try
                        methods(obj)
                    catch
                        s = getString(message('MATLAB:graphicsDisplayText:FooterLinkFailureUnknown',varname));
                        disp(s)
                    end
                end
            end
        end
        
        function getForDisplay(varname, varargin)
            % Helper function for object display
            
            if nargin == 1
                s = getString(message('MATLAB:graphicsDisplayText:FooterLinkFailureMissingVariable',varname));
                disp(s)
            else
                obj = varargin{1};
                classname = varargin{2};
                if ~isa(obj,classname)
                    dots = strfind(classname,'.');
                    if ~isempty(dots)
                        classname = classname(dots(end)+1:end);
                    end
                    s = getString(message('MATLAB:graphicsDisplayText:FooterLinkFailureClassMismatch',varname, classname));
                    disp(s)
                else
                    try
                        %get(obj)
                        showAllProperties(obj)
                    catch
                        s = getString(message('MATLAB:graphicsDisplayText:FooterLinkFailureUnknown',varname));
                        disp(s)
                    end
                end
            end
        end
    end
    
    % Numerics on angles
    %
    methods (Static)
        function cw = isCW(c1,c2)
            %isCW True if two complex points are in clockwise order.
            %   isCW(A,B) returns true if the complex points A and B are in
            %   clockwise order when the angle traversed between A and B is
            %   180 degrees or less.
            %
            %   Both A and B may be scalars, both the same length vectors,
            %   or one may be a scalar and the other a vector.
            
            % compute z-coord of cross-product,
            %     c = cross(a,b)
            % But here, a and b are complex, not 2-element vectors.
            % Also, we're only interested in the sign of c(3),
            %     c = a(1).*b(2)-a(2).*b(1)
            % If c is +ve, it's CW; otherwise it's CCW
            
            cw = imag(c1).*real(c2) - real(c1).*imag(c2) >= 0;
        end
        
        function cw = isCWangle(a1,a2)
            c1 = complex(cos(a1),sin(a1));
            c2 = complex(cos(a2),sin(a2));
            cw = internal.polariCommon.isCW(c1,c2);
        end
        
        function [i1,i2] = cswapIdxIfCW(c,i1,i2)
            % Swap indices i1 and i2 if they are in CW order.
            if internal.polariCommon.isCW(c(i1),c(i2))
                t = i2; i2 = i1; i1 = t;
            end
        end
        
        function [dCCW,dCW] = cangleDiff(A,B)  % cangleDiffOriented
            %cangleDiff Oriented angle difference for complex values.
            %   cangleDiff(A,B) returns counter-clockwise angle differences,
            %   in radians, from the angle of complex scalar A to the angle
            %   of each element of complex vector B.
            %
            %   [CCW,CW] = cangleDiff(A,B) returns both counter-clockwise
            %   and clockwise angle differences.
            %
            %   CCW = cangleDiff(A) returns counter-clockwise angle
            %   differences, in radians, from adjacent pairs of angles in
            %   complex vector A.
            %
            %   [CCW,CW] = cangleDiff(A) returns both counter-clockwise and
            %   clockwise angle differences.
            %
            %   See also angleDiff, angleAbsDiff, cangleAbsDiff.
            
            % ALGORITHM
            %
            % Rotate axes to align with angle of first scalar. Project
            % vector into that space, and compute 4-quadrant arctan.
            % Rotation requires no trig.
            
            % ALGORITHM
            %
            % Express each complex value b in vector B as [bx by 0], and A
            % as [ax ay 0].  The cross product AxB yields C,
            %     C = |A||B||sin(th)| u,
            % where u is a unit vector, C is [cx cy cz] and angle
            % measurement th <= pi.  Assuming A and B are on the unit
            % circle, C reduces to |sin(th)| u.
            %
            % cz gives the orientation of the angle measurement returned by
            % the cross product.  If cz > 0, angle measurement from A to B
            % is CCW, which follows from the right-hand rule: 0-deg is
            % toward +x and increasing angle is CCW. If we desire CCW sense
            % for the angle difference, we keep th when cz > 0, and we
            % record 2*pi-th when cz < 0.  If we desire CW sense, this
            % logic is reversed.
            %
            % With this procedure, we get oriented angle difference with a
            % cross product and a trig computation, and no angle
            % unwrapping is needed.
            %
            % |sin(th)| yields single-quadrant angle resolution, and we use
            % signs to transform results into 4-quadrant angle differences.
            
            twoPi = 2*pi;
            if nargin == 1
                % Single-input case:
                %
                % Use conjugate-delay algorithm which is equivalent to
                % rotation, but cheaper:
                dCCW = angle(A(2:end) .* conj(A(1:end-1)));
            else
                % Two-input case:
                %
                if ~isscalar(A)
                    error('Input A must be a scalar value.');
                end
                dCCW = angle(B .* conj(A));
            end
            sel = dCCW < 0;
            dCCW(sel) = dCCW(sel) + twoPi;
            if nargout > 1
                dCW = twoPi - dCCW;
            end
            
            %{
            % Theory for comparison:
                % Difference between adjacent pairs of angles in A
                %
                % Comments for the 2-input case apply here
                N = numel(A)-1;
                if size(A,1) > 1
                    dCCW = zeros(N,1);
                else
                    dCCW = zeros(1,N);
                end
                Ax = real(A);
                Ay = imag(A);
                for i = 1:N
                    % Project next point A(i+1) onto rotated set of axes
                    % aligned with A(i), then compute angle of A(i+1):
                    %
                    % Signal processing conj delay formulation:
                    % = A(i+1) * conj(A(i))
                    % =     Ax(i+1)*Ax(i) + Ay(i+1)*Ay(i)
                    %   j  -Ax(i+1)*Ay(i) + Ax(i)*Ay(i+1)
                    %
                    % Matrix rotation formulation:
                    % Bpq = [Ax(i+1) Ay(i+1)] * [Ax(i) -Ay(i); Ay(i) Ax(i)];
                    Bp =  Ax(i+1)*Ax(i) + Ay(i+1)*Ay(i);
                    Bq = -Ax(i+1)*Ay(i) + Ay(i+1)*Ax(i);
                    d_i = atan2(Bq,Bp);
                    if d_i < 0
                        dCCW(i) = d_i + twoPi;
                    else
                        dCCW(i) = d_i;
                    end
                end
            
            % Difference between scalar A and vector B
            %
            % Form orthonormal basis vectors on rotated coord system
            % P: perpendicular unit vector, CW to A
            %    P & Q form an orthonormal basis, rotated to A
            % Q: unit vector based on A
            %    (A is on unit circle, so no normalization needed)
            %
            % Taken in this order, we get a right-hand rule system.
            %
            % P=<Ax,Ay>, so Px=Ax and Py=Ay
            % Q=<Ay,-Ax>, so Qx=Ay and Qy=-Ax
            
            dCCW = zeros(size(B));
            if nargout > 1
                dCW = zeros(size(B));
            end
            Ax = trigFixup(real(A));
            Ay = trigFixup(imag(A));
            Bx = trigFixup(real(B));
            By = trigFixup(imag(B));
            
            for i = 1:numel(B)
                % Project B into the "PQ" rotated coord system, where P
                % is aligned along the A vector.
                %
                % P = < Ax,Ay>, so Px =  Ax and Py = Ay
                % Q = <-Ay,Ax>, so Qx = -Ay and Qy = Ax
                %
                % By definition, projecting A onto PQ axes yields
                % (p,q)=(1,0), i.e., one unit along P, no units along Q
                %
                % Bp = Bx*Px + By*Py ->  Bx*Ax + By*Ay
                % Bq = Bx*Qx + By*Qy -> -Bx*Ay + By*Ax
                Bp =  Bx(i)*Ax + By(i)*Ay;
                Bq = -Bx(i)*Ay + By(i)*Ax;
                d_i = atan2(Bq,Bp);
                if d_i < 0
                    dCCW(i) = d_i + twoPi;
                else
                    dCCW(i) = d_i;
                end
                if nargout > 1
                    if d_i < 0
                        dCW(i) = -d_i;
                    else
                        dCW(i) = twoPi - d_i;
                    end
                end
            end
            %}
        end
        
        function d = cangleDiffRel(A,B)
            %cangleDiffRel Relative angle difference for complex values.
            %   cangleDiffRel(A,B) returns relative angle differences, in
            %   radians, from the angle of complex scalar A to the angle of
            %   each element of complex vector B.  Angle differences are in
            %   the range -pi to +pi.
            %
            %   cangleDiffRel(A) returns relative angle difference between
            %   adjacent pairs of angles in complex vector A.
            
            if nargin==1
                % Single-input case:
                %
                % Use conjugate-delay algorithm which is equivalent to
                % rotation, but cheaper:
                d = angle(A(2:end) .* conj(A(1:end-1)));
            else
                % Two-input case:
                %
                if ~isscalar(A)
                    error('Input A must be a scalar value.');
                end
                d = angle(B .* conj(A));
            end
        end
        
        function d = cangleAbsDiff(A,B)
            %cangleAbsDiff Unoriented angle difference for complex values.
            %   cangleAbsDiff(A,B) returns absolute differences, in radians,
            %   from the angle of complex scalar A to the angle of each
            %   element of complex vector B.  Angle differences are in the
            %   range 0 to pi.
            %
            %   cangleAbsDiff(A) returns absolute differences between
            %   adjacent elements of complex vector A.
            %
            %   See also angleAbsDiff, angleDiff, cangleDiff.
            
            % Naive algorithm: return minimum angle for CW and CCW:
            %  [dCCW,dCW] = internal.polariCommon.cangleDiff(A,B);
            %  d = min(dCCW,dCW);
            %
            % Better algorithm: largely follows cangleDiff, but returns the
            % minimum angle for rotated quadrants II and III. This is
            % faster than naive algorithm, and uses less memory since naive
            % requires two return vectors and a min operation.
            
            if nargin == 1
                % Single-input case:
                %
                % Use conjugate-delay algorithm which is equivalent to
                % rotation, but cheaper:
                d = abs(angle(A(2:end) .* conj(A(1:end-1))));
            else
                % Two-input case:
                %
                if ~isscalar(A)
                    error('Input A must be a scalar value.');
                end
                d = abs(angle(B .* conj(A)));
            end
        end
        
        function varargout = angleDiff(angA,angB)
            %angleDiff Oriented angle difference.
            %   angleDiff(A,B) returns counter-clockwise angle differences,
            %   in radians, from scalar angle A to each element of angle
            %   vector B, in radians.
            %
            %   [CCW,CW] = angleDiff(A,B) returns both counter-clockwise
            %   and clockwise angle differences.
            %
            %   angleDiff(A) returns angle differences, in radians, between
            %   adjacent elements in vector A.
            %
            %   See also cangleDiff, angleAbsDiff, cangleAbsDiff.
            
            if nargin==1
                % Diff between adjacent elements
                [varargout{1:nargout}] = internal.polariCommon.cangleDiff(complex(cos(angA),sin(angA)));
            else
                % scalar angA to vector angB
                [varargout{1:nargout}] = internal.polariCommon.cangleDiff( ...
                    complex(cos(angA),sin(angA)), ...
                    complex(cos(angB),sin(angB)) );
            end
        end
        
        function d = angleAbsDiff(angA,angB)
            %angleAbsDiff Unoriented angle difference.
            %   angleAbsDiff(A,B) returns absolute differences, in radians,
            %   from the scalar angle A to each element of angle vector B.
            %   Angle differences are in the range 0 to pi.
            %
            %   angleAbsDiff(A) returns absolute differences between
            %   adjacent elements of angle vector A.
            %
            %   See also cangleAbsDiff, angleDiff, cangleDiff.
            
            if nargin==1
                % Diff between adjacent elements
                d = internal.polariCommon.cangleAbsDiff(complex(cos(angA),sin(angA)));
            else
                % scalar angA to vector angB
                d = internal.polariCommon.cangleAbsDiff( ...
                    complex(cos(angA),sin(angA)), ...
                    complex(cos(angB),sin(angB)) );
            end
        end
        
        function d = angleDiffRel(angA,angB)
            %angleDiffRel Relative angle difference.
            %   angleDiffRel(A,B) returns relative angle difference, in
            %   radians, from the scalar angle A to each element of angle
            %   vector B.  Angle differences are in the range -pi to +pi.
            %
            %   angleDiffRel(A) returns relative angle difference between
            %   adjacent pairs of angles in A.
            %
            %   See also cangleAbsDiff, angleDiff, cangleDiff.
            
            if nargin==1
                % Diff between adjacent elements
                d = internal.polariCommon.cangleDiffRel(complex(cos(angA),sin(angA)));
            else
                % scalar angA to vector angB
                d = internal.polariCommon.cangleDiffRel( ...
                    complex(cos(angA),sin(angA)), ...
                    complex(cos(angB),sin(angB)) );
            end
        end
        
        function sel = anglesWithinSpan(ang,span)
            % Returns true for angles in vector ANG that fall on or within
            % the angle-span SPAN defined by [minAng maxAng].  All angles
            % must be in radians.
            
            if numel(span) ~= 2
                error('SPAN must be a 2-element vector.');
            end
            s1 = span(1);
            s2 = span(2);
            sel = false(size(ang));
            TwoPi = 2*pi;
            for i = 1:numel(ang)
                d1 = internal.polariCommon.angleDiff(s1,ang(i));
                d2 = internal.polariCommon.angleDiff(ang(i),s2);
                sel(i) = d1+d2 <= TwoPi; % ang(i) is in interval
            end
        end
        
        function idx = spanContainingAngle(spans,ref)
            % Returns the row index of the Nx2 angle matrix SPANS that
            % includes angle REF.  All angles must be in radians.
            %
            % SPANS = [span1_startAngle span1_endAngle;
            %          span2_startAngle span2_endAngle;
            %          ... ]
            %
            % If REF falls exactly on a span endpoint and the endpoint is
            % repeated in the SPANS matrix, the lowest row index with a 
            % matching endpoint is returned.  If REF does not appear within
            % any span, or SPANS is empty, empty is returned.
            
            if ~ismatrix(spans) || (~isempty(spans) && size(spans,2)~=2)
                error('SPANS must be an Nx2 matrix.');
            end
            N = size(spans,1); % # spans = # index pairs
            for i = 1:N
                d1 = internal.polariCommon.angleDiff(spans(i,1),ref);
                d2 = internal.polariCommon.angleDiff(ref,spans(i,2));
                if d1+d2 <= 2*pi
                    idx = i; % ref is in interval
                    return % EARLY EXIT
                end
            end
            idx = []; % ref is not in any interval
        end
        
        function [mag,idx] = polarInterp(ang,angVec,magVec)
            % Linearly interpolate between closest two angle samples.
            % Assumes angle in degrees.
            
            % Now in radians:
            ang = internal.polariCommon.constrainAngleTwoPi(ang*pi/180);
            [angVecSort,sortIdx] = sort(internal.polariCommon.constrainAngleTwoPi(angVec*pi/180));
            
            % find closest entries
            adif = angVecSort-ang;
            i1 = find(adif<=0,1,'last');
            i2 = find(adif>=0,1,'first');
            if isempty(i1) || isempty(i2)
                % Need extrapolation
                mag = [];
                idx = [];
                return
            end
            if i1>=i2
                % i1==i2: Exact hit on sample
                % i1>i2: Multiple exact-hit samples
                % Use i2, the FIRST exact hit sample:
                mag = magVec(sortIdx(i2));
                idx = sortIdx(i2);
                return
            end
            a1 = angVecSort(i1);
            a2 = angVecSort(i2);
            if abs(a1-a2) <= eps
                % close enough to be identical
                mag = magVec(sortIdx(i1));
                idx = sortIdx(i1);
                return
            end
            frac = (ang-a1)/(a2-a1);
            assert(frac>0 && frac<1)
            mag = (1-frac)*magVec(sortIdx(i1)) + frac*magVec(sortIdx(i2));
            
            % Return interpolation index (fractional here).
            %
            % We guarantee that this isn't extrapolation, so no need to
            % check if idx exceeds num points.
            idx = sortIdx(i1)+frac;
        end
        
        %{
        function th = angleDiffDir(A,B,Dstr)
            %angleDiffDir Oriented angle difference.
            %   angleDistance(A,B) returns angle difference in radians from
            %   complex scalar A to each element of complex vector B, given
            %   a counter-clockwise measurement orientation in the plane.
            %
            %   Clockwise angle differences may be obtained by computing
            %   2*pi-angleDistance(A,B).
            %
            %   angleDistance(A,B,DIR) specifies angle measurement
            %   orientation to be either 'cw' for clockwise or 'ccw' for
            %   counter-clockwise.
            
            % ALGORITHM
            %
            % Express each complex value b in vector B as [bx by 0], and A
            % as [ax ay 0].  The cross product AxB yields C,
            %     C = |A||B||sin(th)| u,
            % where u is a unit vector, C is [cx cy cz] and angle
            % measurement th <= pi.  Assuming A and B on the unit circle, C
            % is |sin(th)| u.
            %
            % cz gives the orientation of the angle measurement returned by
            % the cross product.  If cz > 0, angle measurement from A to B
            % is CCW, which follows from the right-hand rule: 0-deg is
            % toward +x and increasing angle is CCW. If we desire CCW sense
            % for the angle difference, we keep th when cz > 0, and we
            % record 2*pi-th when cz < 0.  If we desire CW sense, this
            % logic is reversed.
            %
            % With this procedure, we get oriented angle difference with a
            % cross product and a trig computation, and no angle
            % unwrapping is needed.
            %
            % |sin(th)| yields single-quadrant angle resolution, and we use
            % sign to rotate result into a 4-quadrant angle difference.
            
            if ~isscalar(A)
                error('Input must be a scalar value.');
            end
            dirCCW = nargin<3 || strcmpi(Dstr,'ccw');
            sCW = (-1).^dirCCW; % -1 if CCW, +1 if CW
            Ax = real(A);
            Ay = imag(A);
            sAx = Ax < 0;
            sAy = Ay < 0;
            th = zeros(size(B));
            for i = 1:numel(B)
                Bx = real(B(i));
                By = imag(B(i));
                
                % cross product specialized for Az=0 and Bz=0:
                Cz = Ax.*By - Ay.*Bx;
                th_i = asin(abs(Cz)); % norm -> abs on vector
                
                % th_i is a one-quadrant difference, i.e., cross product
                % returns abs(sin(th)) for points on unit circle.  Use
                % sign differences to update for two-quadrant diff:
                if (Bx<0)~=sAx && (By<0)~=sAy
                    th_i = pi-th_i;
                end
                
                % if CCW, Cz >= 0 is natural and Cz < 0 needs correction.
                % if  CW, Cz <= 0 is natural and Cz > 0 needs correction.
                %
                % Straigtforward conditional:
                %    dirCCW && (Cz<0) || ~dirCCW && (Cz>0)
                % Reduced conditional:
                if sCW*Cz > 0
                    th_i = 2*pi-th_i;
                end
                th(i) = th_i;
            end
        end
        %}
        
        function ar = angleRecenter(a,b)
            %angleRecenter Shift domain of angle in radians.
            %   angleRecenter(A,B) shifts the radian angle A by
            %   2*pi*N, N=...,-2,-1,0,1,2,..., such that A lies in the
            %   domain [B-pi, B+pi].
            %
            %   This enables angle computations such as diff, without
            %   concern of "wrapped" angle math, i.e., angles being in
            %   different principal 2*pi domains.
            
            % Alt 1:
            %a_p = mod(a+pi,2*pi)-pi; % a is in range [-pi,+pi]
            %b_p = mod(b+pi,2*pi)-pi; % b is in range [-pi,+pi]
            %ar = a_p + b - b_p;
            % Alt 2:
            %ar = b + diff(mod([b a]+pi,pi+pi));
            % Alt 3:
            ar =  mod(a+pi,pi+pi) - mod(b+pi,pi+pi) + b;
        end
        
        function a = constrainAnglePi(a)
            a = mod(a+pi,2*pi)-pi;
        end
        
        function a = constrainAngleTwoPi(a)
            a = mod(a,2*pi);
        end
        
        
        % Map [-inf,inf] to [-360,360]
        % userDeg = userDeg-fix(userDeg/360)*360
        % Map [-360,360] to [-180,180]
        % userDeg = mod(userDeg+180,360)-180;
        % Map [-360,360] to [0,360]
        % userDeg = mod(userDeg,360);
        
        function [ar,f] = findAngleSpan(th,Nd)
            %findAngleSpan Span of angles represented by angle vector.
            %   S = findAngleSpan(A) returns the contiguous span of angles
            %   represented in angle vector A as a two-element vector S =
            %   [Amin Amax].
            %
            %   The span is not simply the minimum and maximum of the
            %   vector A in general, since angles in A may have wrapped
            %   around the unit circle, and the data set may represent a
            %   contiguous angular span less than a full circle. Angles may
            %   be nonuniformly spaced, and must be in radians.
            %
            %   S = findAngleSpan(A,Nd) sets the minimum angle span
            %   recognized as a contiguous region with no data points to be
            %   Nd*min(delta), where min(delta) is the smallest difference
            %   between successive data points.  The span returned is
            %   reduced by the largest region found with no data points.
            %   If omitted, Nd = 2.
            %
            %   [S,F] = findAngleSpan(...) returns logical flag F that is
            %   true when the span is the full circle.
            
            if nargin < 2
                Nd = 2;
            end
            if isempty(th)
                ar = [];
                f = false;
                return
            end
            
            % To find deltas around first/last bin, repeat first value at
            % end and last value at start.  But only repeat the value if
            % they are not the same.
            %
            %th = th(:);
            %th = [th(end); th; th(1)];
            %
            % watch for the repeated sample later
            %
            th2 = internal.polariCommon.constrainAngleTwoPi(th); % takes/returns radians
            if abs(th2(1) - th2(end)) > 0.001*pi/180
                th2 = [th2(:); th2(1)];
            end
            
            d = internal.polariCommon.angleAbsDiff(th2);
            minClearSpan = Nd * min(d); % min span to declare a gap
            [maxClearSpan,exclIdx] = max(d); % biggest span found
            if maxClearSpan >= minClearSpan
                % We have a gap. We react to one gap - the largest one.
                % Set inclIdx to a 1x2 vector of indices: [i1 i2],
                % where a gap is implied after i2.  Circular continuity is
                % assumed.
                N = numel(d);
                % N = numel(B)-1, but B has 1 extra sample
                if exclIdx==N
                    inclIdx = [1 N];
                elseif exclIdx==1
                    inclIdx = [2 1];
                else
                    inclIdx = [exclIdx+1 exclIdx];
                end
                ar = th2(inclIdx);
                f = false;
            else
                % No gaps
                ar = th2([1 1]);
                f = true;
            end
        end
        
        function y = isBetweenAngles(a,s1,s2)
            %isBetweenAngles True if angle lies in/on CW angle span.
            %  isBetweenAngles(A,S1,S2) returns true if angle A is on or
            %  between angles S1 and S2 interpreted as an inclusive span in
            %  the increasing angle direction (including spans exceeding
            %  180 degrees), where all angles are in DEGREES.
            
            y = internal.polariCommon.isBetweenAnglesRad(a*pi/180,s1*pi/180,s2*pi/180);
        end
        
        function y = isBetweenAnglesRad(a,s1,s2)
            %isBetweenAnglesRad True if angle lies in/on angle span.
            %  isBetweenAnglesRad(A,SPAN) returns true if angle A is on or
            %  between angles SPAN(1) and SPAN(2) interpreted as an
            %  inclusive span in the increasing angle direction (including
            %  wrapping past 2*pi radians), where all angles are in
            %  RADIANS.
            
            % Algorithm:
            % - Add 2pi to second span element until it is no more than
            %   2pi higher than first span element
            % - Add 2pi to a until it is no more than 2pi greater than
            %   first span element
            % - return true if tspan(1) <= ta <= tspan(2)
            TwoPi = 2*pi;
            if s2 > s1
                % see if we need to reduce s2 to be within 2pi > s1
                s2 = s2 - floor((s2-s1)/TwoPi) * TwoPi;
            elseif s1 > s2
                s2 = s2 + ceil((s1-s2)/TwoPi) * TwoPi;
            end
            
            % If s1=0 and s2=360, we will have s1==s2 now.
            % If so, assume full 2pi and return true:
            if abs(s1-s2) < 100*eps
                % close enough to equal:
                y = true;
            else
                if a > s1
                    % see if we need to reduce a to be within 2pi > s1
                    a = a - floor((a-s1)/TwoPi) * TwoPi;
                elseif s1 > a
                    a = a + ceil((s1-a)/TwoPi) * TwoPi;
                end
                y = a >= s1 && a <= s2;
            end
        end
        
        function th = linspaceIncrRad(s1,s2,N)
            % Generate radian angles from s1 to s2 by increasing angle.
            % Angles are interpreted modulo 2*pi.
            
            s1 = mod(s1,2*pi);
            s2 = mod(s2,2*pi);
            if s1 < s2
                th = linspace(s1,s2,N);
            else % s1 > s2
                th = linspace(s1,s2+2*pi,N);
            end
        end
    end
    
    % Peaks and Valleys
    %
    methods (Static)
        function valleysIdx = findPolarValleys(data,peaksIdx,isCircle,units,maxHeight)
            %Return indices of valleys in data defined on a circle.
            %   findPolarValleys(Y,PeakIndices,MaxHeight) returns data
            %   indices of at most one valley between each peak index
            %   specified in vector PeakIndices.  Valleys at or above
            %   MaxHeight are rejected; if omitted, MaxHeight = inf. Data
            %   must be a vector.
            %   If MaxHeight is omitted we use 3db less than max peak
            %   height.
            
            if nargin < 3
                isCircle = false;
            end
            if nargin < 5
                maxHeight = inf;
            end
            cnt = 0;
            valleysIdx = zeros(size(data));
            Np = numel(peaksIdx);
            Nd = numel(data);
            peaksIdx = sort(peaksIdx);
            for i = 1:Np
                if i < Np
                    iR = peaksIdx(i) : peaksIdx(i+1);
                    if nargin < 5 && strcmp(units,'dB')
                        maxHeight = max(data([peaksIdx(i),  peaksIdx(i+1)]))-3;
                    end
                else
                    if ~isCircle
                        break
                    end
                    iR = [peaksIdx(i):Nd, 1:peaksIdx(1)];
                    if nargin < 5 && strcmp(units,'dB')
                        maxHeight = max(data([peaksIdx(i),  peaksIdx(1)]))-3;
                    end
                end
                [~,vIdx] = internal.antfindpeaks(-data(iR), ...
                    'NPeaks',1, ...
                    'SortStr','descend', ...
                    'MinPeakHeight',-maxHeight);
                if isempty(vIdx)
                    [valValley, vIdx] = max(-data(iR));
                    if -valValley > maxHeight && strcmp(units,'dB')
                        continue;
                    end
                end
                cnt = cnt+1;
                valleysIdx(cnt) = iR(vIdx);
            end
            valleysIdx = valleysIdx(1:cnt);
        end
        
        function peakIdxVec = getPeakMidpoints(mag,peakIdxVec)
            % Updates the index of the peak data points to the true
            % "midpoints" of each peak, if there are runs of
            % identical-valued maxima after the peak locations.
            %
            % pidxVec is a vector of peak locations found by findpeaks.
            % Preserves the shape (row or col) of the peakIdxVec vector.
            %
            % In the presence of multiple identical contiguous peak values,
            % findpeaks returns the index of the first maxima.  However, we
            % sometimes require the "middle" value of such spans. For
            % example, the midpoint index best approximates the midpoint of
            % an "antenna lobe" pattern in polar data.
            %
            % The length of the identical-valued data may be the entire
            % dataset, and that could take time to evaluate.  So we do a
            % first test to see if the very NEXT index (in wrapped index
            % space) is identical in value.  We pursue a more complete test
            % only if this necessary precursor test passes.
            %
            % When useful, ang is used to choose the "best possible" peak
            % position that aligns with a 90-degree cardinal point.  Only
            % done when the choice between two adjacent points is
            % equivalent and one has a "nicer angle."
            
            thresh = 10*eps; % value differences <= thresh are "identical"
            
            Nr = numel(mag); 
            Ni = numel(peakIdxVec);  % a vector of peak indices
            for i = 1:Ni
                pidx = peakIdxVec(i);
                pval = mag(pidx);
                
                pnextIdx = 1 + rem(pidx,Nr); % +1 index offset implicit
                if abs(mag(pnextIdx) - pval) <= thresh
                    % Test data for contiguous span of max value.
                    % rotated index 1 is a global min, so we don't bother
                    % wrapping the test.  Just test from pidx to the end.
                    % That's overkill already.
                    %
                    % Since we now know there are at least TWO identical max
                    % values, we search from the first (pidx) to the end of the
                    % rotated data, and identify the first index where this is
                    % NOT true.
                    p_rel_idx2 = find(abs(mag(pidx:end) - pval) > thresh,1,'first');
                    if isempty(p_rel_idx2)
                        % same value keeps going to end of dataset
                        %
                        % Use midpoint between pidx and last point in
                        % dataset
                        plast = numel(mag);
                        pidx = floor((pidx + plast)/2);
                    else
                        % use midpoint as the "best position for the max"
                        %
                        % remember, p_rel_idx2 is PAST the contiguous
                        % stretch of values.  If it's 5, that means the
                        % last same-value entry is the 4th value in
                        % r2(pidx:end), which is pidx+3.  So it's off by
                        % TWO in a 1-based world.  Surprising, but that is
                        % what it is.
                        %
                        % If it's 1, it means r2 is the
                        % only value to consider.  Should never happen due
                        % to our entry test, but let's be safe.
                        if p_rel_idx2 > 1
                            if rem(p_rel_idx2,2) == 1
                                % odd # ... choice of 2 indices to use.
                                % Default is the first of the two:
                                pidx = pidx + floor((p_rel_idx2-2)/2);
                            else
                                pidx = pidx + (p_rel_idx2-2)/2;
                            end
                        end
                    end
                    
                    peakIdxVec(i) = pidx; % Record updated index
                end
            end
        end
        
        function peakIdxVec = getPeakMidpointsWithAngles(mag,ang,peakIdxVec)
            % Updates the index of the peak data points to the true
            % "midpoints" of each peak, if there are runs of
            % identical-valued maxima after the peak locations.
            %
            % pidxVec is a vector of peak locations found by findpeaks.
            %
            % In the presence of multiple identical contiguous peak values,
            % findpeaks returns the index of the first maxima.  However, we
            % sometimes require the "middle" value of such spans. For
            % example, the midpoint index best approximates the midpoint of
            % an "antenna lobe" pattern in polar data.
            %
            % The length of the identical-valued data may be the entire
            % dataset, and that could take time to evaluate.  So we do a
            % first test to see if the very NEXT index (in wrapped index
            % space) is identical in value.  We pursue a more complete test
            % only if this necessary precursor test passes.
            %
            % When useful, ang is used to choose the "best possible" peak
            % position that aligns with a 90-degree cardinal point.  Only
            % done when the choice between two adjacent points is
            % equivalent and one has a "nicer angle."
            
            thresh = 10*eps; % value differences <= thresh are "identical"
            
            Nr = numel(mag); 
            Ni = numel(peakIdxVec);  % a vector of peak indices
            for i = 1:Ni
                pidx = peakIdxVec(i);
                pval = mag(pidx);
                
                pnextIdx = 1 + rem(pidx,Nr); % +1 index offset implicit
                if abs(mag(pnextIdx) - pval) <= thresh
                    % Test data for contiguous span of max value.
                    % rotated index 1 is a global min, so we don't bother
                    % wrapping the test.  Just test from pidx to the end.
                    % That's overkill already.
                    %
                    % Since we now know there are at least TWO identical max
                    % values, we search from the first (pidx) to the end of the
                    % rotated data, and identify the first index where this is
                    % NOT true.
                    p_rel_idx2 = find(abs(mag(pidx:end) - pval) > thresh,1,'first');
                    if isempty(p_rel_idx2)
                        % same value keeps going to end of dataset
                        %
                        % Use midpoint between pidx and last point in
                        % dataset
                        plast = numel(mag);
                        pidx = floor((pidx + plast)/2);
                        
                        % Test if pidx is a multiple of 90-degrees. If so,
                        % use it - it's a "better" answer. Use pidx+1 only
                        % if pidx fails and pidx+1 is a 45-degree multiple.
                        if rem(ang(pidx),45)~=0 && rem(ang(pidx+1),45)==0
                            % this is a better answer
                            pidx = pidx+1;
                        end
                    else
                        % use midpoint as the "best position for the max"
                        %
                        % remember, p_rel_idx2 is PAST the contiguous
                        % stretch of values.  If it's 5, that means the
                        % last same-value entry is the 4th value in
                        % r2(pidx:end), which is pidx+3.  So it's off by
                        % TWO in a 1-based world.  Surprising, but that is
                        % what it is.
                        %
                        % If it's 1, it means r2 is the
                        % only value to consider.  Should never happen due
                        % to our entry test, but let's be safe.
                        if p_rel_idx2 > 1
                            if rem(p_rel_idx2,2) == 1
                                % odd # ... choice of 2 indices to use.
                                % Default is the first of the two:
                                pidx = pidx + floor((p_rel_idx2-2)/2);
                                
                                % Test if pidx is a multiple of 90-degrees.
                                % If so, use it - it's a "better" answer.
                                % Use pidx+1 only if pidx fails and pidx+1
                                % is a 45-degree multiple.
                                if rem(ang(pidx),45)~=0 && rem(ang(pidx+1),45)==0
                                    % this is a better answer
                                    pidx = pidx+1;
                                end
                            else
                                pidx = pidx + (p_rel_idx2-2)/2;
                            end
                        end
                    end
                    
                    peakIdxVec(i) = pidx; % Record updated index
                end
            end
        end
        
        function peaksIdx = findPolarPeaks(varargin)
            %findPolarPeaks Indices of peaks in polar magnitude data.
            %  findPolarPeaks(MAG) returns indices of peaks identified in
            %  polar magnitude data vector MAG, assuming MAG is uniformly
            %  distributed in angle over a complete circle, i.e., from
            %  (0:N-1)/N*360 degrees, where N = numel(MAG).
            %
            %  findPolarPeaks(MAG,CIRC) is used to specify logical flag
            %  CIRC as FALSE if MAG does not span a full circle. If
            %  omitted, CIRC is assumed to be TRUE.
            %
            %  findPolarPeaks(MAG,CIRC,N) specifies the maximum number N of
            %  peak indices to return. If omitted or set to inf, all peaks
            %  are returned.
            %
            %  findPolarPeaks(MAG,CIRC,N,LIM) specifies magnitude limits
            %  LIM as a 2-element vector [minMag maxMag], and removes
            %  peaks lower than minMag or higher than maxMag from the
            %  returned index vector.  If omitted, LIM = [min(M) max(M)].
            %
            %  findPolarPeaks(MAG,CIRC,N,LIM,OPTS) specifies a cell-vector
            %  of name/parameter properties that are passed directly to the
            %  FINDPEAKS function, used internally to find polar peaks. Two
            %  name/value pairs are automatically set and ignored if passed
            %  in OPTS: 'SortStr' (which is always set to 'descend') and
            %  'NPeaks' (which is set by the NPeaks input argument).
            %
            %  Requires Signal Processing Toolbox
            
            args = local_parsePolarPeaks(varargin{:});
            
            % findpeaks operates on a linear (i.e., non-circular,
            % non-wrapped) vector of data.  Our data wraps around a circle
            % (when CIRC=true) and is therefore contiguous at its start
            % and end points.  If we want to detect peaks that occur at the
            % first index, we must precondition the data.
            %
            % We "rotate" the mag data so the first datapoint is the global
            % min.  That way all peaks can be found.
            %
            % findpeaks is shape-preserving for its output vector: returns
            % a column vector if the input is a column vector.
            %
            if args.circ % mag is contiguous around a circle
                r2 = args.mag([args.minidx:end 1:args.minidx-1]);
            else
                r2 = args.mag;
            end
            [~,pidx] = internal.antfindpeaks(r2,args.findpeaks{:});
            
            %a  = args.angd;
            %a2 = a([idx:end 1:idx-1]);
            %pidx = internal.polariCommon.getPeakMidpointsWithAngles(r2,a2,pidx);
            
            pidx = internal.polariCommon.getPeakMidpoints(r2,pidx);
            
            if args.circ
                % Restore "wrapped" index values to reflect the original,
                % unwrapped data indices.  The -2 is correct: idx-1 is the
                % offset to add, with another -1 for the modulo
                % computation.
                peaksIdx = 1 + mod(pidx+args.minidx-2,numel(args.mag));
            else
                peaksIdx = pidx;
            end
        end
    end
    
    % Graphical display of angles/sectors
    %
    methods (Static)
        function [verts,faces] = sectorsPatchRounded(r,th,z0)
            % Compute polygons for sectors.
            %
            % dth is a vector of angle differences, computed using
            % angleDiff() or cangleDiff() by the caller, with a CCW
            % measurement, and normalized such that +ve increasing angle is
            % CCW in normalized coords.
            
            % Always traverses CCW using delta-angle increments.
            %
            % Uses round-sided pie-shaped wedges.
            %
            % faces: one per magnitude (Nr faces)
            %
            % vertices: 3 per face, but ...
            %    - origin point is shared across all faces
            %    - 2*Nr (mag,theta) points arc (span)
            %      Use Np points to smoothly draw an arc for each span
            %      Np = 2 minimum
            %         = one per degree for smooth
            %     -> 1 + Nr*Np total vertices
            
            % Determine total # vertices
            %
            % Use one theta-point every degree, minimum of 2 points.
            % Ex: th = [20.5 25.8]
            %     6 points equispaced over 5.3 degrees
            %     20.5 + 5.3/5*(0:5)
            %
            % All span polygons must include origin, in verts(1,1)
            
            % Compute the angle diff.
            %
            % Make -10 be +350, giving us the "always CCW" perspective
            % +ve increasing angle is CCW in normalized coords
            %
            % NOTE: Incoming angles could be a multiple of 360 degrees
            % apart already, which implies a full circle patch. Detect that
            % case, making sure the diff isn't truly zero.
            quick_dth = abs(th(2)-th(1));
            % quick_dth > pi: check that diff isn't close to zero.
            % and use thresh to test "close enough" to 2*pi*N
            if (quick_dth > pi) && rem(quick_dth,2*pi) < 1000*eps
                dth = 2*pi;  % one full circle
            else
                dth = internal.polariCommon.angleDiff(th);
            end
            
            % degrees per point xxx put in CONSTANT property?
            angle_density = 1.0 * pi/180;
            
            % dth_pts: #pts to compute per span
            %   - doesn't include +1, which is for origin point and needed to
            %     compute total vertices per span polygon
            dth_pts = max(2, ceil(dth/angle_density));
            
            Nv = 1 + sum(dth_pts); % 1 is for origin point, counted once in vert list
            verts = zeros(Nv,3);   % x,y,z coords (normalized)
            
            % # faces = # magnitudes
            Nr = numel(r);
            Nf = Nr;
            
            % faces matrix is trickier
            % - Nr faces
            % - each span has a different # of vertices, so max() # cols
            % - highest vertex count of any face is 1 (origin) plus max:
            Nv_max = 1 + max(dth_pts);
            faces = NaN(Nf,Nv_max);
            faces(:,1) = 1; % first vertex of every face is origin index
            
            % Clamp lower limit of radial distance.
            %
            % Do NOT go to zero, as the filled-patch no longer understands
            % direction and will not fill properly.
            r(r<=0) = eps;
            
            i_v = 2; % next write-index into verts (j=1 is origin)
            for i = 1:Nr
                r_i = r(i);
                
                % Compute the Nj vertices unique to this span
                %  (i.e., don't recompute origin!)
                Nj = dth_pts(i);
                dth_i = dth(i);
                
                % +ve theta is CCW
                th_vec = th(i) + dth_i*(0:Nj-1)'/(Nj-1);
                
                x = r_i.*cos(th_vec); % use last cos_th, sin_th
                y = r_i.*sin(th_vec);
                z = z0*ones(size(x));
                
                verts(i_v:i_v+Nj-1,:) = [x y z];
                
                % faces is simpler, but remember we're padding for the max-sized face.
                % So use indexing to overwrite these face indices into NaN-padding.
                %
                % Also, faces includes the origin point, so Nj+1 vertices form face.
                % We've written "1" index into (i,1) already.
                % So only write the remaining Nj vertices for this face:
                faces(i,2:Nj+1) = i_v:i_v+Nj-1;
                
                % Offset for next set of vertices
                i_v = i_v+Nj;
            end
            
            % Must set the "z" position for the origin point.
            %
            % Origin x,y is 0,0, but z must be in the same plane as rest of
            % the span.  So we set the full x,y,z of the origin here:
            verts(1,:) = [0 0 z0];
        end
    end
    
    % Formatted string conversion
    %
    methods (Static)
        function s = sprintfMaxNumTotalDigits(x,Ntot,forceCell)
            %sprintfMaxNumTotalDigits Convert value to string with at most N digits.
            %  sprintfMaxNumTotalDigits(X,N) converts scalar X to a string
            %  containing at most N total digits.  More digits will be
            %  added in order to fully display integer part of X.
            %
            % No need for negative N (which would imply integers can be
            % expressed without a fraction), since here there is an
            % explicit contract to suppress the fraction if the value is an
            % integer.
            %
            %  For vector X, returns a cell-vector with one string in each
            %  element.
            %
            % If forceCell is true, always returns a cell-vector, even for
            % a scalar input.
            
            if isempty(x)
                s = '[]';
            elseif ~isscalar(x) || nargin>2 && forceCell
                % Vector conversion
                
                N = numel(x);
                s = cell(1,N);
                for i = 1:N
                    s{i} = internal.polariCommon.sprintfMaxNumTotalDigits(x(i),Ntot);
                end
            else
                % Scalar conversion
                
                % ignore -ve value: we always suppress trailing fractional
                % 0's
                Ntot = abs(Ntot);
                
                % Nints is the # of digits total, including decimals. Display
                % all integer digits, plus a maximum of Nfrac decimals.
                if x==0
                    Nints = 1; % show one char for zero
                else
                    Nints = ceil(log10(abs(x)));
                end
                if Nints >= Ntot
                    % No additional digits to spare on fractional part
                    Nfrac = 0;
                else
                    Nfrac = Ntot - Nints; % max #
                end
                th = 1000*eps;

                if Nfrac==0
                    s = sprintf('%d',round(x));
                    
                elseif (abs(x)-fix(abs(x))<th)
                    % If within "threshold" of an int, it's an int
                    s = sprintf('%g',round(x));
                    
                else
                    % General case: number with int and fraction
                    
                    % Scale: shift left Nfrac places
                    sc = 10.^Nfrac;
                    x = round(x*sc);
                    
                    % Test trailing integer-part 0's, at most Nfrac # could
                    % be 0 that we care about.  Subtract up to Nfrac number
                    % of fraction digits if up to Nfrac number of trailing
                    % digits are zero.
                    y = x;
                    NIntDigitsAreZero = 0;
                    for j = 1:Nfrac
                        p10 = 10.^j;
                        y = round(y./p10)*p10;
                        if abs(y-x)>th
                            break;
                        end
                        NIntDigitsAreZero = NIntDigitsAreZero+1;
                    end
                    Nfrac = Nfrac-NIntDigitsAreZero;
                    
                    % scale it back after rounding:
                    x = x./sc;
                    s = sprintf('%*.*f',Nints,Nfrac,x);
                end
                
                % Suppress negative zero ('-0')
                %
                % Can occur when a small negative value is expressed
                % using just 2 fractional digits.
                if strcmpi(s,'-0')
                    s = '0';
                end
            end
        end
        
        function s = sprintfNumTotalDigits(x,Ntot,forceCell)
            %sprintfNumTotalDigits Convert value to string with N digits.
            %  sprintfNumTotalDigits(X,N) converts scalar X to a string
            %  containing N total digits.  More digits will be added in
            %  order to fully display integer part of X.
            %
            %  For vector X, returns a cell-vector with one string in each
            %  element.
            %
            % If forceCell is true, always returns a cell-vector, even for
            % a scalar input.
            
            if isempty(x)
                s = '[]';
            elseif ~isscalar(x) || nargin>2 && forceCell
                % Vector conversion
                
                N = numel(x);
                s = cell(1,N);
                for i = 1:N
                    s{i} = internal.polariCommon.sprintfNumTotalDigits(x(i),Ntot);
                end
            else
                % Scalar conversion
                
                suppressFracIfInt = Ntot < 0;
                Ntot = abs(Ntot);
                
                % Nints is the # of digits total, including decimals. Display
                % all integer digits, plus a maximum of Nfrac decimals.
                if x==0
                    Nints = 1; % show one char for zero
                else
                    Nints = ceil(log10(abs(x)));
                end
                if Nints >= Ntot
                    % No additional digits to spare on fractional part
                    Nfrac = 0;
                else
                    Nfrac = Ntot - Nints;
                end
                
                if Nfrac==0
                    s = sprintf('%d',round(x));
                    
                elseif suppressFracIfInt && (abs(x)-fix(abs(x))<1000*eps)
                    % If within 1000*eps of an int, it's an int
                    s = sprintf('%g',round(x));
                    
                else
                    %s = sprintf('%.*g',Nints+Nfrac,x);
                    sc = 10.^Nfrac;
                    x = round(x*sc)./sc;
                    s = sprintf('%*.*f',Nints,Nfrac,x);
                end
                
                % Suppress negative zero ('-0')
                %
                % Can occur when a small negative value is expressed
                % using just 2 fractional digits.
                if strcmpi(s,'-0')
                    s = '0';
                end
            end
        end
        
        function s = sprintfNumTotalDigitsAsVector(x,Ntot)
            %sprintfNumTotalDigitsAsVector Convert vector of values to
            %  a string representing a MATLAB vector, with N total
            %  digits for each value.
            %
            %  sprintfNumTotalDigitsAsVector (X,N) converts vector X to a
            %  string containing all N total digits.  Returns a string with
            %  MATLAB vector notation, suitable for evaluation using feval.
            %
            % Example:
            %    sprintfNumTotalDigitsAsVector([pi 1+pi 2+pi],3)
            %       returns '[3.14 4.14 5.14]'
            %
            %    sprintfNumTotalDigitsAsVector([pi 1+pi 2+pi],3)
            %       returns {'3.14' '4.14' '5.14'}
            
            % Get individual value conversions
            t = internal.polariCommon.sprintfNumTotalDigits(x,Ntot);
            if ~iscell(t)
                % was empty or a scalar
                 s = t;
            else
                % Assemble into a new string
                Nt = numel(t);
                s = '[';
                for i = 1:Nt
                    s = [s t{i}]; %#ok<AGROW>
                    if i < Nt
                        s = [s ', ']; %#ok<AGROW>
                    end
                end
                s = [s ']'];
            end
        end
        
        function s = sprintfMaxNumFracDigits(x,Nfrac,forceCell)
            %sprintfMaxNumFracDigits Convert value to string with at most N frac digits.
            %  sprintfMaxNumFracDigits(X,N) converts scalar X to a string
            %  containing all integer digits and at most N fractional
            %  digits.
            %
            %  For vector X, returns a cell-vector with one string in each
            %  element.
            %
            %  Example:
            %    sprintfMaxNumFracDigits([pi 1+pi 2+pi],2)
            %       returns {'3.14' '4.14' '5.14'}
            %    sprintfMaxNumFracDigitsAsVector([pi 1+pi 2+pi],2)
            %       returns '[3.14 4.14 5.14]'
            %
            %  If forceCell is true, always returns a cell-vector, even for
            %  a scalar input.
            
            if ~isscalar(x) || nargin>2 && forceCell
                % Vector conversion
                
                N = numel(x);
                s = cell(1,N);
                for i = 1:N
                    s{i} = internal.polariCommon.sprintfMaxNumFracDigits(x(i),Nfrac);
                end
            else
                % Scalar conversion
                
                % Ignore negative Nfrac, always prune trailing 0's
                Nfrac = abs(Nfrac);
                
                % Nints is the # of digits total, including decimals. Display
                % all integer digits, plus a maximum of Nfrac decimals.
                if x==0
                    Nints = 1; % show one char for zero
                else
                    % Nints can be negative, for x<1.
                    % Leave one int digit for zero.
                    Nints = max(1,ceil(log10(abs(x))));
                end
                %s = sprintf('%.*g',Nints+Nfrac,x);
                
                th = 1000*eps;
                if Nfrac==0
                    s = sprintf('%d',round(x));
                    
                elseif (abs(x)-fix(abs(x)) < th)
                    % If within "threshold" of an int, it's an int:
                    %
                    % No frac part
                    s = sprintf('%*g',Nints,round(x));
                else
                    % General case:
                    
                    sc = 10.^Nfrac;
                    x = round(x*sc);
                    
                    % Test trailing integer-part 0's, at most Nfrac # could
                    % be 0 that we care about.  Subtract up to Nfrac number
                    % of fraction digits if up to Nfrac number of trailing
                    % digits are zero.
                    y = x;
                    NIntDigitsAreZero = 0;
                    for j = 1:Nfrac
                        p10 = 10.^j;
                        y = round(y./p10)*p10;
                        if abs(y-x)>th
                            break;
                        end
                        NIntDigitsAreZero = NIntDigitsAreZero+1;
                    end
                    Nfrac = Nfrac-NIntDigitsAreZero;
                    
                    % scale it back after rounding:
                    x = x./sc;
                    
                    s = sprintf('%*.*f',Nints,Nfrac,x);
                end
                
                % Suppress negative zero ('-0')
                %
                % Can occur when a small negative value is expressed
                % using just 2 fractional digits.
                if strcmpi(s,'-0')
                    s = '0';
                end
            end
        end
        
        function s = sprintfNumFracDigits(x,Nfrac,forceCell)
            %sprintfNumFracDigits Convert value to string with N fractional digits.
            %  sprintfNumFracDigits(X,N) converts scalar X to a string
            %  containing all integer digits and N fractional digits.
            %
            %  If N is negative, suppresses fraction when value is an
            %  integer.
            %
            %  For vector X, returns a cell-vector with one string in each
            %  element.
            %
            %  Example:
            %    sprintfNumFracDigits([pi 1+pi 2+pi],2)
            %       returns {'3.14' '4.14' '5.14'}
            %    sprintfNumFracDigitsAsVector([pi 1+pi 2+pi],2)
            %       returns '[3.14 4.14 5.14]'
            %
            %  If forceCell is true, always returns a cell-vector, even for
            %  a scalar input.
            
            if ~isscalar(x) || nargin>2 && forceCell
                % Vector conversion
                
                N = numel(x);
                s = cell(1,N);
                for i = 1:N
                    s{i} = internal.polariCommon.sprintfNumFracDigits(x(i),Nfrac);
                end
            else
                % Scalar conversion
                
                % Detect negative Nfrac:
                suppressFracIfInt = Nfrac<0;
                Nfrac = abs(Nfrac);
                
                % Nints is the # of digits total, including decimals. Display
                % all integer digits, plus a maximum of Nfrac decimals.
                if x==0
                    Nints = 1; % show one char for zero
                else
                    % Nints can be negative, for x<1.
                    % Leave one int digit for zero.
                    Nints = max(1,ceil(log10(abs(x))));
                end
                %s = sprintf('%.*g',Nints+Nfrac,x);
                
                if suppressFracIfInt && (abs(x)-fix(abs(x)) < 1000*eps)
                    % If within 1000*eps of an int, it's an int
                    %
                    % No frac part
                    s = sprintf('%*g',Nints,round(x));
                else
                    s = sprintf('%*.*f',Nints,Nfrac,x);
                end
                
                % Suppress negative zero ('-0')
                %
                % Can occur when a small negative value is expressed
                % using just 2 fractional digits.
                if strcmpi(s,'-0')
                    s = '0';
                end
            end
        end
        
        function s = sprintfNumFracDigitsAsVector(x,Nfrac)
            %sprintfNumFracDigitsAsVector Convert vector of values to
            %  a string representing a MATLAB vector, with N fractional
            %  digits for each value.
            %
            %  If N is negative, suppresses fraction when value is an
            %  integer.
            %
            %  sprintfNumFracDigitsAsVector (X,N) converts vector X to a
            %  string containing all integer digits and N fractional
            %  digits.  Returns a string with MATLAB vector notation,
            %  suitable for evaluation using feval.
            %
            % Example:
            %    sprintfNumFracDigitsAsVector([pi 1+pi 2+pi],2)
            %       returns '[3.14 4.14 5.14]'
            %
            %    sprintfNumFracDigits([pi 1+pi 2+pi],2)
            %       returns {'3.14' '4.14' '5.14'}
            
            % Get individual value conversions
            t = internal.polariCommon.sprintfNumFracDigits(x,Nfrac);
            
            % Assemble into a new string
            Nt = numel(t);
            s = '[';
            for i = 1:Nt
                s = [s t{i}]; %#ok<AGROW>
                if i < Nt
                    s = [s ', ']; %#ok<AGROW>
                end
            end
            s = [s ']'];
        end
        
        function y = roundDigitsNoRescale(x,N)
            % Reduce precision of X such that only 3 digits are returned if
            % possible.  No change in units/exponent takes place.
            %
            % Examples:
            %   x = [234.56 -4.5678 0.00234 34567.8];
            %     234.56  -> 235
            %     -4.5678 -> -4.57
            %     0.00234 -> 0.00
            %     34567.8 -> 34568  (no solution)
            
            y = x;
            
            % To get 3 digits, we add 1 until we get a number > 2.  It
            % can't be =2, it must be > 2. Each add is equivalent to x10 of
            % the original number.  Then ceil() tells us how many integer
            % (non-fractional) digits are present, which will be 3.
            sgn = sign(x);
            lx = log10(abs(x));
            
            % Deal with small values: less than 3 mantissa digits
            Nl = N-1;
            delta = lx-Nl; % log values need to be > 2 for 3 digits
            selSm = delta <= 0;
            ySm = delta(selSm);
            fSm = ceil(-ySm); % exponent to scale
            % Impose max scale of 10^2 -> 0.00234 -> 0.00, not 0.002 etc
            fSm = min(fSm,Nl);
            y(selSm) = 10.^fSm .* x(selSm);
            
            % Deal with large values: more than 3 mantissa digits
            %{
            delta = lx-3;
            selLg = delta >= 1;
            yLg = delta(selLg);
            fLg = floor(-yLg);
            % Impose min scale of 10^-2 -> 34567 -> 34567, not 34500 etc
            fLg = max(fLg,-2);
            y(selLg) = 10.^fLg .* x(selLg);
            %}
            
            y = round(y);
            y(selSm) = y(selSm).*10.^(-fSm);
            %y(selLg) = y(selLg).*10.^(-fLg);
            y = y.*sgn;
        end
        
        function h = figRectNorm(fig,pos)
            % Draw a rect on a fig with its own full-figure norm axes.
            % figRectNorm(pos)
            % figRectNorm(fig,pos)
            %
            % Example:
            % internal.polariCommon.figRectNorm(getExtent(polari.getCurrentPlot))
            
            if nargin==1
                pos = fig;
                fig = gcf;
            end
            ax = axes('parent',fig, ...
                'Units','norm', ...
                'XLim',[0 1], ...
                'YLim',[0 1], ...
                'Position',[0 0 1 1], ...
                'Visible','off');
            h = line('parent',ax, ...
                'XData', [pos(1) pos(1) pos(1)+pos(3) pos(1)+pos(3) pos(1)], ...
                'YData', [pos(2) pos(2)+pos(4) pos(2)+pos(4) pos(2) pos(2)]);
        end
    end
    
    % String handlers
    %
    methods (Static)
        function arg = xlatExtendedASCII(arg,dir)
            %xlatExtendedASCII Translate extended ASCII codes within a string.
            %   xlatExtendedASCII(STR) translates character-code keywords
            %   identified within a string or a cell-array of strings to
            %   single-character symbols.
            %
            %  Symbols:
            %     '#deg', '#infin' (degree sign, infinity)
            %     '#plusmn'        (plus-or-minus)
            %     '#copy', '#reg'  (copyright, registered signs)
            %     '#sup1', '#sup2', '#sup3' (numeric superscripts)
            %     '#micro', '#dagger'
            %     '#nabla'
            %     '#ohm'
            %
            %  Greek alphabet:
            %     '#alpha'      '#Alpha'
            %     '#beta'       '#Beta'
            %     '#gamma'      '#Gamma'
            %     '#delta'      '#Delta'
            %     '#epsilon'    '#Epsilon'
            %     '#zeta'       '#Zeta'
            %     '#eta'        '#Eta'
            %     '#theta'      '#Theta'
            %     '#iota'       '#Iota'
            %     '#kappa'      '#Kappa'
            %     '#lambda'     '#Lambda'
            %     '#mu'         '#Mu'
            %     '#nu'         '#Nu'
            %     '#xi'         '#Xi'
            %     '#omicron'    '#Omicron'
            %     '#pi'         '#Pi'
            %     '#rho'        '#Rho'
            %     '#sigma'      '#Sigma'       '#sigmaf'
            %     '#tau'        '#Tau'
            %     '#upsilon'    '#Upsilon'
            %     '#phi'        '#Phi'
            %     '#chi'        '#Chi'
            %     '#psi'        '#Psi'
            %     '#omega'      '#Omega'
            %
            % Character codes are HTML names, without the ampersand prefix
            % and semicolon postfix characters.  Ex: '&deg;' is the HTML
            % name for degree sign, while '#deg' is the code here.
            %
            % xlatExtendedASCII(STR,DIR) specifies translation direction
            % DIR to be 'forward' or 'reverse'.  If omitted, DIR is
            % 'forward'.
            
            patterns = { ...
                '#copy',  char(169);
                '#reg',   char(174);
                '#dagger',char(134);
                '#copy',  char(169);
                '#deg',   char(176);
                '#plusmn',char(177);
                '#infin', char(8734);
                '#micro', char(181);
                '#sup1',  char(185);
                '#sup2',  char(178);
                '#sup3',  char(179);
                '#nabla', char(8711);
                '#ohm',   char(937);
                '#Alpha', char(913);
                '#Beta',  char(914);
                '#Gamma', char(915);
                '#Delta', char(916);
                '#Epsilon',char(917);
                '#Zeta',  char(918);
                '#Eta',   char(919);
                '#Theta', char(920);
                '#Iota',  char(921);
                '#Kappa', char(922);
                '#Lambda',char(923);
                '#Mu',    char(924);
                '#Nu',    char(925);
                '#Xi',    char(926);
                '#Omicron',char(927);
                '#Pi',    char(928);
                '#Rho',   char(929);
                '#Sigma', char(931);
                '#Tau',   char(932);
                '#Upsilon',char(933);
                '#Phi',   char(934);
                '#Chi',   char(935);
                '#Psi',   char(936);
                '#Omega', char(937);
                '#alpha', char(945);
                '#beta',  char(946);
                '#gamma', char(947);
                '#delta', char(948);
                '#epsilon',char(949);
                '#zeta',  char(950);
                '#eta',   char(951);
                '#theta', char(952);
                '#iota',  char(953);
                '#kappa', char(954);
                '#lambda',char(955);
                '#mu',    char(956);
                '#nu',    char(957);
                '#xi',    char(958);
                '#omicron',char(959);
                '#pi',    char(960);
                '#rho',   char(961);
                '#sigmaf', char(962);
                '#sigma', char(963);
                '#tau',   char(964);
                '#upsilon',char(965);
                '#phi',   char(966);
                '#chi',   char(967);
                '#psi',   char(968);
                '#omega', char(969)};
            
            if nargin == 0
                % Just return list of supported codes
                if nargout > 0
                    arg = patterns(:,1);
                else
                    disp('Text markup symbols supported:');
                    disp(patterns(:,1));
                end
                return % EARLY EXIT
            end
            
            if ~ischar(arg) && ~(isstring(arg) && isscalar(arg)) && ~iscellstr(arg) && ~isstring(arg)
                error('Input must be a string or a cell-array of strings');
            end
            if nargin<2 || strcmpi(dir,'forward')
                % forward conversion: keyword to char-code
                i_src = 1;
                i_dst = 2;
            else
                % reverse convesion: char-code to keyword
                i_src = 2;
                i_dst = 1;
            end
            
            % Translate between "#keyword" and char(code)
            if ischar(arg) || (isstring(arg) && isscalar(arg))
                % string to use as label for all datasets
                arg = xlatExtendedASCII_convertOneCharMatrix(arg,patterns,i_src,i_dst);
            else
                % cell string, one label per dataset
                for i = 1:numel(arg)
                    arg{i} = xlatExtendedASCII_convertOneCharMatrix(arg{i},patterns,i_src,i_dst);
                end
            end
        end
        
        function s = removeExtendedASCII(s)
            % Remove ASCII chars > 255.
            %
            % Accepts char vector and char matrix only (not cellstr).
            
            assert(ischar(s));
            t = cellstr(s);
            for i = 1:numel(t)
                t_i = t{i};
                t_i(fix(t_i) > 255) = '';
                t{i} = t_i;
            end
            s = char(t);
        end
        
        function strs = convertEmbeddedCRsToCharMatrices(strs)
            % input: string vector or cell-string, possibly containing CR
            % chars.
            %
            % returns: string or cell-string, containing char vectors or
            %   char matrices.  No CR chars remain.  Any string with a CR
            %   is converted to a multi-row string matrix.
            
            is_cell = iscellstr(strs);
            if ~is_cell
                strs = {strs};
            end
            CR = sprintf('\n');
            for i = 1:numel(strs)
                % Convert each entry in cell-string strs to a cellstring.
                % Each entry could be:
                %   - vector string -> now a scalar cell in c_i
                %   - char matrix -> now a cell-vector in c_i
                c_i = cellstr(strs{i});
                
                Nc = numel(c_i);
                for j = 1:Nc % Nc>1 only if top-level entry was a char matrix
                    c_j = c_i{j}; % vector or string matrix
                    idx = find(c_j == CR);
                    Ncr = numel(idx);
                    if Ncr > 0
                        t = cell(Ncr+1,1);
                        i2 = 0;
                        for k = 1:Ncr
                            i1 = i2;
                            i2 = idx(k);
                            t{k} = c_j(i1+1 : i2-1);
                        end
                        t{end} = c_j(i2+1:end);
                        c_i{j} = char(t); % convert to string matrix
                    end
                end
                strs{i} = char(c_i);
            end
            if ~is_cell
                strs = strs{1};
            end
        end
        
        function strs = convertEmbeddedCellsToCharMatrices(strs)
            % Remove cell-strings within the top-level cell-string.
            %
            % Interactive use of Legend enables these to be created.
            % iscellstr() fails on those, and they aren't supported here.
            
            % iscellstr returns false if contents of top-level cell vector contains
            % anything other than a string.  It supports vector and matrix chars.
            %
            % But it returns false if cellstrs are put inside the cell-string.
            % We don't want to support that, either.
            %
            % So this offers quick top-level detection of this issue:
            if ~iscellstr(strs) && ~isstring(strs) && ~ischar(strs) && ~(isstring(strs) && isscalar(strs))
                for i = 1:numel(strs)
                    if iscell(strs{i})
                        % i'th label has an embedded cell.
                        % Convert it to a char matrix:
                        strs{i} = char(strs{i});
                    end
                end
            end
        end
        
        function strs = convertStringMatrixToCR(strs)
            % Convert string matrices to vector strings with embedded CR chars.
            % Accept cell-vector at top level.
            
            CR = sprintf('\n');
            if iscell(strs)
                Ns = numel(strs);
                for i = 1:Ns
                    s_i = strs{i};
                    Nr = size(s_i,1);
                    if ischar(s_i) || (isstring(s_i) && isscalar(s_i)) && Nr > 1 % string matrix
                        t = [s_i'; repmat(CR,1,Nr)];
                        strs{i} = t(:)';
                    end
                end
            else
                Nr = size(strs,1);
                if Nr > 1
                    t = [strs'; repmat(CR,1,Nr)];
                    strs = t(:)';
                    strs = strs(1:end-1); % don't add trailing CR
                end
            end
        end
    end
    
    % String support
    %
    methods (Static)
        function c = getUTFCircleChar(ch)
            %getUTFCircleChar Return UTF-8 character for circled character.
            %   getUTFCircleNumber(ch) returns circled characters, encoded as UTF-8
            %   unicode characters.  Accepts A-Z, a-z, 0-9.
            
            % UniCode: UTF-8
            %   http://www.utf8-chartable.de/unicode-utf8-table.pl
            %
            %        -> circled A   -> 226 146 182
            %        -> circled K   -> 226 147 128
            %        -> circled a   -> 226 147 144
            %        -> circled 1-9 -> 226 145 160
            
            if isstrprop(ch,'upper')
                if fix(ch) >= fix('K')
                    dec = [226 147 128]; % 'K'
                    ch = ch-fix('K');
                else % 'A'
                    dec = [226 146 182]; % 'A'
                    ch = ch-fix('A');
                end
                
            elseif isstrprop(ch,'lower')
                dec = [226 147 144]; % 'a'
                ch = ch-fix('a');
                
            elseif isstrprop(ch,'digit')
                if strcmp(ch,'0')
                    dec = [226 147 170]; % '0'
                    ch = 0;
                else
                    dec = [226 145 160]; % '1'
                    ch = ch-fix('1');
                end
            end
            c = native2unicode(dec+[0 0 fix(ch)],'UTF-8');
        end
        
        function c = getUTFCircleNumber(nums)
            %getUTFCircleNumber Return UTF-8 character for circled integer.
            %   getUTFCircleNumber(N) returns characters representing circled numbers,
            %   encoded as UTF-8 unicode characters.  Accepts a vector of integers in
            %   range 0 to 20.
            
            % UniCode: UTF-8
            %   http://www.utf8-chartable.de/unicode-utf8-table.pl
            %
            % Superscripts
            % U+     -> circled 0  -> dec 226 147 170
            % U+2460 -> circled 1  -> dec 226 145 160
            % U+2461 -> circled 2  -> dec 226 145 161
            %  ...
            % U+2473 -> circled 20  -> dec 226 145 179
            
            validateattributes(nums,{'double','single'}, ...
                {'real','finite','integer','>=',-20,'<=',20});
            
            for i = numel(nums):-1:1
                n_i = nums(i);
                if n_i == 0
                    dec = [226 147 170];
                elseif n_i > 0 % positive
                    dec = [226 145 159+n_i]; % 1-20
                    
                    % DINGBAT CIRCLED SANS-SERIF DIGIT ONE
                    %dec = [226 158 127+n_i]; % dingbat, 1-10
                else % negative
                    
                    % DINGBAT NEGATIVE CIRCLED DIGIT ONE
                    %dec = [226 157 181-n_i];
                    
                    % DINGBAT NEGATIVE CIRCLED SANS-SERIF DIGIT ONE
                    %dec = [226 158 137-n_i]; % dingbat, 1-10
                end
                c(i) = native2unicode(dec,'UTF-8');
            end
        end
        
        function c = getUTFSubscriptNumber(num,style)
            % Return UTF-8 encoded string with subscript integers.
            % Pass scalar integer, can be multiple digits.
            %
            % getUTFSubscriptNumber(N) returns subscripted string containing digits of
            % N encoded as UTF-8 unicode characters.
            %
            % getUTFSubscriptNumber(N,STYLE) where STYLE is either 'superscript' or
            % 'subscript', returns string using either superscript or subscript
            % characters.
            
            % UniCode: UTF-8
            %   http://www.utf8-chartable.de/unicode-utf8-table.pl
            %
            % Superscripts
            % U+2070 -> superscript 0  -> dec 226 129 176
            % U+00B9 -> superscript 1  -> dec 194 185
            % U+00B2 -> superscript 2  -> dec 194 178
            % U+00B3 -> superscript 3  -> dec 194 179
            % U+2074 -> superscript 4  -> dec 226 129 180
            %  ...
            % U+2079 -> superscript 9  -> dec 226 129 185
            %
            % Subscripts
            % U+2080 -> subscript 0    -> dec 226 130 128
            % ...
            % U+2089 -> subscript 9    -> dec 226 130 137
            
            validateattributes(num,{'double','single'}, ...
                {'real','finite','nonnegative','integer'});
            if num==0
                digits = 0;
            else
                Nd = ceil(log10(num+1));
                digits = zeros(1,Nd);
                for i = 1:Nd
                    hi = floor(num/10);
                    lo = num-hi*10;
                    digits(Nd+1-i) = lo;
                    num = hi;
                end
            end
            
            % Convert each digit to unicode
            if nargin<2
                str = 'subscript';
            else
                str = validatestring(style,{'superscript','subscript'},mfilename,'STYLE',2);
            end
            if strcmpi(str,'subscript')
                % Subscript
                for i = numel(digits):-1:1
                    c(i) = native2unicode([226 130 128+digits(i)],'UTF-8');
                end
            else
                % Superscript
                for i = numel(digits):-1:1
                    d_i = digits(i);
                    if d_i == 0
                        code = [226 129 176];
                    elseif d_i == 1
                        code = [194 185];
                    elseif d_i == 2
                        code = [194 178];
                    elseif d_i == 3
                        code = [194 179];
                    else
                        code = [226 129 176+d_i];
                    end
                    c(i) = native2unicode(code,'UTF-8');
                end
            end
        end
    end
    
    % Input-argument parameter sweeps
    %
    methods (Static)
        function y = fevalArgSets(fcn,argSlices)
            %fevalArgSets Evaluate function at each set of argument slices.
            %   Y = fevalSlices(@FUN,SLICES) returns a cell-vector Y, where each
            %   element of Y is the result of evaluating function FUN at each set of
            %   arguments defined by the cell-vector SLICES.
            
            N = numel(argSlices);
            y = cell(1,N);
            for i = 1:N
                y{i} = feval(fcn,argSlices{i}{:});
            end
            
        end
        
        function allSlices = createArgSets(varargin)
            %createArgSets Argument slicer for vector argument sets.
            %   Return a list (cell-vector) of cells, where each cell
            %   contains a set of scalar values.  Each set is created from
            %   a combination of scalars derived from the input argument
            %   list.  Scalar input arguments are replicated such that they
            %   appear in each output set, while vector input arguments are
            %   sliced into scalars.  Inputs must therefore be either
            %   scalars or be vectors with identical lengths.
            %
            % Ex: argSlice(1:3)
            %      y = { {1},{2},{3} }
            % Ex: argSlice(2,1:3)
            %      y = { {2,1},{2,2},{2,3} }
            % Ex: argSlice(1:3,11:13)
            %      y = { {1,11},{2,12},{3,13} }
            % Ex: argSlice(1:3,{'x','y','z'},100)
            %      y = { {1,'x',100},{2,'y',100},{3,'z',100} }
            
            % Get # of entries in each input
            % strings: call a string "1" entry
            % numeric vectors: use actual numel
            Nargs = numel(varargin);  % each slice has this many args
            Nel = zeros(1,Nargs);
            for i = 1:Nargs
                v_i = varargin{i};
                if ischar(v_i) || (isstring(v_i) && isscalar(v_i))
                    Nel(i) = 1;
                else
                    Nel(i) = numel(v_i);
                end
            end
            %Nel = cellfun(@numel,varargin); % # elements in each arg, Nargs # of elements
            
            % Check that each arg has 1 or Nel elements
            Nel_nonscalar = Nel;
            Nel_nonscalar(Nel==1) = [];
            if numel(unique(Nel_nonscalar)) > 1
                error('All arguments must be scalar or have the same length.');
            end
            if isempty(Nel_nonscalar)
                Nmaxel = 1;
            else
                Nmaxel = Nel_nonscalar(1); % all are the same non-scalar
            end
            
            allSlices = cell(1,Nmaxel); % # slices to return
            oneSlice = cell(1,Nargs); % each slice
            for i = 1:Nmaxel
                for j = 1:Nargs
                    t = varargin{j};
                    if Nel(j)==1
                        val = t;
                    else
                        if iscell(t)
                            val = t{i};
                        else
                            val = t(i);
                        end
                    end
                    oneSlice{j} = val;
                end
                allSlices{i} = oneSlice;
            end
            
        end
    end
end

function arg = xlatExtendedASCII_convertOneCharMatrix(arg,patterns,i_src,i_dst)
% Translate one string.
%
% May be a char vector or a char matrix

assert(~iscell(arg))

if size(arg,1) < 2
    % Allow empty, or single-row string
    % -- careful not to allow an Nx1 string, which is multi-row
    %
    % char vector (simple string)
    for i = 1:size(patterns,1)
        arg = strrep(arg,patterns{i,i_src},patterns{i,i_dst});
    end
else
    % Char matrix
    t = cellstr(arg);
    N = numel(t);
    for i = 1:size(patterns,1)
        % Replace each row separately.
        % strrep fails on char matrices.
        for k = 1:N
            t{k} = strrep(t{k},patterns{i,i_src},patterns{i,i_dst});
        end
    end
    arg = char(t); % convert back to char matrix
end
end

function x = trigFixup(x)
% Derived from makehgtform local function "fixup"
vals = [0 1 -1];
for i = 1:numel(x)
    sel = abs(x(i)-vals) < 2*eps;
    if any(sel)
        x(i) = vals(sel);
    end
end

end

function args = local_parsePolarPeaks(mag,varargin)
% Parses input arg syntax:
%  (mag,circ,Npeaks,magLim,userOpts)
%
%
%  All args after mag are optional but must appear in order if specified.
%  mag - vector, may be empty but MUST be present
%  circ - scalar logical
%  Npeaks - scalar numeric, may be inf
%  magLim - 2-element vector
%  userOpts - userOpts is a cell-vector, may be empty
%
% Example inputs arg lists:
%   (mag)
%   (mag,circ)
%   (mag,Npeaks)
%   (mag,Npeaks,magLim,userOpts)
%   (mag,Npeaks,userOpts)
%   (mag,userOpts)
%
% Returns struct:
%   args.mag
%   args.minidx
%   args.circ
%   args.findpeaks_args
            
% 0: PARSE DATA
%
narginchk(1,5); % need at least mag input, others are optional
if ~isempty(mag)
    validateattributes(mag, {'numeric'}, ...
        {'vector','real'},'findPolarPeaks');
end
args.mag = mag;

% 1: PARSE OPTIONAL ARGS
%
% Define validateattributes for each optional arg
%   { validateattributes_args, field_name, default_value, canBeEmpty; ...}
argAttrs = { ...
    {{'logical'},{'scalar','real'},'findPolarPeaks'},'circ',true,false; ...
    {{'numeric'},{'scalar','real','nonnegative'},'findPolarPeaks'},'Npeaks',inf,false; ...
    {{'numeric'},{'vector','real','numel',2},'findPolarPeaks'},'magLim',[],false; ...
    {{'cell'},{'vector','even'},'findPolarPeaks'},'userOpts',{},true };

Nattrs = size(argAttrs,1); % # rows
Nargs = numel(varargin);
argInputIdx = 1;
argAttrIdx = 1; % next arg type to check

% Create DEFAULT values in fields of return struct
% Ex: args.mag = <defaultValue>
for i = 1:Nattrs
    args.(argAttrs{i,2}) = argAttrs{i,3};
end

while (argInputIdx <= Nargs) && (argAttrIdx <= Nattrs)
    try %#ok<TRYNC>
        val = varargin{argInputIdx};
        
        % Can arg be empty as an exception?
        % - If arg cannot be empty, just do validate
        % - If arg can be empty, and it's NOT empty, do validate
        % - Otherwise (can be and is empty), SKIP validate
        t = argAttrs{argAttrIdx,4};
        if ~t || ~isempty(val)
            % Process via validateattributes
            t = argAttrs{argAttrIdx,1}; % t is a cell-vector
            validateattributes(val,t{:}); % may throw an error
        end
        
        % Create field in args return struct
        % Ex: args.mag = <userValue>
        args.(argAttrs{argAttrIdx,2}) = val;
        
        % Successfully parsed arg - move to next input arg:
        argInputIdx = argInputIdx+1;
    end
    
    % Error or not, bump argAttr index
    argAttrIdx = argAttrIdx+1;
end

% If we parsed all user args, argInputIdx should now be Nmax+1
if argInputIdx <= Nargs
    error('Invalid input argument #%d.\n',argInputIdx-1);
end

% 2: Post-process: override defaults, add other fields
%

% magLim:
%
mag = args.mag;
if isempty(args.magLim)
    % Magnitude limits of data - can be overriden to match a
    % "displayed" range of magnitudes versus actual data
    % magnitudes
    %
    % We need index of true mag min as well.  Since the min/max
    % vector wasn't passed by the caller, we compute that here,
    % and find the true mag min index at the same time:
    [minmag,minidx] = min(mag);
    args.magLim = [minmag max(mag)];
else
    % min/max vector passed in by caller
    %
    % We can't rely on magLim containing the actual min, and it
    % doesn't include its index in any case, so we compute this
    % index now:
    [~,minidx] = min(mag); % true min, regardless of magLim value
end
args.minidx = minidx;

% NPeaks:
%
if isinf(args.Npeaks)
    args.Npeaks = []; % means "all" for findpeaks()
end

% UserOpts:
%
% - options that CAN be overridden by user opts:
%   Default minimum peak prominence is 5% of the specified dynamic range
magLim = args.magLim;
if isempty(magLim)
    mpp = 0;
    mph = 0;
else
    mpp = (magLim(2)-magLim(1))*0.05;
    mph = magLim(1);
end
t = {'MinPeakHeight',mph, 'MinPeakProminence',mpp};
if ~isempty(args.userOpts)
    t = [t args.userOpts(:).'];
end
% - options that CANNOT be overridden by user opts:
args.findpeaks = [t {'SortStr','descend', 'NPeaks',args.Npeaks}];

end
