11  Implementation of the TALWANI approach

We define a Python class Polygon which takes as argument a list of coordinates and a parameter. The purpose of this class is twofold:

• it stores all vertices an a convenient manner
• it sorts the vertices in CCW order.

class Polygon:
    def __init__(self, x, z, density):
        xc = np.mean(x)
        zc = np.mean(z)
        angles = np.array([np.arctan2(p[1] - zc, p[0] - xc) for p in zip(x, z)])
        ind = np.argsort(angles) # if reverse order is required, flip sign of angles 
        x = x[ind]
        z = z[ind]
        x = np.append(x, x[0])
        z = np.append(z, z[0])
        x[isclose(x, 0.0)] = 0.0
        z[isclose(z, 0.0)] = 0.0
        self.x = x
        self.z = z
        self.density = density
        self.n = len(x) - 1

    def __str__(self):
        return f"n: {self.n} \nx:\n {self.x} \nz:\n {self.z} \nDensity:\n {self.density}"

To compare the numerical results, we provide here a simple function which computes the gravitational attraction of a horizontal cylinder of mass m that intersects the vertical plane at the coordinate given by the 2-D array axis. The point of observation is given by the coordinates in the 2-D array obs:

g = cylinder(obs, mass, axis)

def cylinder(obs, mass, axis):
    g = -2 * 6.674e-11 * mass / np.linalg.norm(obs - axis)**2 * (obs - axis)
    return g

The following code implements the Talwani approach:

ATOL = 1e-8
RTOL = 1e-8

def talwani(obs, polygon):
    gz = 0.0
    density = polygon.density
    x = polygon.x
    z = polygon.z
    n = polygon.n
    xp = obs[0]
    zp = obs[1]
    for k in range(n):   
        xv = x[k] - xp
        zv = z[k] - zp
        xvp1 = x[k+1] - xp
        zvp1 = z[k+1] - zp
        if not isclose(zv, zvp1, atol=ATOL, rtol=RTOL):
            a_i = xvp1 + zvp1 * (xvp1 - xv) / (zv - zvp1)
        theta_v = np.arctan2(zv, xv)
        theta_vp1 = np.arctan2(zvp1, xvp1)
        phi_i = np.arctan2(zvp1 - zv, xvp1 - xv)
        
        if theta_v < 0:
            theta_v += np.pi
        if theta_vp1 < 0:
            theta_vp1 += np.pi
            
        if isclose(xv, 0.0, atol=ATOL, rtol=RTOL): # A
            Zi = -a_i * np.sin(phi_i) * np.cos(phi_i) * (
                theta_vp1 - np.pi/2 + np.tan(phi_i) * np.log(np.cos(theta_vp1) * (np.tan(theta_vp1) - np.tan(phi_i))))

        elif isclose(xvp1, 0.0, atol=ATOL, rtol=RTOL): # B
            Zi = a_i * np.sin(phi_i) * np.cos(phi_i) * (
                theta_v - np.pi/2 + np.tan(phi_i) * np.log(np.cos(theta_v) * (np.tan(theta_v) - np.tan(phi_i))))

        elif isclose(zv, zvp1, atol=ATOL, rtol=RTOL): # C
            Zi = zv * (theta_vp1 - theta_v)

        elif isclose(xv, xvp1, atol=ATOL, rtol=RTOL): # D
            Zi = xv * np.log(np.cos(theta_v) / np.cos(theta_vp1))

        elif (
            np.isclose(theta_v, theta_vp1, atol=ATOL, rtol=RTOL)   # E
            or np.isclose(xv, zv, atol=ATOL, rtol=RTOL)           # F
            or np.isclose(xvp1, zvp1, atol=ATOL, rtol=RTOL)       # G
            or np.isclose(theta_vp1, phi_i, atol=ATOL, rtol=RTOL)):
            Zi = 0.0

        else:
            a_i = xvp1 + zvp1 * (xvp1 - xv) / (zv - zvp1)
            Zi = a_i * np.sin(phi_i) * np.cos(phi_i) * ( 
            theta_v - theta_vp1 + np.tan(phi_i) * 
            np.log( 
                np.cos(theta_v) * ( np.tan(theta_v) - np.tan(phi_i) ) /
                np.cos(theta_vp1) / ( np.tan(theta_vp1) - np.tan(phi_i) )))

        gz += Zi

    gz = gz * 2 * 6.674e-11 * density

    return gz

\[ Z_{i} =a_{i} \sin \phi_{i} \cos \phi_{i} \left[ \theta_{i}-\theta_{i+1} +\tan \phi_{i} \log _{e} \frac{\cos \theta_{i}\left(\tan \theta_{i}-\tan \phi_{i}\right)}{\cos \theta_{i+1}\left(\tan \theta_{i+1}-\tan \phi_{i}\right)}\right] \]

11.1 Numerical experiment

We define a cylinder of radius \(r=8\) m at a depth of \(z=20\) m that will be approximated in the \(x-z\)-plane by a sufficiently large number of line segments around its enclosing circle.

n = 32
xc = 0.0
zc = 20.0
radius = 8.0
x = np.array([xc + radius * np.cos(1e-6 + 2 * np.pi * k / n) for k in np.linspace(0, n-1, n)])
z = np.array([zc + radius * np.sin(1e-6 + 2 * np.pi * k / n) for k in np.linspace(0, n-1, n)])
poly = Polygon(x, z, 500)

fig, ax = plt.subplots()
plt.plot(poly.x, poly.z, marker='.')
plt.ylim((40, 0))
plt.xlim((-100, 100))
ax.set_aspect('equal', 'box');

We define a profile perpendicular to the strike of the cylinder.

The arrays gtal and gcyl store the values obtained by the Talwani method and the analytical solution for the cylinder, respectively.

profile = np.arange(-100, 101, 4)
npts = profile.size
gtal = np.zeros(npts)
gcyl = np.zeros(npts)

A = n * 0.5 * radius**2 * np.sin(2 * np.pi / n)

mass = poly.density * A # np.pi * radius**2
for k in range(npts):
    gtal[k] = talwani(np.array([profile[k], 0]), poly)
    gcyl[k] = cylinder(np.array([profile[k], 0]), mass, np.array([xc, zc]))[1]

The following plots illustrate the accuracy of the Talwani approach.

fig, (ax1, ax2) = plt.subplots(1, 2)
ax1.plot(profile, gtal, marker='.', label="Talwani")
ax1.plot(profile, gcyl, marker='.', label="cylinder")
ax1.legend()
ax2.plot(profile, gcyl - gtal, marker='.')
ax1.set_title('anomaly')
ax2.set_title('difference');

Comparison of analytical and numerical solutions (left) and absolute difference (right).

The relative error is

print(norm(gtal - gcyl) / norm(gcyl))
Ac = np.pi * radius**2
print(f"Area ratio is {A/Ac:.4e}") 

5.711301342967773e-15
Area ratio is 9.9359e-01

CautionNote

We have seen that in Python a comparison of floats to zero is dangerous. This is why we use np.isclose() with appropriate tolerances (ATOL and RTOL) throughout the talwani function instead of direct equality checks. Direct comparisons like xv == 0.0 can fail due to floating-point arithmetic precision issues, leading to incorrect results in the special case handling (cases A through G in the implementation).