Skip to content

gravwonbev: vertical polygon sides are counted twice #4

Description

@Advik-B

In gravwonbev, the branch for a vertical side (grav_2dpolybodies.jl line 409) returns

g = factor*x1*lor21

If you take the general expression as x21 goes to 0 (so x1 = x2, x1*z2 - x2*z1 = x1*z21 and R = z21^2), the θdiff term drops out and you're left with 0.5*factor*x1*lor21. So every exactly vertical side counts double.

My guess is the 0.5 got lost when this was adapted from the magnetic version, since in tmagwonbev it's already inside lor21 = 0.5*(log(r2)-log(r1)) while in gravwonbev it's lor21 = log(r2) - log(r1) (line 376).

Repro on b233034:

using MagGravPoly.MG2D

xz = [3.0 0.0]

rect = [-10.0 30.0; 10.0 30.0; 10.0 10.0; -10.0 10.0]
body = GravPolygBodies2D([[1, 2, 3, 4]], rect, [400.0]; ylatext=nothing)
println("talwani: ", tgravpolybodies2Dgen(xz, body, "talwani"))
println("wonbev:  ", tgravpolybodies2Dgen(xz, body, "wonbev"))

# same box with two corners nudged by 1e-9 so neither side is exactly vertical
rect2 = [-10.0 30.0; 10.0 + 1e-9 30.0; 10.0 10.0; -10.0 + 1e-9 10.0]
body2 = GravPolygBodies2D([[1, 2, 3, 4]], rect2, [400.0]; ylatext=nothing)
println("wonbev, nudged: ", tgravpolybodies2Dgen(xz, body2, "wonbev"))
talwani: [0.10323558366559354]
wonbev:  [0.1857224974488919]
wonbev, nudged: [0.10323558366580585]

So a plain rectangle comes out about 1.8x too big, and 1e-9 away from vertical it's fine. g = 0.5*factor*x1*lor21 should fix it. The magnetic x21 == 0 branch agrees with talwani on the same box, so as far as I can tell it's only gravity.

Activity

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type

    Projects

    No projects

      Milestone

      No milestone

      Relationships

      None yet

      Development

      No branches or pull requests

      Issue actions