function outputs = calcContour(function_name, fcn, numargout, varargin)
%CALCCONTOUR Calculate a contour plot of a data grid.
%
%  OUTPUTS = CALCCONTOUR(FUNCTION_NAME, FCN, NUMARGOUT, VARARGIN) displays
%  a contour plot of the data in VARARGIN.  FUNCION_NAME is the name of the
%  contour function, generally the name of the calling function. FCN is a
%  function handle to the contour function. NUMARGOUT is a scalar defining
%  the number of output arguments of the calling function. VARARGIN is a
%  cell array of the input data of the calling function. OUTPUTS is a cell
%  array with NUMARGOUT outputs.
%
%  % Example
%  ---------
%  % Contour the geoid height with 10 levels and set the edge color to red.
%  output=calcContour('contourm',@contour, 2, geoid, geoidrefvec, ...
%                     10,'edgecolor','red');
%
%  See also CONTOUR, CONTOUR3, CONTOURM, CONTOUR3M.

%  Copyright 2004 The MathWorks, Inc.
%  $Revision: 1.1.4.2 $    $Date: 2004/12/18 07:42:55 $

if numargout > 2
   eid = sprintf('%s%s%s','map:',function_name, ':tooManyOutputs');
   msg = sprintf('%s%s%s\n%s%d%s', ...
      'Function ', upper(function_name), ...
      ' expected at most 2 output arguments', ...
      'but was called instead with ', numargout, ' output arguments.');
   error(eid,msg);
end

% If no inputs, then execute as a GUI.
if numel(varargin) == 0
   contor3mui;
   switch numargout
      case 1
         outputs{1} = {};
      case 2
         outputs{1} = {};
         outputs{2} = {};
      otherwise
         outputs = {};
   end
   return
end

% Check the number of input arguments
checknargin(2,inf,numel(varargin),function_name);

% Parse the contour inputs.
[Z, lat, lon, levels, lineStyle, props] = ...
   parseContourInputs(function_name, varargin{:});

% Check if using a map axes.
ismapped = ismap(gca);
if ismapped

   %  Erase existing map if NextPlot property is replace
   nextmap;
   mstruct = gcm;

   % Create a graticule grid if lat,lon are vectors
   if isvector(lat) || isvector(lon)
     [lat,lon] = meshgrat(lat,lon);
   end

   % Check if using CONTOURFM
   if isequal(function_name, 'contourfm')

      if isfield(mstruct,'origin') && mstruct.origin(1) ~=0
         % CONTOURFM does not support a non-zero latitude origin
         eid = sprintf('%s%s%s','map:',function_name, ':latOriginNonZero');
         msg = sprintf('%s%d%s\n%s', ...
                       'The latitude coordinate of the origin vector is set at ', ...
                       mstruct.origin(1), ' degrees.', ...
                      'CONTOURF will fail using this rotated coordinate system.');
         error(eid, msg);
      end

      % Compute the forward transform.
      % Do not trim the data 
      maplatlimit = mstruct.maplatlimit;
      maplonlimit = mstruct.maplonlimit;
      mstruct.maplatlimit = [-90 90];
      mstruct.maplonlimit = [-180 180];

      % Compute the forward transform, without clipping
      % in order to sort the lat,lon, Z arrays
      [x,y,z] = mfwdtran(mstruct,lat,lon,Z,'none');

      % Sort by X
      [null, index] = sortrows(x');
      lat = lat(:,index);
      lon = lon(:,index);
      Z = Z(:,index);

      % Compute the forward transform on the sorted lat,lon, and Z arrays
      [x,y,Z,S] = mfwdtran(mstruct,lat,lon,Z,'surface');
      %x(lon >= 180) = abs(x(lon>= 180));

      % Return the mstruct to previous values
      mstruct.maplatlimit = maplatlimit;
      mstruct.maplonlimit = maplonlimit;

   else
      % Using CONTOURM or CONTOUR3M
      [x,y,null,S] = mfwdtran(mstruct,lat,lon,Z,'surface');
   end

else
   % The axes is not a map axes.
   % Set x and y from the longitude and latitude values.
   x = lon;
   y = lat;
end

% Create the argument cell array to pass to the contour function.
args = {x, y, Z, levels, lineStyle};
args = args(~cellfun('isempty',args));

% Execute the contour function
[c,h] = fcn(args{:});

% Add the MATLAB graphics properties.
if ~isempty(props)
   for i=length(props)-1:-2:1
      if ~isprop(h,props{i})
         props(i:i+1) = [];
      end
   end
end

% Set the HG properties
if ~isempty(props)
   set(h,props{:});
end

% Post execution for contour*m functions
if ismapped
   % Set the Tag to the contour level and 
   % add the ButtonDownFcn for clicking contour levels
   setUserData(h,S)

   if isequal(function_name, 'contourfm') 
      set(h,'Tag','FillContour');
      % Restack to bottom
      restack(h,'bottom');
      tightmap
   elseif isequal(function_name, 'contourm') 
      set(h,'Tag','Contour');
   else
      set(h,'Tag','Contour3d');
   end
end

% Set the output arguments and 
% inverse transform the countour matrix, if needed
switch numargout
   case 1
      outputs{1} = invtranContourMatrix(ismapped,c);
   case 2
      outputs{1} = invtranContourMatrix(ismapped, c);
      outputs{2} = h;
   otherwise
      outputs = {};
end

%--------------------------------------------------------------------------
function c = invtranContourMatrix(ismapped,c)
% Inverse transform the contour matrix from X,Y to Lat,Lon
% Inputs:  
%    ismapped - True if axes is created with AXESM
%    c        - ContourMatrix
% Outputs: 
%    c        - unprojected ContourMatrix

if ismapped

   i=0;
   startpt = 1;
   while (startpt < size(c,2))

      i=i+1;
      npoints = c(2,startpt);
      x = c(1,startpt+1:startpt+npoints)';
      y = c(2,startpt+1:startpt+npoints)';

      [lat,lon] = minvtran(x,y);
      c(1,startpt+1:startpt+npoints) = lon';
      c(2,startpt+1:startpt+npoints) = lat';
      startpt = startpt + npoints + 1;

   end
end

%-------------------------------------------------------------------------
function setUserData(h,S)
% Set the UserData and the ButtonDownFcn
% Inputs:
%    h - contourgroup handle
%    S - UserData structure

% Save and set the renderer to painters
%renderer = get(gcf,'renderer');
set(gcf,'renderer','painters');

% Set the contourgroup handle to accept button click and the ButtonDownFcn
set(h,'ButtonDownFcn','uimaptbx');
props = {'HitTest','HitTestArea'};
for i=1:length(props)
   if isprop(h,props{i})
     set(h,props{i},'on');
   end
end

% Set the childern (patch) objects with the UserData and Tag
if isequal(get(h,'type'),'hggroup')
   handle = get(h,'children');
   set(h,'UserData',S);
else
   handle = h;
end
for linecount=1:numel(handle)
   z_level = get(handle(linecount),'UserData');
   set(handle(linecount), ...
      'HitTest','on', ...
      'ButtonDownFcn','uimaptbx', ...
      'Tag',num2str(z_level))
   if handle == h
      set(handle(linecount), 'UserData',S);
   end
end

% Reset the renderer 
%set(gcf,'renderer',renderer);