function [cc, hh, bb] = altitudeEnvelopeContour(varargin)
%ALTITUDEENVELOPECONTOUR Draw altitude envelope contour plot
%
%   ALTITUDEENVELOPECONTOUR( LOADFACTOR ) draws a contour plot of the LOADFACTOR
%   matrix in the x-y plane. The x-coordinates of the vertices correspond
%   to the column indices of LOADFACTOR and the y-coordinates correspond to the row
%   indices of LOADFACTOR. The contour levels are chosen automatically by contour.
%
%   ALTITUDEENVELOPECONTOUR( AIRSPEED, ALTITUDE, LOADFACTOR ) draws a contour plot
%   of LOADFACTOR using vertices from the mesh defined by AIRSPEED and ALTITUDE.
%   AIRSPEED and ALTITUDE can be vectors or matrices.
%
%   ALTITUDEENVELOPECONTOUR(___, LEVELS) plots an altitude envelope contour
%   specified by the desired levels LEVELS.
%
%   ALTITUDEENVELOPECONTOUR(___, LINESPEC) plots an altitude envelope contour
%   specified by the desired line spec LINESPEC.
%
%   ALTITUDEENVELOPECONTOUR(___, Name,Value ) plots an altitude envelope contour
%   specified by one or more Name,Value pairs.
%
%   ALTITUDEENVELOPECONTOUR(AX, ___) draws an altitude contour plot onto the axes AX.
%
%   [C,H,B] = ALTITUDEENVELOPECONTOUR(...) returns contour matrix C, a contour object H,
%   and a vector of boundarylines objects B.
%   C and H can be used as inputs to CLABEL. The structure of a contour matrix
%   is described in the documentation for CONTOURC.
%
% Inputs:
%   AX         :Optional scalar handle to a valid axes. By default, ALTITUDEENVELOPECONTOUR
%               plots to the current axes. 
%               To get the current axes, use the gca function.
%
%   AIRSPEED   :Optional numeric vector or matrix specifying the airspeed for each
%               corresponding index of LOADFACTOR. By default, AIRSPEED is the
%               column indices of LOADFACTOR.
%
%   ALTITUDE   :Optional numeric vector or matrix specifying the altitude for each
%               corresponding index of LOADFACTOR. By default, ALTITUDE is the
%               row indices of LOADFACTOR.
%
%   LOADFACTOR :Numeric matrix specifying the load factor for each AIRSPEED
%               and ALTITUDE. LOADFACTOR is usually defined as Lift / Weight with units of g's.
%               However, LOADFACTOR can be any matrix of numeric data.
%
%   LEVELS     :Optional numeric scalar or vector specifying the levels to contour.
%               If LEVELS is scalar, LEVELS specifies the number of contour lines
%               and the contour levels are chosen automatically by contour.
%               If LEVELS is a vector, LEVELS specifies the number and levels of
%               contour lines to plot.
%               To draw the contours at one height (k), specify levels as a
%               two-element row vector [k k].
%               By default, LEVELS is chosen automatically by contour.
%
%   LINESPEC :Scalar string specifying the line specification for the
%             plotted boundaryline.
%             LINESPEC is made of one element from any or all the following 3 columns. 
%          Line Color:         Marker Style:           Line Style:
%          b     blue          .     point              -     solid
%          g     green         o     circle             :     dotted
%          r     red           x     x-mark             -.    dashdot 
%          c     cyan          +     plus               --    dashed   
%          m     magenta       *     star             (none)  no line
%          y     yellow        s     square
%          k     black         d     diamond
%          w     white         v     triangle (down)
%                              ^     triangle (up)
%                              <     triangle (left)
%                              >     triangle (right)
%                              p     pentagram
%                              h     hexagram
%
% Name,Values:
%   MinimumAltitude  :Numeric scalar or N-by-2 matrix specifying the minimum
%                     altitude boundary.
%                     If MinumumAltitude is scalar, MinumumAltitude is assumed to
%                     be a horizontal line that intersects with the X limits
%                     of the axes or with the intersection of MinimumSpeed and
%                     MaximumSpeed, if specified.
%                     If MinimumAltitude is an N-by-2 matrix, each row is an
%                     AIRSPEED and ALTITUDE point of the boundary.
%                     LOADFACTOR values below MinimumAltitude are marked as NaN.
%
%   MaximumAltitude  :Numeric scalar or N-by-2 matrix specifying the maximum
%                     altitude boundary.
%                     If MaximumAltitude is scalar, MaximumAltitude is assumed to
%                     be a horizontal line that intersects with the X limits
%                     of the axes or with the intersection of MinimumSpeed and
%                     MaximumSpeed, if specified.
%                     If MaximumAltitude is an N-by-2 matrix, each row is an
%                     AIRSPEED and ALTITUDE point of the boundary.
%                     LOADFACTOR values above MaximumAltitude are marked as NaN.
%
%   MinimumSpeed     :Numeric scalar or N-by-2 matrix specifying the minimum
%                     speed boundary.
%                     If MinimumSpeed is scalar, MinimumSpeed is assumed to
%                     be a vertical line that intersects with the Y limits
%                     of the axes or with the intersection of MinimumAltitude and
%                     MaximumAltitude, if specified.
%                     If MinimumSpeed is an N-by-2 matrix, each row is an
%                     AIRSPEED and ALTITUDE point of the boundary.
%                     LOADFACTOR values behind MinimumSpeed are marked as NaN.
%
%   MaximumSpeed     :Numeric scalar or N-by-2 matrix specifying the maximum
%                     speed boundary.
%                     If MaximumSpeed is scalar, MaximumSpeed is assumed to
%                     be a vertical line that intersects with the Y limits
%                     of the axes or with the intersection of MinimumAltitude and
%                     MaximumAltitude, if specified.
%                     If MaximumSpeed is an N-by-2 matrix, each row is an
%                     AIRSPEED and ALTITUDE point of the boundary.
%                     LOADFACTOR values in front of MaximumSpeed are marked as NaN.
%
%   BoundaryXData    :Numeric vector specifying the X data for the boundary line.
%   BoundaryYData    :Numeric vector specifying the Y data for the boundary line.
%
%   ClipContour   :Scalar on-off switch that controls the contour data display behavior
%                     Must be one of the following:
%                       'on'  :Remove the contour data outside the specified boundary lines.
%                       'off' :Show all contour data and specified boundary lines.
%                     By default, ClipContour is 'on'.
%
%   ResolveBoundary  :Scalar on-off switch that controls the boundary line
%                     intersection and enclosing behavior.
%                     Must be one of the following:
%                       'on'  :Resolve the boundary line segments to form a
%                              closed boundary around the line intersection points.
%                              Boundary data points outside the resolved boundary
%                              are removed. 
%                              This method produces a well-formed boundary,
%                              but does not allow infinite limits or unconnected
%                              boundary lines.
%                       'off' :Do not resolve the boundary line segments. This method
%                              leaves the boundary line data unmodified to
%                              allow infinite limits and unconnected boundary lines.
%                              However, this method might produce a malformed
%                              boundary that affects the ClipContour behavior.
%                     By default, ResolveBoundary is 'on'.
%
% The MinimumAltitude, MaximumAltitude, MininumSpeed, and MaximumSpeed Name,Values
% work together to form the sides of a complete boundary line around the
% LOADFACTOR data. For more control of the boundary line, use the BoundaryXData
% and BoundaryYData Name,Values.
%
% LOADFACTOR values outside of the connected polygon formed by BoundaryXData and
% BoundaryYData are marked as NaN. Specifying BoundaryXData and BoundaryYData
% overwrites any of the specified MinimumAltitude, MaximumAltitude,
% MininumSpeed, or MaximumSpeed Name,Values.
%
% Additionally, all Name,Values of contour can be specified.
% For a list of these Name,Values, see the documentation for contour.
%
% Outputs:
%   C :Numeric contour matrix.
%   H :Scalar contour graphics object.
%   B :Vector of one or more boundaryline objects.
%
% Examples:
%   Plot an altitude envelope contour with 20 levels and no boundaries.
%     [speed,alt,loadfactor] = peaks();
%     altitudeEnvelopeContour(speed,alt,loadfactor,20)
%
%   Plot an altitude envelope contour specifying constant boundary lines
%   and turn off boundary line intersection clipping.
%     [speed,alt,loadfactor] = peaks();
%     altitudeEnvelopeContour(speed,alt,loadfactor,"MinimumSpeed",0,"MaximumAltitude",2,"ResolveBoundary","off")
%
%   Plot an altitude envelope contour specifying custom boundary data.
%     x = linspace(-15,15);
%     y = linspace(-15,15);
%     [X,Y] = meshgrid(x,y);
%     Z = sin(X)+cos(Y);
%     t = linspace(2*pi,0);
%     boundaryX = 16*sin(t).^3;
%     boundaryY = 13*cos(t) - 5*cos(2*t) - 2*cos(3*t) - cos(4*t);
%     altitudeEnvelopeContour(X,Y,Z,"BoundaryXData",boundaryX,"BoundaryYData",boundaryY)

