function point_in_polygon ( n, x, y, x0, y0 )

!*****************************************************************************80
!
!! point_in_polygon() determines if a point is inside a polygon
!
!  Discussion:
!
!    If the points ( x(i), y(i) ) ( i = 1, 2, ..., n ) are,
!    in this cyclic order, the vertices of a simple closed polygon and
!    (x0,y0) is a point not on any side of the polygon, then the
!    procedure determines, by setting "point_in_polygon" to TRUE or FALSE,
!    whether (x0,y0) lies in the interior of the polygon.
!
!  Licensing:
!
!    This code is distributed under the MIT license. 
!
!  Modified:
!
!    07 November 2016
!
!  Author:
!
!    Original Fortran77 code by Moshe Shimrat;
!    This version by John Burkardt.
!
!  Reference:
!
!    Moshe Shimrat,
!    ACM Algorithm 112,
!    Position of Point Relative to Polygon,
!    Communications of the ACM,
!    Volume 5, Number 8, page 434, August 1962.
!
!    Richard Hacker,
!    Certification of Algorithm 112,
!    Communications of the ACM,
!    Volume 5, Number 12, page  606, December 1962.
!
!  Parameters:
!
!    Input, integer N, the number of nodes or vertices in 
!    the polygon.  N must be at least 3.
!
!    Input, real ( kind = rk ) V(2,N), the vertices of the polygon.
!
!    Input, real ( kind = rk ) P(2), the coordinates of the point to be tested.
!
!    Output, logical INSIDE, is TRUE if the point is
!    inside the polygon.
!
  implicit none

  integer, parameter :: rk = kind ( 1.0D+00 )

  integer n

  logical b
  integer i
  integer ip1
  logical point_in_polygon
  real ( kind = rk ) t
  real ( kind = rk ) x(n)
  real ( kind = rk ) x0
  real ( kind = rk ) y(n)
  real ( kind = rk ) y0
  
  b = .false.

  do i = 1, n
    ip1 = mod ( i, n ) + 1
!
!   if ( ( y(ip1) < y0 .and. y0 <= y(i)   ) .or. &
!        ( y(i)   < y0 .and. y0 <= y(ip1) ) ) then
!
    if ( y(ip1) < y0 .eqv. y0 <= y(i) ) then
      t = x0 - x(i) - ( y0 - y(i) ) * ( x(ip1) - x(i) ) / ( y(ip1) - y(i) )
      if ( t < 0.0D+00 ) then
        b = .not. b
      end if
    end if
  end do

  point_in_polygon = b

  return
end