anchor.diff:        call geo2xyz(src % lat, src % lng, x1(1), x1(2), x1(3))
         call geo2xyz(      lat,       lng, x2(1), x2(2), x2(3))
   
  -      call cross(x1, x2, xi)
  -      call cross(xi / sqrt(sum(xi**2)), x1, xw)
  +      xi = cross(x1, x2)
  +      xi = xi / sqrt(sum(xi**2))
  +      xw = cross(xi, x1)
   
         interpolate_terrain: do i = 1, nn
           xi(:) = x1 * cos(t(i) * adist) + xw * sin(t(i) * adist)
  
positive.code: 
    contains
  
      pure subroutine cross(x,y,z)
        real(dp), dimension(3), intent(in) :: x,y
        real(dp), dimension(3), intent(out) :: z
        z(1) =   x(2)*y(3) - x(3)*y(2)
        z(2) = - x(1)*y(3) + x(3)*y(1)
        z(3) =   x(1)*y(2) - x(2)*y(1)
      end subroutine
  
      elemental real(dp) function f(t, k)
        real(dp), intent(in) :: t, k
  
negative.code: 
      trace_ray: block
  
        real(dp), dimension(nn) :: t, lat_ray, lng_ray, h_ray, hterr
        real(dp), dimension(3) :: x1, x2, xi, xw, xr
        integer :: i, j
  
  