function [ cen_lon, cen_lat ] = usa_centroid_geo ( gen_lon, gen_lat, sample_num )

%*****************************************************************************80
%
%% usa_centroid_geo() estimates geometric centroids of Voronoi regions.
%
%  Licensing:
%
%    This code is distributed under the MIT license.
%
%  Modified:
%
%    06 December 2016
%
%  Input:
%
%    real GEN_LON(*), GEN_LAT(*), the longitude and latitude 
%    of generators to be displayed on a map of Florida.
%
%    integer SAMPLE_NUM, the number of samples to use.
%
%  Output:
%
%    real CEN_LON(*), CEN_LAT(*), the longitude and latitude 
%    of the centroids of the Voronoi regions.
%
  [ usa_lon, usa_lat ] = usa_shape_read ( );
%
%  Extract the LL bounding box for the Florida polygon.
%
  usa_lon_min = min ( usa_lon );
  usa_lon_max = max ( usa_lon );
  usa_lat_min = min ( usa_lat );
  usa_lat_max = max ( usa_lat );
%
%  Initialize the counts for each generator.
%
  n = length ( gen_lon );
  cen_lat = zeros ( n, 1 );
  cen_lon = zeros ( n, 1 );
  cen_count = zeros ( n, 1 );
%
%  Choose random LL points in the bounding box.
%  Accept them only if they are contained in the Florida polygon.
%
  for s = 1 : sample_num
    inside = 0;
    while ( ~ inside )
      lat = usa_lat_min + ( usa_lat_max - usa_lat_min ) * rand ( 1, 1 );
      lon = usa_lon_min + ( usa_lon_max - usa_lon_min ) * rand ( 1, 1 );
%
%  I ASSUME I can use LON, LAT in that order!!!
%  If so, go back and fix USA_SAMPLE too.
%
      inside = inpolygon ( lat, lon, usa_lat, usa_lon );

      if ( inside )
        d_min = Inf;
        i_min = 0;
        for i = 1 : n
          d = ( lat - gen_lat(i) ) ^ 2 + ( lon - gen_lon(i) ) ^ 2;
          if ( d < d_min )
            d_min = d;
            i_min = i;
          end
        end
        cen_lat(i_min) = cen_lat(i_min) + lat;
        cen_lon(i_min) = cen_lon(i_min) + lon;
        cen_count(i_min) = cen_count(i_min) + 1;
        break
      end

    end
  end
%
%  Set the centroid estimates.
%
  for i = 1 : n
    if ( cen_count(i) == 0 )
      cen_lon(i) = gen_lon(i);
      cen_lat(i) = gen_lat(i);
    else
      cen_lon(i) = cen_lon(i) / cen_count(i);
      cen_lat(i) = cen_lat(i) / cen_count(i);
    end
  end

  return
end

