Comparison between the magnetic field produced by a triaxial ellipsoid and an oblate ellipsoid

Import the required modules and functions



In [1]:

    
%matplotlib inline
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import BoundaryNorm
from matplotlib.ticker import MaxNLocator
from fatiando import gridder, utils
import triaxial_ellipsoid, oblate_ellipsoid
from mesher import TriaxialEllipsoid, OblateEllipsoid
from mesher import _coord_transf_matrix_triaxial, _coord_transf_matrix_oblate
import plot_functions as pf



In [2]:

    
# Set some plot parameters
from matplotlib import rcParams
rcParams['figure.dpi'] = 300.
rcParams['font.size'] = 6
rcParams['xtick.labelsize'] = 'medium'
rcParams['ytick.labelsize'] = 'medium'
rcParams['axes.labelsize'] = 'large'
rcParams['legend.fontsize'] = 'medium'
rcParams['savefig.dpi'] = 300.

Set some parameters for modelling



In [3]:

    
# The local-geomagnetic field
F, inc, dec = 60000, 50, 0

# Create a regular grid at z = 0 m
shape = (200, 200)
area = [-5000, 5000, -4000, 6000]
xp, yp, zp = gridder.regular(area, shape, z=0)

Triaxial ellipsoid versus prolate ellipsoid

This test compares the total-field anomalies produced by a triaxial ellipsoid with that produced by a sphere. The triaxial ellipsoid has semi-axes $a$, $b$, and $c$ equal to 700 m, 400.1 m, and 399.9 m, respectively, and the prolate ellipsoid has semi-axes $a$ and $b$ equal to 700 m and 400 m. Both bodies are centered at the point (0, 0, 1000).

Triaxial ellipsoid



In [5]:

    
triaxial = TriaxialEllipsoid(0, 0, 1000, 600.1, 599.9, 300, 50, 36, 0,
                             {'principal susceptibilities': [0.03, 0.02, 0.01],
                              'susceptibility angles': [-20, 20, 9],
                              'remanent magnetization': [0.7, -7, 10]})



In [6]:

    
magnetization_triaxial = triaxial_ellipsoid.magnetization(triaxial, F, inc, dec, demag=True)



In [7]:

    
magnetization_triaxial









    Out[7]:





array([ 1.63120758,  0.15733409,  0.36991097])



In [8]:

    
plt.close('all')
fig = plt.figure(figsize=(10,8))
ax = fig.gca(projection='3d')

# triaxial body
pf.draw_ellipsoid(ax, triaxial, body_color=(1,1,0), body_alpha=0.3)

ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_zlabel('z (m)')

ax.view_init(215, 20)

plt.tight_layout(True)

plt.show()

Oblate



In [20]:

    
oblate = OblateEllipsoid(triaxial.x, triaxial.y, triaxial.z, 
                         triaxial.small_axis, 0.5*(triaxial.large_axis + triaxial.intermediate_axis),
                         triaxial.strike, triaxial.dip, triaxial.rake,
                         {'principal susceptibilities': triaxial.props['principal susceptibilities'],
                          'susceptibility angles': triaxial.props['susceptibility angles'],
                          'remanent magnetization': triaxial.props['remanent magnetization']})



In [21]:

    
magnetization_oblate = oblate_ellipsoid.magnetization(oblate, F, inc, dec, demag=True)



In [22]:

    
magnetization_oblate









    Out[22]:





array([ 1.22359088, -0.13520746,  0.86163125])



In [23]:

    
plt.close('all')
fig = plt.figure(figsize=(10,8))
ax = fig.gca(projection='3d')

# oblate body
pf.draw_ellipsoid(ax, oblate, body_color=(0,1,1), body_alpha=0.3)

ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_zlabel('z (m)')

ax.view_init(215, 20)

plt.tight_layout(True)

plt.show()



In [24]:

    
plt.close('all')
fig = plt.figure(figsize=(10,8))
ax = fig.gca(projection='3d')

# triaxial body
pf.draw_ellipsoid(ax, triaxial, body_color=(1,1,0), body_alpha=0.3)

# oblate body
pf.draw_ellipsoid(ax, oblate, body_color=(0,1,1), body_alpha=0.3)

ax.set_xlabel('x (m)')
ax.set_ylabel('y (m)')
ax.set_zlabel('z (m)')

ax.view_init(215, 20)

plt.tight_layout(True)

plt.show()

Total-field anomalies



In [25]:

    
# total-field anomaly produced by the triaxial ellipsoid (in nT)
tf_t = triaxial_ellipsoid.tf(xp, yp, zp, [triaxial],
                             F, inc, dec)

