function cvt_corn ( alg, it_max, nb, ni, np, ns, pr_bud, pr_bud_angular, r )

%*****************************************************************************80
%
%% cvt_corn() models the growth of a corn kernel.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    18 December 2020
%
%  Author:
%
%    John Burkardt
%
%  Input:
%
%    integer ALG: specifies the algorithm to use.
%    1 = update circle + disk together.
%    2 = update circle only, then update disk only.
%    3 = update circle only, then update circle + disk together.
%
%    integer IT_MAX: the maximum number of iterations.
%
%    integer NB, the number of boundary points.
%
%    integer NI, the number of interior points.
%
%    integer NP, the total number of points.
%
%    integer NS, the number of sample points.
%
%    real PR_BUD, the probability that a boundary point will bud.
%
%    real PR_BUD_ANGULAR, the probability that a boundary point bud
%    will be constrained to the boundary.
%
%    real R: the radius of the circle that represents the corn kernel.
%
%  Local:
%
%    integer IT: the current iteration index.
%
%    integer P_TYPE(NP): the "type" of each generator.
%    1, the point is constrained to the boundary of the disk.
%    2, the point is constrained to the disk.
%
%    real P_XY(2,NP): the location of the generators.
%
%    real R_MAX: the radius of the plot region.
%
  fprintf ( 1, '\n' );
  fprintf ( 1, 'cvt_corn():\n' );
  fprintf ( 1, '  MATLAB/Octave version %s\n', version ( ) );
  fprintf ( 1, '  Simulate the growth of a corn kernel.\n' );

  parameter_print ( alg, nb, ni, np, ns, pr_bud, pr_bud_angular, r );

  r_max = 1.5 * r;
%
%  Carry out steps 0 through IT_MAX of the iteration.
%
  for it = 0 : it_max

    if ( it == 0 )

      [ p_xy, p_type ] = state_initialize ( nb, ni, np, r );

      header = sprintf ( ...
        'CVT step %d, with %d boundary, %d interior, R = %g', it, nb, ni, r );
      figure ( 1 );
      cvt_display ( r, r_max, p_xy, p_type, header );

    else

      [ nb, ni, np, p_xy, p_type ] = bud ( pr_bud, pr_bud_angular, r, ...
        nb, ni, np, p_xy, p_type );

      [ r, p_xy, p_type ] = radial_growth ( r, np, p_xy, p_type );
      if ( r_max <= r )
        r_max = 1.5 * r;
      end

      if ( alg == 1 )
        p_xy = cvt_disk ( r, ns, np, p_xy, p_type );
        header = sprintf ( ...
          'CVT step %d, %d boundary, %d interior, R = %g', it, nb, ni, r );
        figure ( 1 );
        cvt_display ( r, r_max, p_xy, p_type, header );
        pause ( 1 );
      elseif ( alg == 2 )
        p_xy = cvt_circle ( r, np, p_xy, p_type );
        header = sprintf ( ...
          'CVT step %d, %d boundary, %d interior, R = %g', it, nb, ni, r );
        figure ( 1 );
        cvt_display ( r, r_max, p_xy, p_type, header );
        pause ( 1 );
        p_xy = cvt_disk_interior ( r, ns, np, p_xy, p_type );
        header = sprintf ( ...
          'CVT %d, %d boundary, %d interior, R = %g', it, nb, ni, r );
        figure ( 1 );
        cvt_display ( r, r_max, p_xy, p_type, header );
        pause ( 1 );
      elseif ( alg == 3 )
        p_xy = cvt_circle ( r, np, p_xy, p_type );
        header = sprintf ( ...
          'CVT step %d, %d boundary, %d interior, R = %g', it, nb, ni, r );
        figure ( 1 );
        cvt_display ( r, r_max, p_xy, p_type, header );
        pause ( 1 );
        p_xy = cvt_disk ( r, ns, np, p_xy, p_type );
        header = sprintf ( ...
          'CVT step %d, %d boundary, %d interior, R = %g', it, nb, ni, r );
        figure ( 1 );
        cvt_display ( r, r_max, p_xy, p_type, header );
        pause ( 1 );
      end

    end

    if ( false )
      figure ( 2 );
      header = sprintf ( ...
        'CVT step %d, %d boundary, %d interior, R = %g', it, nb, ni, r );
      state_display ( nb, ni, np, r, r_max, p_xy, p_type, header );
      pause ( 1 );
    end

  end

  return
end

