anchor.diff:      !! Note additional discussion here: https://wrf.ecse.rpi.edu//Research/Short_Notes/pnpoly.html
       !!
       !!------------------------------------------------------------
  -    function point_in_poly(x, y, poly) result(inside)
  +    function point_in_poly(x, y, poly, precision) result(inside)
           implicit none
           real, intent(in) :: x, y
           real, intent(in) :: poly(:,:)
  +        real, intent(in), optional :: precision
           logical :: inside
  -
  +        logical :: top_vertex, bottom_vertex
  +
  +        real :: internal_precision
           integer :: n, i
           real :: x0,y0,x1,y1
           double precision :: slope, x_line
   
  +        internal_precision = 1e-7
  +        if (present(precision)) internal_precision = precision
           n = size(poly,2)
   
           ! by default we assume we are not in the polygon
           inside = .False.
  +        top_vertex = .False.
  +        bottom_vertex = .False.
   
           ! check if point is a vertex
           do i=1,n
  
positive.code:                 if (y <= max(y0,y1)) then
                      if (x <= max(x0,x1)) then
  
                          if (point_is_on_line(x,y,x0,y0,x1,y1)) then
                              inside=.True.
                              return
                          endif
  
negative.code:                             slope = (dble(x1)-x0)/(dble(y1)-y0)
                              x_line = (y-dble(y0)) * slope + x0
                          else
                              x_line = -1e10 ! only inside if x0==x1
                          endif
                          ! if we cross one line then we are inside, but if we cross two we are outside, and so on
                          if ((x0 == x1) .or. (x <= x_line)) then
                              inside = .not. inside
                              call check_vertex_hits(y,y0,y1,inside, top_vertex, bottom_vertex, internal_precision)
                          else
                              top_vertex = .False.
                              bottom_vertex = .False.
                          endif
                      endif
                  endif
  