This notebook aims at reproducing the demagnetizing factors values presented by Stoner (1945).
In [1]:
from __future__ import division
%matplotlib inline
import numpy as np
import os
from matplotlib import pyplot as plt
from fatiando import utils
import mesher
import prolate_ellipsoid, oblate_ellipsoid, 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 [3]:
e2 = np.linspace(0.001, 0.99, 100)
In [4]:
# m = a/b
m = np.sqrt(1 - e2)
print m
In [5]:
# semi-axes (in m)
a0 = 100
b = a0*m
# demagnetizing factors
n11_prolate = []
n22_prolate = []
for bi in b:
ellipsoid = mesher.ProlateEllipsoid(0, 0, 0, a0, bi, 29, -50, 180)
N1, N2 = prolate_ellipsoid.demag_factors(ellipsoid)
n11_prolate.append(N1)
n22_prolate.append(N2)
In [6]:
# semi-axes (in m)
b0 = 100
a = b0*m
# demagnetizing factors
n11_oblate = []
n22_oblate = []
for ai in a:
ellipsoid = mesher.OblateEllipsoid(0, 0, 0, ai, b0, -13, 10, 180)
N1, N2 = oblate_ellipsoid.demag_factors(ellipsoid)
n11_oblate.append(N1)
n22_oblate.append(N2)
In [7]:
x = np.linspace(-1, 1, 11)
xlables = x.copy()
xlables[xlables < 0] *= -1
print x
print xlables
In [8]:
lenght = 3.27
plt.close('all')
plt.figure(figsize=(lenght, 0.6*lenght))
label_fontsize = 10
label_tex_fontsize = 10
plt.plot([0, 0], [0, 1], 'k-')
plt.plot([-1, 1], [0.5, 0.5], 'k-')
plt.plot([-1, 1], [1/3 + 2/15, 1/3 - 2/15], color=(0.6, 0.6, 0.6))
plt.plot([-1, 1], [1/3 - 1/15, 1/3 + 1/15], color=(0.6, 0.6, 0.6))
plt.plot(-e2, n11_oblate, color='r')
plt.plot(-e2, n22_oblate, color='g')
plt.plot(e2, n11_prolate, color='r')
plt.plot(e2, n22_prolate, color='g')
plt.ylabel('demagnetizing factor')
#plt.xlabel('$\epsilon^{2}$', fontsize=label_tex_fontsize)
plt.xlim(-1, 1)
plt.ylim(0, 1)
plt.xticks(x, xlables)
plt.yticks(np.linspace(0, 1, 6))
plt.annotate(s='Oblate', xy=(0.17,0.85),
xycoords = 'axes fraction', color='k',
fontsize=label_fontsize, backgroundcolor='w')
plt.annotate(s='Prolate', xy=(0.66,0.85),
xycoords = 'axes fraction', color='k',
fontsize=label_fontsize, backgroundcolor='w')
plt.annotate(s='$1 - m^{2}$', xy=(0.29,0.02), color='k',
xycoords='figure fraction', fontsize=label_fontsize)
plt.annotate(s='$(m^{2} - 1)/m^{2}$', xy=(0.64,0.02), color='k',
xycoords='figure fraction', fontsize=label_fontsize)
plt.grid()
plt.tight_layout()
plt.show()
Demagnetizing factors $\tilde{n}^{\dagger}_{11}$ (in red) and $\tilde{n}^{\dagger}_{22}$ (in green) for prolate (Eqs. 34 and 35) and oblate (Eqs. 36 and 37) ellipsoids as functions of eccentricity $\epsilon$. For prolate ellipsoids, $\epsilon^{2} = (m^{2} - 1)/m^{2}$, and for oblate ellipsoids, $\epsilon^{2} = 1 - m^{2}$, where $m = a/b$. The straight lines of slopes $-2/15$ and $+1/15$ (in grey) are tangential to the curves at $\epsilon^{2} = 0$.
In [9]:
fname = 'Stoner1945_TableI.txt'
file_path = os.path.join(os.path.pardir, 'data', fname)
m, Da = np.loadtxt(file_path, unpack=True)
In [10]:
# semi-axes (in m)
b0 = 100.1
c0 = 99.9
a = 100*m[m > 1]
# demagnetizing factors
n11_triaxial1 = []
n22_triaxial1 = []
n33_triaxial1 = []
for ai in a:
ellipsoid = mesher.TriaxialEllipsoid(0, 0, 0, ai, b0, c0, -43, 10, -1.7)
N1, N2, N3 = triaxial_ellipsoid.demag_factors(ellipsoid)
n11_triaxial1.append(N1)
n22_triaxial1.append(N2)
n33_triaxial1.append(N3)
In [11]:
# semi-axes (in m)
b0 = 100
a = b0*m[m > 1]
# demagnetizing factors
n11_prolate = []
n22_prolate = []
for ai in a:
ellipsoid = mesher.ProlateEllipsoid(0, 0, 0, ai, b0, 0, 40, 1)
N1, N2 = prolate_ellipsoid.demag_factors(ellipsoid)
n11_prolate.append(N1)
n22_prolate.append(N2)
In [12]:
# semi-axes (in m)
a0 = 100.1
b0 = 100
c = b0*m[m < 1]
# demagnetizing factors
n11_triaxial2 = []
n22_triaxial2 = []
n33_triaxial2 = []
for ci in c:
ellipsoid = mesher.TriaxialEllipsoid(0, 0, 0, a0, b0, ci, -43, 10, -1.7)
N1, N2, N3 = triaxial_ellipsoid.demag_factors(ellipsoid)
n11_triaxial2.append(N1)
n22_triaxial2.append(N2)
n33_triaxial2.append(N3)
In [13]:
# semi-axes (in m)
b0 = 100
a = b0*m[m < 1]
# demagnetizing factors
n11_oblate = []
n22_oblate = []
for ai in a:
ellipsoid = mesher.OblateEllipsoid(0, 0, 0, ai, b0, 15, 100, -7)
N1, N2 = oblate_ellipsoid.demag_factors(ellipsoid)
n11_oblate.append(N1)
n22_oblate.append(N2)
In [21]:
lenght = 3.27
plt.close('all')
plt.figure(figsize=(lenght, 4))
label_fontsize = 10
label_tex_fontsize = 10
plt.subplot(2, 1, 1)
plt.semilogy(m[m > 1], Da[m > 1], 'k-', linewidth=2)
plt.semilogy(m[m > 1], n11_prolate, 'r--', linewidth=1.5)
plt.semilogy(m[m > 1], n11_triaxial1, 'b-', linewidth=0.5)
plt.ylabel('$\\tilde{n}^{\dagger}_{11}$')
#plt.xlabel('$m$', fontsize=label_tex_fontsize)
plt.xlim(m[m > 1].min(), m[m > 1].max())
plt.grid()
plt.annotate(s='(a)', xy=(0.92,0.90),
xycoords = 'axes fraction', color='k',
fontsize=label_fontsize)
plt.subplot(2, 1, 2)
plt.plot(m[m < 1], Da[m < 1], 'k-', linewidth=2)
plt.plot(m[m < 1], n11_oblate, 'r--', linewidth=1.5)
plt.plot(m[m < 1], n33_triaxial2, 'b-', linewidth=0.5)
plt.ylabel('$\\tilde{n}^{\dagger}_{11}$')
plt.xlabel('$m$', fontsize=label_tex_fontsize)
plt.xlim(m[m < 1].min(), m[m < 1].max())
plt.grid()
plt.annotate(s='(b)', xy=(0.92,0.90),
xycoords = 'axes fraction', color='k',
fontsize=label_fontsize)
plt.tight_layout()
plt.show()
Comparison between the demagnetizing factors calculated here by using our routines and those presented by Stoner (1945, Table I) (continuous black lines). (a) The blue and red lines represent, respectively, the demagnetizing factor $\tilde{n}^{\dagger}_{11}$ of triaxial ellipsoids (Eq. 29) and the demagnetizing factor $\tilde{n}^{\dagger}_{11}$ of prolate ellipsoids (Eq. 34). (b) The blue and red lines represent, respectively, the demagnetizing factor $\tilde{n}^{\dagger}_{33}$ of triaxial ellipsoids (Eq. 31) and the demagnetizing factor $\tilde{n}^{\dagger}_{11}$ of oblate ellipsoids (Eq. 36).
In [ ]: