good_contour.f90

module good_contour_m

  implicit none

contains

  type(polyline) function good_contour(corner, step, z, level, inside_point, &
       outside_points)

    ! Finds a contour in the plane, if it exists, which is closed,
    ! contains inside_point and does not contain any of
    ! outside_points. There might be several satisfying contours, this
    ! procedure returns the first satisfying contour it finds. If such
    ! a contour does not exist, the procedure returns a null polyline.

    use contour_531, only: polyline, find_contours_reg_grid, null_polyline
    use geometry, only: polygon_contains_point

    real, intent(in):: corner(:) ! (2) [x, y] corresponding to z(1, 1)
    real, intent(in):: step(:) ! (2)

    real, intent(in):: z(:, :)
    ! We are assuming that significant values of z cannot be greater
    ! than zmax and that the missing value flag in z is greater than
    ! zmax.

    real, intent(in):: level
    real, intent(in):: inside_point(:) ! (2)
    real, intent(in):: outside_points(:, :) ! (2, n_out)

    ! Local:
    type(polyline), allocatable:: contours(:)
    logical found_good_contour
    integer i ! contour index
    integer j ! index in outside_points
    integer n_cont ! number of contours
    integer n_out ! number of outside points
    real, parameter:: zmax = 100.

    !--------------------------------------------------------------

    n_out = size(outside_points, 2)
    call find_contours_reg_grid(corner, step, z, level, contours, zmax)
    n_cont = size(contours)
    found_good_contour = .false.
    i = 1

    do while (.not. found_good_contour .and. i <= n_cont)
       if (contours(i)%closed) then
          if (polygon_contains_point(contours(i)%points, inside_point)) then
             found_good_contour = .true. ! maybe
             ! Does contours(i) contain one of outside_points?
             j = 1

             do while (found_good_contour .and. j <= n_out)
                found_good_contour &
                     = .not. polygon_contains_point(contours(i)%points, &
                     outside_points(:, j))
                j = j + 1
             end do
          end if
       end if

       i = i + 1
    end do

    if (found_good_contour) then
       good_contour = contours(i - 1)
    else
       good_contour = null_polyline()
    end if

  end function good_contour

end module good_contour_m