In [1]:

    

"""
GravMag: 3D forward modeling of total-field magnetic anomaly using triaxial
ellipsoids (model with induced and remanent magnetization)
"""
import numpy as np
from matplotlib import pyplot as plt
import matplotlib
%matplotlib inline
from fatiando import mesher, gridder, utils
from fatiando.gravmag import ellipsoid_triaxial,ellipsoid_triaxial2,ellipsoid_prolate,ellipsoid_oblate


    

In [6]:

    

# The regional field
inten,inc, dec = 60000., 30, -15
F_unit = utils.ang2vec(1, inc, dec)

bounds = [-5000, 5000, -5000, 5000, 0, 5000]
# Create a regular grid at 100m height
shape = (200, 200)
area = bounds[:4]
xp, yp, zp = gridder.regular(area, shape, z=0)


    

In [3]:

    

I_l = []
D_l = []
k_l = np.linspace(0.0001, 0.5, 100)

for k in k_l:
    model = [mesher.EllipsoidTriaxial(0.,0.,500.,500.,100.,50.,0.,0.,0., 
                                     {'remanence': [0., 0., 0.],
                                      'k': [k, k, k, 0., 90., 90.]})]

    # Calculate the anomaly for a given regional field
    tf,N1,N2,N3,JRD_ang = ellipsoid_triaxial.tf_c(xp,yp,zp,inten,inc,dec,model)
    I_l.append(JRD_ang[1])
    D_l.append(JRD_ang[2])

cte_inc = np.zeros_like(I_l) + inc
cte_dec = np.zeros_like(I_l) + dec
k_l = k_l


    

In [4]:

    

(D_l[15]+D_l[17])/2.+15


    

    
Out[4]:





array([ 0.68445636])

In [5]:

    

chi_max = 0.05/N3
chi_max


    

    
Out[5]:





0.077088922153742925

In [4]:

    

model = [mesher.EllipsoidTriaxial(0.,0.,500.,500.,100.,50.,0.,0.,0., 
                                     {'remanence': [0., 0., 0.],
                                      'k': [0.077088922153742925, 0.077088922153742925, 0.077088922153742925, 0., 90., 90.]})]

tf,N1,N2,N3,JRD_ang = ellipsoid_triaxial.tf_c(xp,yp,zp,inten,inc,dec,model)


    

In [14]:

    

JRD_ang


    

    
Out[14]:





[3.6453945132060279, array([ 29.55497116]), array([-14.34616156])]

In [15]:

    

mag_v = utils.ang2vec(3.6453945132060279, 29.55497116, -14.34616156)


    

In [18]:

    

mag_a = 0.077088922153742925*(inten/(4*np.pi*100))*F_unit


    

In [19]:

    

# relative error
mag_max_norm = np.linalg.norm(mag_v, ord = 2)
mag_max_approx_norm = np.linalg.norm(mag_a, ord = 2)
delta_mag_max_norm = np.linalg.norm(mag_v - mag_a, ord = 2)

epsilon_calculated = delta_mag_max_norm/mag_max_norm

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


    

    




epsilon = 1.593 percent

In [15]:

    

plt.figure(figsize=(15,10))
axes = plt.gca()
plt.plot([0.1,0.1], [-20.,35.], 'k-', linewidth=3.0)
plt.plot([0.077088922153742925,0.077088922153742925], [-20.,35.], 'y-', linewidth=3.0)
plt.plot(k_l,I_l, label='$I$', linewidth=3.0)
plt.plot(k_l,D_l, label='$D$', linewidth=3.0)
plt.plot(k_l,cte_inc, '--k', label='$I_{0}$', linewidth=3.0 )
plt.plot(k_l,cte_dec, '--r', label='$D_{0}$', linewidth=3.0)
plt.legend(loc='best', ncol=4, fontsize=28)
plt.ylabel('($^{\circ}$)', fontsize=30)
plt.xlabel('$\chi$', fontsize=30)
axes.set_ylim(-20.,35.)
#axes.set_xlim([0,5])
matplotlib.rcParams.update({'font.size': 18})
plt.tight_layout()
#plt.savefig('..\\manuscript\\figures\\test_k_triaxial.pdf', dpi = 600, facecolor='w', bbox_inches='tight')


    

In [ ]: