Verify equation 31

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.


    

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)


    

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 [ ]: