Magnetic Models Comparison¶
In this tutorial, we'll compare various Magnetic Models.
Setup¶
Let's import the necessary libraries:
import math
import numpy as np
%matplotlib widget
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from ostk.physics.time import Scale
from ostk.physics.time import Instant
from ostk.physics.time import DateTime
from ostk.physics.environment.object.celestial import Earth
from ostk.physics.environment.magnetic import Earth as EarthMagneticModel
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[1], line 5
1 import math
3 import numpy as np
----> 5 get_ipython().run_line_magic('matplotlib', 'widget')
6 import matplotlib.pyplot as plt
7 from matplotlib.patches import Circle
File /usr/local/lib/python3.11/dist-packages/IPython/core/interactiveshell.py:2480, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
2478 kwargs['local_ns'] = self.get_local_scope(stack_depth)
2479 with self.builtin_trap:
-> 2480 result = fn(*args, **kwargs)
2482 # The code below prevents the output from being displayed
2483 # when using magics with decorator @output_can_be_silenced
2484 # when the last Python token in the expression is a ';'.
2485 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):
File /usr/local/lib/python3.11/dist-packages/IPython/core/magics/pylab.py:103, in PylabMagics.matplotlib(self, line)
98 print(
99 "Available matplotlib backends: %s"
100 % _list_matplotlib_backends_and_gui_loops()
101 )
102 else:
--> 103 gui, backend = self.shell.enable_matplotlib(args.gui)
104 self._show_matplotlib_backend(args.gui, backend)
File /usr/local/lib/python3.11/dist-packages/IPython/core/interactiveshell.py:3665, in InteractiveShell.enable_matplotlib(self, gui)
3662 import matplotlib_inline.backend_inline
3664 from IPython.core import pylabtools as pt
-> 3665 gui, backend = pt.find_gui_and_backend(gui, self.pylab_gui_select)
3667 if gui != None:
3668 # If we have our first gui selection, store it
3669 if self.pylab_gui_select is None:
File /usr/local/lib/python3.11/dist-packages/IPython/core/pylabtools.py:338, in find_gui_and_backend(gui, gui_select)
321 def find_gui_and_backend(gui=None, gui_select=None):
322 """Given a gui string return the gui and mpl backend.
323
324 Parameters
(...)
335 'WXAgg','Qt4Agg','module://matplotlib_inline.backend_inline','agg').
336 """
--> 338 import matplotlib
340 if _matplotlib_manages_backends():
341 backend_registry = matplotlib.backends.registry.backend_registry
ModuleNotFoundError: No module named 'matplotlib'
Comparison¶
Let's set up some magnetic models:
magnetic_models = {
EarthMagneticModel.Type.Dipole: None,
EarthMagneticModel.Type.EMM2010: None,
EarthMagneticModel.Type.EMM2015: None,
EarthMagneticModel.Type.EMM2017: None,
EarthMagneticModel.Type.IGRF11: None,
EarthMagneticModel.Type.IGRF12: None,
EarthMagneticModel.Type.WMM2010: None,
EarthMagneticModel.Type.WMM2015: None,
}
earth = Earth.default()
def init_magnetic_models(magnetic_models):
for key in magnetic_models:
if magnetic_models[key] is None:
magnetic_models[key] = EarthMagneticModel(key)
init_magnetic_models(magnetic_models)
instant = Instant.date_time(DateTime(2015, 1, 1, 0, 0, 0), Scale.UTC)
def calc_B_2d(magnetic_model, X, Z):
size_x = len(X)
size_z = len(X[0])
size = size_x * size_z
Bx = np.zeros(size)
Bz = np.zeros(size)
earth_radius_m = earth.get_equatorial_radius().in_meters()
i = 0
for x, z in np.vstack([X.ravel(), Z.ravel()]).T:
if math.sqrt(pow(abs(x), 2) + pow(abs(z), 2)) < earth_radius_m:
Bx[i] = 0.0
Bz[i] = 0.0
else:
try:
b = magnetic_model.get_field_value_at(np.array((x, 0.0, z)), instant).T
bx = b[0]
bz = b[2]
Bx[i] = bx
Bz[i] = bz
except Exception as e:
print(e)
print("x = {}".format(x))
print("z = {}".format(z))
raise e
i += 1
return Bx.reshape(size_x, size_z), Bz.reshape(size_x, size_z)
def calc_B_3d(magnetic_model, X, Y, Z):
size_x = len(X)
size_y = len(X[0])
size_z = len(X[0][0])
size = size_x * size_y * size_z
Bx = np.zeros(size)
By = np.zeros(size)
Bz = np.zeros(size)
earth_radius_m = earth.get_equatorial_radius().in_meters()
i = 0
for x, y, z in np.vstack([X.ravel(), Y.ravel(), Z.ravel()]).T:
if math.sqrt(pow(abs(x), 2) + pow(abs(y), 2) + pow(abs(z), 2)) < earth_radius_m:
Bx[i] = 0.0
By[i] = 0.0
Bz[i] = 0.0
else:
try:
b = magnetic_model.get_field_value_at(np.array((x, y, z)), instant).T
bx = b[0]
by = b[1]
bz = b[2]
Bx[i] = bx
By[i] = by
Bz[i] = bz
except Exception as e:
print(e)
print("x = {}".format(x))
print("y = {}".format(y))
print("z = {}".format(z))
raise e
i += 1
return (
Bx.reshape(size_x, size_y, size_z),
By.reshape(size_x, size_y, size_z),
Bz.reshape(size_x, size_y, size_z),
)
Display¶
2D plot:
lim = 100000e3
xlim = (-lim, +lim)
zlim = (-lim, +lim)
(nx, nz) = 64, 64
x = np.linspace(xlim[0], xlim[1], nx)
z = np.linspace(zlim[0], zlim[1], nz)
(X, Z) = np.meshgrid(x, z)
(Bx, Bz) = calc_B_2d(magnetic_models[EarthMagneticModel.Type.EMM2015], X, Z)
figure = plt.figure()
ax = figure.add_subplot(111)
# Plot Earth circle
ax.add_artist(Circle((0.0, 0.0), earth.get_equatorial_radius().in_meters(), color="b"))
# Plot B field
color = 2 * np.log(np.hypot(Bx, Bz))
ax.streamplot(
x,
z,
Bx,
Bz,
color=color,
linewidth=1,
cmap=plt.cm.inferno,
density=2,
arrowstyle="->",
arrowsize=1.5,
)
ax.set_xlabel("$x_{ITRF}\ [m]$")
ax.set_ylabel("$z_{ITRF}\ [m]$")
ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(zlim[0], zlim[1])
ax.set_aspect("equal")
plt.show()
/tmp/ipykernel_3907/1284092256.py:11: RuntimeWarning: divide by zero encountered in log
color = 2 * np.log(np.hypot(Bx, Bz))
3D plot:
lim = 10000e3
xlim = (-lim, +lim)
ylim = (-lim, +lim)
zlim = (-lim, +lim)
(nx, ny, nz) = (9, 9, 9)
x = np.linspace(xlim[0], xlim[1], nx)
y = np.linspace(ylim[0], ylim[1], ny)
z = np.linspace(zlim[0], zlim[1], nz)
(X, Y, Z) = np.meshgrid(x, y, z)
(Bx, By, Bz) = calc_B_3d(magnetic_models[EarthMagneticModel.Type.EMM2010], X, Y, Z)
figure = plt.figure()
ax = figure.add_subplot(projection="3d")
# Plot Earth sphere
u = np.linspace(0.0, 2.0 * np.pi, 100)
v = np.linspace(0.0, np.pi, 100)
earth_radius_m = float(earth.get_equatorial_radius().in_meters())
x = earth_radius_m * np.outer(np.cos(u), np.sin(v))
y = earth_radius_m * np.outer(np.sin(u), np.sin(v))
z = earth_radius_m * np.outer(np.ones(np.size(u)), np.cos(v))
ax.plot_surface(x, y, z, rstride=4, cstride=4, color="b", linewidth=0, alpha=1)
# Plot B field
ax.quiver(X, Y, Z, Bx, By, Bz, length=1e11, color="g")
ax.set_xlabel("$x_{ITRF}\ [m]$")
ax.set_ylabel("$y_{ITRF}\ [m]$")
ax.set_zlabel("$z_{ITRF}\ [m]$")
ax.set_xlim(xlim[0], xlim[1])
ax.set_ylim(ylim[0], ylim[1])
ax.set_zlim(zlim[0], zlim[1])
plt.show()