function [m,ok] = leaf_iterate( m, varargin )
%[m,ok] = leaf_iterate( m, numsteps, ... )
%   Run the given number of iterations of the growth process.
%   In each iteration, the following things happen:
%       * Strains are set up in the leaf according to the effects of the
%         morphogens at each point causing the material to grow.
%       * The elasticity problem is solved to produce a new shape for the
%         leaf.
%       * Morphogens which have a nonzero diffusion coefficient are allowed 
%         to diffuse through the mesh.
%       * If dilution by growth is enabled, the amount of every morphogen
%         is diluted at each point according to how much that part of the
%         leaf grew.
%       * A user-supplied routine is invoked to model local interactions
%         between the morphogens.
%       * The layer of biological cells is updated.
%       * If requested by the options, the mesh is plotted after each
%         iteration.
%
%   Results:
%   ok: True if there was no problem (e.g. invalid arguments, a user
%   interrupt, or an error in the interaction function).  False if there
%   was.
%
%   Arguments:
%   numsteps: The number of iterations to perform.  If this is zero, the
%   computation will continue indefinitely or until terminated by the
%   'until' or 'targetarea' options.  The default is 1.
%
%   Options:
%   'until'       Run the simulation until this time has been reached or
%                 exceeded.  A value of zero disables this option.
%   'targetarea'  Run the simulation until the area of the canvas is at
%                 least this number times the initial area.  A value of
%                 zero disables this option.
%   'plot'  An integer n.  The mesh will be plotted after every n
%           iterations.  0 means do not plot the mesh at all; -1 means plot
%           the mesh only after all the iterations have been completed.
%           The default value is 1.
%   Example:
%       m = leaf_iterate( m, 20, 'plot', 4 );
%
%   Equivalent GUI operation: clicking one of the following buttons in the
%   "Simulation" panel: "Run" (do a specified number of steps), "Step" (do
%   one step), or "Run to..." (run until the area has increased by a
%   specified factor).
%
%   Topics: Simulation.

%m = recordframe( m );
%beep;
%return;

