wonbev: wrong result when the station is exactly level with a polygon corner

[Advik-B]

Hi, I've been comparing the 2D formulations in MG2D against each other and found a case where wonbev gives a wrong answer, for both magnetics and gravity: when the observation point has exactly the same z as one of the polygon's corners.

Here's a small repro on b233034:

using MagGravPoly.MG2D

tri = [20.0 120.0; 80.0 120.0; 50.0 80.0]
gbody = GravPolygBodies2D([[1, 2, 3]], tri, [400.0]; ylatext=nothing)
jind = MagnetizVector(mod=[0.76], Ideg=[35.0], Ddeg=[-8.0])
jrem = MagnetizVector(mod=[0.0], Ideg=[0.0], Ddeg=[0.0])
mbody = MagPolygBodies2D([[1, 2, 3]], tri, jind, jrem; ylatext=nothing)

for z in (80.0, 80.0 + 1e-9)
    xz = [30.0 z]
    println("z = $z")
    println("  mag  talwani: ", tmagpolybodies2Dgen(xz, 70.0, mbody, "talwani"), "  wonbev: ", tmagpolybodies2Dgen(xz, 70.0, mbody, "wonbev"))
    println("  grav talwani: ", tgravpolybodies2Dgen(xz, gbody, "talwani"), "  wonbev: ", tgravpolybodies2Dgen(xz, gbody, "wonbev"))
end

Output:

z = 80.0
  mag  talwani: [35.139302018304384]  wonbev: [214.2656270859927]
  grav talwani: [0.1461228781933853]  wonbev: [-0.4980119881727003]
z = 80.000000001
  mag  talwani: [35.139302018122024]  wonbev: [35.13930201812213]
  grav talwani: [0.14612287819432265]  wonbev: [0.14612287819432263]

If I move the station 1e-9 up or down they agree again to ~12 digits. At exactly z = 80 you also get the "A polygon side is too close to an observation point" warning, even though the station is 20 m away from the triangle.

I think it comes from the crossing check in tmagwonbev (mag_2dpolybodies.jl line 835) and gravwonbev (grav_2dpolybodies.jl line 383):

if sign(z1) != sign(z2)

With z1 exactly 0, sign(z1) is 0, so a side that only touches the station's horizontal at its end counts as a crossing. For the side from (50, 80) to (20, 120), relative to the station that's x1 = 20, z1 = 0, x2 = -10, z2 = 40, so test = x1*z2 - x2*z1 = 800 > 0, z1 >= 0 passes and θ2 gets 2π added. The side never crosses the negative x axis, so it shouldn't wrap. the other side that ends on that corner gets the same treatment on θ1.

Something like if (z1 < 0) != (z2 < 0) might be enough, or taking the angle directly as atan(x1*z2 - x2*z1, x1*x2 + z1*z2) avoids the branch cut entirely. I haven't run either against your tests though.

Thanks for putting the package out, it's been really handy for cross-checking.