function m = leaf_lobes( m, varargin )
%m = leaf_lobes( m, ... )
%   Create a new mesh in the form of one or more lobes joined together in a
%   row.  A lobe is a semicircle on top of a rectangle.
%
%   Arguments:
%       M is either empty or an existing mesh.  If it is empty, then an
%       entirely new mesh is created, with the default set of morphogens
%       If M is an existing mesh, then its geometry is replaced by the new
%       mesh.  It retains the same set of morphogens (all set to zero
%       everywhere on the new mesh), interaction function, and all other
%       properties not depending on the specific geometry of the mesh.
%
%   Options:
%       'radius'        The radius of the semicircle.  Default 1.
%       'rings'         The number of circular rings of triangles to divide
%                       it into. Default 4.
%       'height'        The height of the rectangle, as a multiple of the
%                       semicircle's diameter.  Default 0.7.
%       'strips'        The number of strips of triangles to divide the
%                       rectangular part into.  If 0 (the default), this will
%                       be calculated from the height so as to make the
%                       triangles similar in size to those in the lobes.
%       'lobes'         The number of lobes.  The default is 1.
%       'base'          Half the number of divisions along the base of a
%                       lobe.  Defaults to rings.
%       'cylinder'      The series of lobes is to behave as if wrapped
%                       round a cylinder and the two ends stitched
%                       together.  This is implemented by constraining the
%                       nodes on the outer edges in such a way that the
%                       outer edges remain parallel to the y axis.
%
%   Example:
%       m = leaf_lobes( 'radius', 2, 'rings', 4, 'lobes', 3, 'base', 2 );
%   See also: LEAF_CIRCLE, LEAF_CYLINDER, LEAF_LUNE, LEAF_ONECELL,
%   LEAF_RECTANGLE.
%
%   Equivalent GUI operation: selecting "Lobes" in the pulldown menu in the
%   "Mesh editor" panel and clicking the "Generate mesh" button.
%
%   Topics: Mesh creation.

    [s,ok] = safemakestruct( mfilename(), varargin );
    if ~ok, return; end
    setGlobals();
    global gMAX_MGENVERSION
    s = defaultfields( s, ...
        'radius', 2, 'rings', 4, 'height', 0.7, 'strips', 0, 'lobes', 1, 'cylinder', 0, ...
        'base', 0, 'version', gMAX_MGENVERSION );
      % 'fixright', false, 'fixfoot', false );
    ok = checkcommandargs( mfilename(), s, ...
        'exact', 'radius', 'rings', 'height', 'strips', 'lobes', 'base', 'cylinder', ...
        'version' );
    if ~ok, return; end
    
    if s.strips <= 0
        s.strips = max(1,ceil(s.rings*s.height/s.radius));
    end
    if s.lobes==1
        newm = lobemesh( s.radius, s.rings, s.height, s.strips, s.base ); % , s.fixright, s.fixfoot );
    else
        newm = lobemesh( s.radius, s.rings, s.height, s.strips, s.base ); % , s.fixright, s.fixfoot );
        leftBorder = newm.borders.left;
        rightBorder = newm.borders.right;
        singleLobe = lobemesh( s.radius, s.rings, s.height, s.strips, s.base ); % , s.fixright, s.fixfoot );
        for i=2:s.lobes-1
            newm = stitchmeshes( singleLobe, newm, rightBorder, leftBorder );
        end
        [newm,renumber] = stitchmeshes( singleLobe, newm, rightBorder, leftBorder );
    end
    newm = centremesh( newm );
    m = setmeshfromnodes( newm, m, s.version );
    
    % Now set the fixed degrees of freedom.  The nodes along the bottom
    % edge of the rectangle should be fixed in the y dimension.  The
    % nodes along the right-hand side (positive x) of the rectangular part
    % should be fixed in the x dimension.
    
    totalnumnodes = size(m.nodes);

    if false
        if isfield(s,'fixright') && s.fixright
            % The nodes on the right edge are the nth node in the semicircular
            % part, and s.strips nodes at intervals of 2nrings+1, ending at the last node.
            rightedge = [ s.nrings; ((totalnumnodes - (s.strips-1)*(s.nrings+s.nrings+1)) : (s.nrings+s.nrings+1) : totalnumnodes)' ];
            m = leaf_constrain( m, 'vertexes', rightedge, 'dfs', 'x' );
        end
    end
    
    if s.cylinder
        fprintf( 1, '%s: Option ''cylinder'' not implemented.\n', mfilename() );
    end

    ymin = min(m.nodes(:,2));
    ymax = max(m.nodes(:,2));
    delta = (ymax-ymin)/1000;
    bottomedge = find( m.nodes(:,2) < ymin+delta );
    m = leaf_fix_vertex( m, 'vertex', bottomedge, 'dfs', 'y' );
    m.meshparams = s;
    m.meshparams.randomness = 0;
    m.meshparams.type = regexprep( mfilename(), '^leaf_', '' );
end