% Steps in each cycle:
%   * Invoke interaction function.
%   * Perform diffusion (taking account of clamped values).
%   * Perform decay.  (This would be better built into the diffusion
%         calculation.)
%   * Restore clamped values.
%   * Calculate the polarisation gradients.
%   * Perform a FEM calculation.
%   * Generate the FEM cell normal vectors.
%   * Split large FEM cells.
%   * Recalculate FEM cell areas.
%   * Recalculate the second layer global positions and split large
%     second layer cells.

    if isempty(m), return; end
    [ok, numsteps, args] = getTypedArg( mfilename(), 'numeric', varargin );
    if ~ok, return; end
    numsteps = floor(numsteps);
    [s,ok] = safemakestruct( mfilename(), args );
    if ~ok, return; end
    s = defaultfields( s, 'plot', 1, 'until', 0, 'targetarea', 0 );
    ok = checkType( mfilename(), 'numeric', s.plot );
    if ~ok, return; end
    s.plot = floor(s.plot);
    ok = checkcommandargs( mfilename(), s, 'only', ...
            'plot', 'until', 'targetarea', 'handles' );
    if ~ok, return; end
    haveHandles = isfield( s, 'handles' );

    recordedStep = m.globalDynamicProps.currentIter;
    
    step = 1;
    interrupted = false;
    m.stop = false;
    plotted = false;
    while ok && ~interrupted && ~finished( m, step, numsteps, s.until, s.targetarea )
        fprintf( 1, 'Starting iteration %d, time %.3f, at %s.\n', ...
        	m.globalDynamicProps.currentIter+1, m.globalDynamicProps.currenttime, datestring() );
        startTime = cputime;
        tic;
        if ~m.globalProps.flatten
            if m.globalProps.allowInteraction
                if haveHandles
                    [m,ok] = attemptInteractionFunction( m, s.handles );
                  % updateGUIFromMesh( s.handles );
                    updateGUIFromMesh( m );
                else
                    [m,ok] = attemptInteractionFunction( m );
                end
                if ~ok
                    break;
                end
            else
                fprintf( 1, 'Interaction function disabled.\n' );
            end
        end
        m = disallowNegativeGrowth( m );
        if ~m.globalProps.flatten
            if m.globalProps.diffusionEnabled
                m = diffusegrowth( m );
                nondiffusibles = find( all( m.absKvector(:,:) == 0, 1 ) );
            else
                nondiffusibles = 1:size(m.morphogens,2);
            end
            for i=nondiffusibles(:)'
                nonclamped = m.morphogenclamp(:,i) < 1;
                m.morphogens(nonclamped,i) = m.morphogens(nonclamped,i) ...
                    + m.globalProps.timestep * ( ...
                        m.mgen_production(nonclamped,i) ...
                        - m.morphogens(nonclamped,i) * m.mgen_absorption(i) ...
                      );
            end
        end
        if m.globalProps.performinternalrotation
            badfixedNodes = find( m.fixedDFmap(:,1) ~= m.fixedDFmap(:,2) );
            if isempty(badfixedNodes)
                m = stepinternalrotation( m );
            else
                % Can't rotate.
                fprintf( 1, [ 'leaf_iterate: cannot perform internal rotation due to nodes\n', ...
                    '    fixed in X but not Y, or Y but not X:\n    ' ] );
                fprintf( 1, ' %d', badfixedNodes );
                fprintf( 1, '\n' );
                m = derotate( m );
                m.globalProps.performinternalrotation = false;
                if haveHandles
                    set( s.handles.performinternalrotation, 'Value', false );
                end
            end
        end
        if m.stop
            ok = false;
            break;
        end
        [m,ok] = onestep( m, m.globalProps.useGrowthTensors );
        m.globalDynamicProps.currentIter = m.globalDynamicProps.currentIter + 1;
        m.globalDynamicProps.currenttime = ...
            m.globalDynamicProps.currenttime + m.globalProps.timestep;
        if haveHandles
            if m.globalDynamicProps.currentIter==1
                enableMutations( s.handles, 'off' );
            end
            if testAndClear( s.handles.plotFlag )
                fprintf( 1, 'Change (plot) flag detected during iteration %d.\n', ...
                    m.globalDynamicProps.currentIter );
                m = getPlotOptionsFromDialog( m, s.handles );
            end
            c = get( s.handles.commandFlag, 'UserData' );
            set( s.handles.commandFlag, 'UserData', struct([]) );
            if ~isempty(c)
                fprintf( 1, 'User commands detected during iteration %d.\n', ...
                    m.globalDynamicProps.currentIter );
                if recordedStep < m.globalDynamicProps.currentIter
                    m = recordcommand( m, false, 'iterate', ...
                            m.globalDynamicProps.currentIter - recordedStep, ...
                            'plot', s.plot );
                    recordedStep = m.globalDynamicProps.currentIter;
                end
                m = executeCommands( m, c, false, s.handles );
            end
            if testAndClear( s.handles.stopButton )
                interrupted = true;
            end
        end
        if m.globalProps.performinternalrotation && (~m.globalProps.showinternalrotation)
        	m = derotate( m );
        end
        m.timeForIter = cputime() - startTime;
        m.ticForIter = toc;
        m = calculateOutputs( m );
        if haveHandles
            s.handles.mesh = m;
            guidata( s.handles.output, s.handles );
            announceSimStatus( s.handles );
        end
        if (s.plot > 0) && (mod(step,s.plot)==0)
            if haveHandles
                s.handles = processPendingData( s.handles );
                m = s.handles.mesh;
            else
                m = leaf_plot( m );
            end
            wasMakingMovie = m.globalProps.makemovie;
            m = recordframe( m );
            if haveHandles && (wasMakingMovie && ~m.globalProps.makemovie)
                % Movie was closed for some reason.  Reset the label on
                % the movieButton.
                set( s.handles.movieButton, 'String', 'Record movie...' );
            end
            plotted = true;
        else
            plotted = false;
        end
        step = step+1;
        fprintf( 1, 'Completed iteration %d at %s.\n', ...
        	m.globalDynamicProps.currentIter, datestring() );
    end
    if recordedStep < m.globalDynamicProps.currentIter
        % Some iterations have been performed since an iterate command was
        % added to the script.
        if interrupted || (numsteps > 0)
            % Record a command to iterate for a number of steps.
            m = recordcommand( m, false, 'iterate', ...
                    m.globalDynamicProps.currentIter - recordedStep, ...
                    'plot', s.plot );
        else
            % Record a command to iterate until a target area is reached.
            m = recordcommand( m, false, 'iterate', ...
                    0, ...
                    'targetarea', s.targetarea, ...
                    'plot', s.plot );
        end
    end
    if m.globalProps.maxIters ~= m.globalDynamicProps.currentIter
        m.globalProps.maxIters = m.globalDynamicProps.currentIter;
        if haveHandles
            announceSimStatus( s.handles, m );
        end
    end
    if (plotted == false) && ...
            ((s.plot == -1) || ((s.plot > 0) && (mod(step,s.plot) ~= 0)))
        if haveHandles && testAndClear( s.handles.plotFlag )
            m = getPlotOptionsFromDialog( m, s.handles );
        end
        m = leaf_plot( m );
        m = recordframe( m );
    end
    updateGUIFromMesh( m );  % In case the interaction function changed any
                             % of the mesh properties that are displayed in
                             % the GUI.  For efficiency, we could do this
                             % only if the i.f. notifies us in some way,
                             % but since updateGUIFromMesh is just setting
                             % the contents of a couple of dozen text boxes
                             % and checkboxes, it's not worth it.
                             % Note that updateGUIFromMesh is safe to call
                             % even if this function is not being called
                             % from within GFtbox.
    ok = ok && ~interrupted;
end

function f = finished( m, step, numsteps, targettime, targetarea )
    if isfield( m, 'stop' ) && m.stop
        fprintf( 1, 'Simulation terminated by interaction function.\n' );
        f = true;
        return;
    end
    if (numsteps > 0) && (step > numsteps)
        fprintf( 1, 'Simulation terminated after %d steps.\n', numsteps );
        f = true;
        return;
    end
    if (targettime > 0) && ((targettime - m.globalDynamicProps.currenttime) <= m.globalProps.timestep/2)
        fprintf( 1, 'Simulation terminated on reaching time %f.\n', targettime );
        f = true;
        return;
    end
    if (targetarea > 0) && (m.globalDynamicProps.currentArea/m.globalProps.initialArea >= targetarea)
        fprintf( 1, 'Simulation terminated on reaching area multiple %f.\n', targetarea );
        f = true;
        return;
    end
    f = false;
end

function isSet = checkFlag( s, flag )
   isSet = isfield( s.handles, flag ) ...
           && (get( s.handles.(flag), 'Value' ) ~= 0);
   if isSet
       set( s.handles.(flag), 'Value', 0 )
   end
end
