function value = lens_half_2d ( func, xc, yc, r, theta1, theta2, norder )

%*****************************************************************************80
%
%% lens_half_2d() approximates an integral in a circular half lens in 2D.
%
%  Discussion:
%
%    A circular half lens is formed by drawing a circular arc,
%    and joining its endpoints.
%
%    This rule for a circular half lens simply views the region as 
%    a product region, with a coordinate "S" that varies along the
%    radial direction, and a coordinate "T" that varies in the perpendicular
%    direction, and whose extent varies as a function of S.
%
%    A Gauss-Legendre rule is used to construct a product rule that is
%    applied to the region.  The accuracy of the Gauss-Legendre rule,
%    which is valid for a rectangular product rule region, does not
%    apply straightforwardly to this region, since the limits in the
%    "T" coordinate are being handled implicitly.
%
%    This is simply an application of the QMULT_2D algorithm of Stroud.
%
%  Licensing:
%
%    This code is distributed under the MIT license. 
%
%  Modified:
%
%    27 November 2004
%
%  Author:
%
%    John Burkardt
%
%  Reference:
%
%    Arthur H Stroud,
%    Approximate Calculation of Multiple Integrals,
%    Prentice Hall, 1971,
%    QA311.S85
%
%  Input:
%
%    external FUNC, the name of the user supplied function of two
%    variables which is to be integrated, of the form:
%      function value = func ( x, y )
%
%    real XC, YC, the center of the circle.
%
%    real R, the radius of the circle.
%
%    real THETA1, THETA2, the angles of the rays
%    that begin and end the arc.
%
%    integer NORDER, the order of the Gauss-Legendre rule
%    to be used.  Legal values include 1 through 20, 32 or 64.
%
%  Output:
%
%    real VALUE, the approximate value
%    of the integral of the function over the half lens.
%

%
%  Retrieve the Legendre rule of the given order.
%
  [ xtab, weight ] = legendre_set ( norder );
%
%  Determine the points A (on the secant) and B (on the circumference)
%  that will form the "S" direction.
%
  ax = xc + r * 0.5 * ( cos ( theta1 ) + cos ( theta2 ) );
  ay = yc + r * 0.5 * ( sin ( theta1 ) + sin ( theta2 ) );
  bx = xc + r * cos ( 0.5 * ( theta1 + theta2 ) );
  by = yc + r * sin ( 0.5 * ( theta1 + theta2 ) );
%
%  Find the length of the line between A and B.
%
  s_length = sqrt ( ( ax - bx )^2 + ( ay - by )^2 );

  if ( s_length == 0.0 )
    value = 0.0;
    return
  end
%
%  Determine the unit vector in the T direction.
%
  tdirx = ( ay - by ) / s_length;
  tdiry = ( bx - ax ) / s_length;

  quad = 0.0;

  for i = 1 : norder

    w1 = 0.5 * s_length * weight(i);
%
%  Map the quadrature point to an S coordinate.
%
    sx = ( ( 1.0 - xtab(i) ) * ax + ( 1.0 + xtab(i) ) * bx ) / 2.0;
    sy = ( ( 1.0 - xtab(i) ) * ay + ( 1.0 + xtab(i) ) * by ) / 2.0;
%
%  Determine the length of the line in the T direction, from the
%  S axis to the circle circumference.
%
    thi = sqrt ( ( r - 0.25 * ( 1.0 - xtab(i) ) * s_length ) ...
                         *    ( 1.0 - xtab(i) ) * s_length );
% 
%  Determine the maximum and minimum T coordinates by going
%  up and down in the T direction from the S axis.
%
    cx = sx + tdirx * thi;
    cy = sy + tdiry * thi;
    dx = sx - tdirx * thi;
    dy = sy - tdiry * thi;
%
%  Find the length of the T direction.
%
    t_length = sqrt ( ( cx - dx )^2 + ( cy - dy )^2 );

    for j = 1 : norder

      w2 = 0.5 * t_length * weight(j);
%
%  Map the quadrature point to a T coordinate.
%
      tx = ( ( 1.0 - xtab(j) ) * cx + ( 1.0 + xtab(j) ) * dx ) / 2.0;
      ty = ( ( 1.0 - xtab(j) ) * cy + ( 1.0 + xtab(j) ) * dy ) / 2.0;

      quad = quad + w1 * w2 * feval ( func, tx, ty );

    end

  end

  value = quad;

  return
end