% See also boundaryline, contour, contourc, clabel

% Copyright 2021 The MathWorks, Inc.

narginchk(1,inf)

args = matlab.graphics.internal.convertStringToCharArgs(varargin);

% get the parent nvp g2447291
[args, cax] = Aero.internal.namevalues.getAxesFromFirstArgumentOrNameValuePair(args);

try
    [contourargs, boundaryargs] = parseAltitudeEnvelopeInputs(args);
catch ME
    throw(ME)
end

[c,h] = contour(cax, contourargs{:});

if isempty(c)
    % if c is empty, there was nothing to contour. So produce an empty boundary
    b = repmat(Aero.graphics.primitive.BoundaryLine, 0);
else
    % we passed contour, so now we can fill in any leftover missing values
    % to complete the box
    xbounds = [min(h.XData,[],'all'), max(h.XData,[],'all')];
    ybounds = [min(h.YData,[],'all'), max(h.YData,[],'all')];

    if isempty(boundaryargs.Groups)
        [boundaryData, groups] = Aero.internal.math.clipBoundingBoxIntersections(boundaryargs.Boundary, xbounds, ybounds);
    else
        boundaryData = Aero.internal.math.fillBoundingBoxLimits(boundaryargs.Boundary, inf, xbounds, ybounds);
        groups = boundaryargs.Groups;
    end

    % plot a boundary line for each group
    b = repmat(Aero.graphics.primitive.BoundaryLine, max(groups),1);
    for i = 1:numel(b)
        idx = groups == i;
        b(i) = boundaryline(cax, boundaryData(idx,1), boundaryData(idx,2));

        % intersect seems to like flipping the points around.
        % so flip it back
        b(i).FlipBoundary = boundaryargs.ResolveBoundary;
    end
end

if boundaryargs.ClipContour && ~isempty(b)
    % to prevent warnings being thrown in "inpolygon", we have to replace
    % the inf's used with "inf", aka, a really big number below 1e150 as that
    % is the absolute limit allowed in inpolygon.
    
    % 1e50 seems to be fairly stable about not throwing the warning.
    idx = ~isfinite(boundaryData);
    boundaryData(idx) = 1e50.*sign(boundaryData(idx));
    
    % sometimes the x and y data are vectors
    if isvector(h.XData)
        [xdata, ydata] = meshgrid(h.XData, h.YData);
    else
        xdata = h.XData;
        ydata = h.YData;
    end
    
    idx = ~inpolygon(xdata,ydata,boundaryData(:,1),boundaryData(:,2));
    h.ZData(idx) = nan;
end

if nargout > 0
    cc = c;
    hh = h;
    bb = b;
end

end
