Verify equation 31

Import the required modules


In [1]:
%matplotlib inline
import numpy as np
from numpy.random import rand
from matplotlib import pyplot as plt
from matplotlib.colors import Normalize
from matplotlib.colorbar import ColorbarBase
from matplotlib.cm import get_cmap
from fatiando import mesher, utils
from fatiando.gravmag import triaxial_ellipsoid

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.

Randomly generated models


In [16]:
# isostropic susceptibility (in SI)
k = 0.077088922153742925

# Ellipsoid semi-axes (in m)
a = 500
b = 100
c = 50

# demagnetizing factors
n11, n22, n33 = triaxial_ellipsoid.demag_factors(a, b, c)

# orietation angles (in degrees)
strike = 180.
dip = 0.
rake = 0.

# auxiliary orientation angles (in radians)
alpha, gamma, delta = triaxial_ellipsoid.structural_angles(strike, dip, rake)

V = triaxial_ellipsoid.V(alpha, gamma, delta)

In [28]:
n33


Out[28]:
0.64860162268557975

In [17]:
F = 60000 # intensity (nT) of the local-geomagnetic field
inc = 30 # inclination (degrees)
dec = -15 # declination (degrees)

In [18]:
# Local-geomagnetic field
F_unit = utils.ang2vec(1, inc, dec)

# isotropic susceptibility tensor
K = k*np.identity(3)

Magnetization


In [19]:
# in the local coordinate system
mag = triaxial_ellipsoid.magnetization(K, F, inc, dec, 0., 0., 0., V, [a,b,c])

# in the main coordinate system
mag = np.dot(V, mag)

In [20]:
aux = np.array([1. - k*n11,
                1. - k*n22,
                1. - k*n33])
Lambda = np.dot(V, np.dot(np.diag(1./aux), V.T))
mag2 = np.dot(Lambda, k*(F/(4*np.pi*100))*F_unit)

In [21]:
np.allclose(mag, mag2)


Out[21]:
False

In [29]:
utils.vec2ang(mag)


Out[29]:
[3.6266306536030735, 28.900711736817904, -14.692171478310154]

In [30]:
utils.vec2ang(mag2)


Out[30]:
[3.7412390716954858, 31.184974996392089, -15.322363506454506]

In [32]:
mag3 = k*(F/(4*np.pi*100))*F_unit

In [35]:
utils.vec2ang(mag3)


Out[35]:
[3.6807249055183511, 30.000000000000004, -14.999999999999998]

In [38]:
# relative error
mag_max_norm = np.linalg.norm(mag2, ord = 2)
mag_max_approx_norm = np.linalg.norm(mag3, ord = 2)
delta_mag_max_norm = np.linalg.norm(mag2 - mag3, ord = 2)

epsilon_calculated = delta_mag_max_norm/mag_max_norm

print 'epsilon = %.3f percent' % (epsilon_calculated*100)


epsilon = 2.656 percent

In [ ]: