function [ n, x, y, w, s ] = hexagon_rule_lyness13 ( )

%*****************************************************************************80
%
%% hexagon_rule_lyness13 sets the Lyness hexagon quadrature rule #13.
%
%  Licensing:
%
%    This code is distributed under the GNU LGPL license. 
%
%  Modified:
%
%    22 April 2019
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    James Lyness, Giovanni Monegato,
%    Quadrature rules for regions having regular hexagonal symmetry,
%    SIAM Journal on Numerical Analysis.
%    Volume 14 Number 2, April 1977, pages 283-295.
%
%  Parameters:
%
%    Output, integer n, the number of points.
%
%    Output, real x(n), y(n), the points.
%
%    Output, real w(n), the weights.
%
%    Output, integer s, the strength.
%
  n = 48;
  x = zeros ( n, 1 );
  y = zeros ( n, 1 );
  w = zeros ( n, 1 );
  
  x(1) = 0.3918730038483139;
  y(1) = 0.0000000000000000;
  w(1) = 0.2764602580376149 / 6.0;
  for i = 2 : 6
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  x(7) = 0.08082128711904690;
  y(7) = 0.0000000000000000;
  w(7) = 0.06018247684485211 / 6.0;
  for i = 8 : 12
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  x(13) = 0.8673936568297751;
  y(13) = 0.0000000000000000;
  w(13) = 0.03162712410615477 / 6.0;
  for i = 14 : 18
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  x(19) = 0.7081424125773276;
  y(19) = 0.0000000000000000;
  w(19) = 0.1865934847620135 / 6.0;
  for i = 20 : 24
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  x(25) = 0.8371347439666993;
  y(25) = 0.4833199697771623;
  w(25) = 0.02136506644258773 / 6.0;
  for i = 26 : 30
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  x(31) = 0.5474220549086223;
  y(31) = 0.3160542707618311;
  w(31) = 0.2187296029611373 / 6.0;
  for i = 32 : 36
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  x(37) = 0.7991748327247630;
  y(37) = 0.2974218977320152;
  w(37) = 0.2050419868456396 / 12.0;
  for i = 38 : 42
    x(i) = cosd ( 60.0 ) * x(i-1) - sind ( 60.0 ) * y(i-1);
    y(i) = sind ( 60.0 ) * x(i-1) + cosd ( 60.0 ) * y(i-1);
    w(i) = w(i-1);
  end

  for i = 43 : 48
    x(i) =  x(i-6);
    y(i) = -y(i-6);
    w(i) =  w(i-6);
  end

  s = 15;

  return
end
