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

%*****************************************************************************80
%
%% hexagon_lyness_rule12() sets the Lyness hexagon quadrature rule #12.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    21 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.
%
%  Output:
%
%    integer n, the number of points.
%
%    real x(n), y(n), the points.
%
%    real w(n), the weights.
%
%    integer s, the strength.
%
  n = 37;
  x = zeros ( n, 1 );
  y = zeros ( n, 1 );
  w = zeros ( n, 1 );
  
  x(1) = 0.0;
  y(1) = 0.0;
  w(1) = 0.06086889959752326;

  x(2) = 0.6537326298921177;
  y(2) = 0.0000000000000000;
  w(2) = 0.2183328883521238 / 6.0;
  for i = 3 : 7
    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(8) = 0.8075955374019075;
  y(8) = 0.4662655009153318;
  w(8) = 0.05842604367827636 / 6.0;
  for i = 9 : 13
    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(14) = 0.3558364317161740;
  y(14) = 0.2054422596388089;
  w(14) = 0.3104788608398542 / 6.0;
  for i = 15 : 19
    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(20) = 0.6365561973937096;
  y(20) = 0.3675158919195828;
  w(20) = 0.1811126315420415 / 6.0;
  for i = 21 : 25
    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(26) = 0.8258707492247416;
  y(26) = 0.1842928212756751;
  w(26) = 0.1707806759901808 / 12.0;
  for i = 27 : 31
    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 = 32 : 37
    x(i) =  x(i-6);
    y(i) = -y(i-6);
    w(i) =  w(i-6);
  end

  s = 13;

  return
end
