anchor.diff:                              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
  +                        if ((x0 == x1) .or. (dble(x) <= x_line)) then
                               inside = .not. inside
                               call check_vertex_hits(y,y0,y1,inside, top_vertex, bottom_vertex, internal_precision)
                           else
  
positive.code:         real,intent(in)::lat,lon
          type(position),intent(in)::pos
          integer,intent(in) :: nx,ny
          integer :: i
  
          ! test the lower left quadrant
          find_surrounding%x = [pos%x, pos%x-1, pos%x-1, pos%x  ]
          find_surrounding%y = [pos%y, pos%y,   pos%y-1, pos%y-1]
  
          if (test_surrounding(lat, lon, lo, find_surrounding, nx, ny, 4)) then
              ! rearrange the order consistent with what bilin_weights is expecting
              find_surrounding%x = [pos%x, pos%x-1, pos%x, pos%x-1]
              return
          endif
  
  
          ! test the upper left quadrant
          find_surrounding%x = [pos%x, pos%x-1, pos%x-1, pos%x  ]
          find_surrounding%y = [pos%y, pos%y,   pos%y+1, pos%y+1]
  
          if (test_surrounding(lat, lon, lo, find_surrounding, nx, ny, 4)) then
              ! rearrange the order consistent with what bilin_weights is expecting
              find_surrounding%x = [pos%x, pos%x-1, pos%x, pos%x-1]
              return
          endif
  
  
          ! test the upper right quadrant
          find_surrounding%x = [pos%x, pos%x+1, pos%x+1, pos%x  ]
          find_surrounding%y = [pos%y, pos%y,   pos%y+1, pos%y+1]
  
          if (test_surrounding(lat, lon, lo, find_surrounding, nx, ny, 4)) then
              ! rearrange the order consistent with what bilin_weights is expecting
              find_surrounding%x = [pos%x, pos%x+1, pos%x, pos%x+1]
              return
          endif
  
  
          ! test the lower right quadrant
          find_surrounding%x = [pos%x, pos%x+1, pos%x+1, pos%x  ]
          find_surrounding%y = [pos%y, pos%y,   pos%y-1, pos%y-1]
  
          if (test_surrounding(lat, lon, lo, find_surrounding, nx, ny, 4)) then
              ! rearrange the order consistent with what bilin_weights is expecting
              find_surrounding%x = [pos%x, pos%x+1, pos%x, pos%x+1]
              return
          endif
  
          write(*,*) "ERROR: Failed to find point", lon, lat, " in a quadrant near:"
          write(*,*) lo%lon(pos%x, pos%y), lo%lat(pos%x, pos%y)
  
negative.code:     !! @retval  real    weights
      !!
      !!------------------------------------------------------------
      function tri_weights(yi,y,xi,x)
          implicit none
          real,intent(in)::yi,y(4),xi,x(4)
          real::x0,y0
          real, dimension(4) ::tri_weights
  
          real :: w1,w2,w3,denom
  
          x0 = sum(x)/4
          y0 = sum(y)/4
  
          associate(y1 => y(1), y2 => (y(2)), y3 => y0, &
                    x1 => x(1), x2 => (x(2)), x3 => x0)
  
          denom = ((y2-y3) * (x1-x3) + (x3-x2) * (y1-y3))
          if (denom==0) then
              write(*,*) "ERROR: Triangular interpolation broken"
              write(*,*) xi, yi
              write(*,*) "Triangle vertices"
              write(*,*) x1, y1
              write(*,*) x2, y2
              write(*,*) x3, y3
              stop "Denominator in triangulation is broken"
          endif
  
          ! This is the core of the algorithm.
          ! compute weights for each x,y point
          w1 = ((y2-y3) * (xi-x3) + (x3-x2) * (yi-y3)) &
                             / denom
  
          w2 = ((y3-y1) * (xi-x3) + (x1-x3) * (yi-y3)) &
                             / denom
  
          w3 = 1 - w1 - w2
  
          if (minval([w1,w2,w3]) < -1e-4) then
              write(*,*) "ERROR: Point not located in bounding triangle"
              write(*,*) xi, yi
              write(*,*) "Triangle vertices"
              write(*,*) x1, y1
              write(*,*) x2, y2
              write(*,*) x3, y3
              write(*,*) w1, w2, w3
              stop "Triangulation is broken"
          endif
          w1=max(0.,w1); w2=max(0.,w2); w3=max(0.,w3);
  
          if (abs((w1+w2+w3)-1)>1e-4) then
              write(*,*) "Error, w1+w2+w3 != 1"
              write(*,*) w1,w2,w3
              write(*,*) "Point"
              write(*,*) xi, yi
              write(*,*) "Triangle vertices"
              write(*,*) x1, y1
              write(*,*) x2, y2
              write(*,*) x3, y3
          endif
          denom = 1.0 / (w1+w2+w3)
          w1 = w1 * denom
          w2 = w2 * denom
          w3 = w3 * denom
  
          end associate
  
          tri_weights = [w1, w2, w3, 0]
  
      end function tri_weights
  
      !>------------------------------------------------------------
      !! Calculate the weights to use for inverse distance interpolation between surrounding points x, y
      !!
      !! Takes a point (yi,xi) and a set of 4 surrounding points (y[4],x[4]) as input
      !!
      !! @param   real    yi  input y location to interpolate to.
      !! @param   real    y   Array of surrounding y-locations to interpolate from.
      !! @param   real    xi  input x location to interpolate to.
      !! @param   real    x   Array of surrounding x-locations to interpolate from.
      !! @retval  real    weights
      !!
      !!------------------------------------------------------------
      function idw_weights(yi,y,xi,x)
          implicit none
          real,intent(in)::yi,y(:),xi,x(:)
  