# total-field anomaly produced by the oblate ellipsoid (in nT)
tf_p = oblate_ellipsoid.tf(xp, yp, zp, [oblate],
                           F, inc, dec)

# residuals
tf_r = tf_t - tf_p



In [26]:

    
plt.figure(figsize=(3.15, 7))
plt.axis('scaled')

ranges = np.max(np.abs([np.min(tf_t), np.max(tf_t),
                        np.min(tf_p), np.max(tf_p)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,1)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             tf_t.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
cbar = plt.colorbar()
plt.annotate(s='(a)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.subplot(3,1,2)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             tf_p.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(b)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

ranges = np.max(np.abs([np.min(tf_r), np.max(tf_r)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,3)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             tf_r.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlabel('y (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(c)', xy=(0.88,0.92), 
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.tight_layout()

plt.show()

Field components



In [27]:

    
# field components produced by the triaxial ellipsoid (in nT)
bx_t = triaxial_ellipsoid.bx(xp, yp, zp, [triaxial],
                             F, inc, dec)
by_t = triaxial_ellipsoid.by(xp, yp, zp, [triaxial],
                             F, inc, dec)
bz_t = triaxial_ellipsoid.bz(xp, yp, zp, [triaxial],
                             F, inc, dec)

bt = [bx_t, by_t, bz_t]

# field components produced by the oblate ellipsoid (in nT)
bx_p = oblate_ellipsoid.bx(xp, yp, zp, [oblate],
                           F, inc, dec)
by_p = oblate_ellipsoid.by(xp, yp, zp, [oblate],
                           F, inc, dec)
bz_p = oblate_ellipsoid.bz(xp, yp, zp, [oblate],
                           F, inc, dec)

bp = [bx_p, by_p, bz_p]

# residuals
bx_r = bx_t - bx_p
by_r = by_t - by_p
bz_r = bz_t - bz_p

br = [bx_r, by_r, bz_r]



In [28]:

    
plt.figure(figsize=(3.15, 7))
plt.axis('scaled')

ranges = np.max(np.abs([np.min(bx_t), np.max(bx_t),
                        np.min(bx_p), np.max(bx_p)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,1)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             bx_t.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
cbar = plt.colorbar()
plt.annotate(s='(a)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.subplot(3,1,2)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             bx_p.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(b)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

ranges = np.max(np.abs([np.min(bx_r), np.max(bx_r)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,3)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             bx_r.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlabel('y (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(c)', xy=(0.88,0.92), 
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.tight_layout()

plt.show()



In [29]:

    
plt.figure(figsize=(3.15, 7))
plt.axis('scaled')

ranges = np.max(np.abs([np.min(by_t), np.max(by_t),
                        np.min(by_p), np.max(by_p)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,1)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             by_t.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
cbar = plt.colorbar()
plt.annotate(s='(a)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.subplot(3,1,2)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             by_p.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(b)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

ranges = np.max(np.abs([np.min(by_r), np.max(by_r)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,3)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             by_r.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlabel('y (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(c)', xy=(0.88,0.92), 
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.tight_layout()

plt.show()



In [30]:

    
plt.figure(figsize=(3.15, 7))
plt.axis('scaled')

ranges = np.max(np.abs([np.min(bz_t), np.max(bz_t),
                        np.min(bz_p), np.max(bz_p)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,1)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             bz_t.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
cbar = plt.colorbar()
plt.annotate(s='(a)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.subplot(3,1,2)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             bz_p.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(b)', xy=(0.88,0.92),
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

ranges = np.max(np.abs([np.min(bz_r), np.max(bz_r)]))
levels = MaxNLocator(nbins=20).tick_values(-ranges, ranges)
cmap = plt.get_cmap('RdBu_r')
norm = BoundaryNorm(levels, ncolors=cmap.N, clip=True)

plt.subplot(3,1,3)
plt.contourf(0.001*yp.reshape(shape), 0.001*xp.reshape(shape),
             bz_r.reshape(shape), levels=levels,
             cmap = cmap, norm=norm)
plt.ylabel('x (km)')
plt.xlabel('y (km)')
plt.xlim(0.001*np.min(yp), 0.001*np.max(yp))
plt.ylim(0.001*np.min(xp), 0.001*np.max(xp))
plt.colorbar()
plt.annotate(s='(c)', xy=(0.88,0.92), 
             xycoords = 'axes fraction', color='k',
             fontsize = 10)

plt.tight_layout()

plt.show()



In [ ]: