Routines to easily access information about magnetic fields of planets in our solar system and visualize them in both 2D and 3D. These require NumPy, Matplotlib and SciPy (pronounced "Sigh Pie"). Other than that, the following external libraries are used for a few different functions:
- 2D plotting for map projections other than Hammer : Cartopy library ( see more under Projections )
- Potential extrapolation: SHTns library
- Writing vts files for 3D visualisation: PyEVTK library
This gives access to all the relevant properties of a planet and has methods to plot
the field and write a vts file for 3D visualization. Usage:
from planetmagfields import *
p = planet(name='earth',datDir='planetmagfields/data/')This displays the some information about the planet
Planet: Earth
l_max = 13
Dipole tilt (degrees) = -9.410531
and gives access to variables associated with the planet such as:
p.lmax: maximum spherical harmonic degree till which data is availablep.glm,p.hlm: the Gauss coefficientsp.Br: computed radial magnetic field at surfacep.dipTheta: dipole tilt with respect to the rotation axisp.dipPhi: dipole longitude ( in case zero longitude is known, applicable to Earth )p.idx: indices to get values of Gauss coefficients
Example:
In [1]: from planetmagfields import *
In [2]: p = planet(name='jupiter')
Planet: Jupiter
l_max = 10
Dipole tilt (degrees) = 10.307870
In [3]: p.glm[p.idx[2,0]] # g20
Out[3]: 11670.4
In [4]: p.hlm[p.idx[4,2]] # h42
Out[4]: 27811.2as well as the functions:
This function plots a 2D surface plot of the radial magnetic field at radius r given in terms of the surface radius.
For example,
from planetmagfields import *
p = planet(name='earth',datDir='planetmagfields/data/')
p.plot(r=1,proj='Mollweide')produces the info mentioned above first and then the following plot of Earth's magnetic field using a Mollweide projection
while
from planetmagfields import *
p = planet(name='jupiter',datDir='planetmagfields/data/')
p.plot(r=0.85,proj='Mollweide')produces the following info about Jupiter and then plot that follows
Planet: Jupiter
l_max = 10
Dipole tilt (degrees) = 10.307870
This can be compared with Fig. 1 g from Moore et al. 2018
This function computes the Lowes spectrum of a planet at a given radius. It adds an array p.emag_spec which contains the energy at different spherical harmonic degrees and two variables p.dipolarity and p.dip_tot which provide the fraction of energies in the axial dipole and the total dipole with respect to the total energy at all degrees. Usage example:
from planetmagfields import *
p = planet(name='jupiter')
p.spec()will provide variables
In [8]: p.dipolarity
Out[8]: 0.7472047129875864
In [9]: p.dip_tot
Out[9]: 0.7719205112704368
In [10]: p.emag_spec
Out[10]:
array([0.00000000e+00, 3.47735422e+11, 2.36340423e+10, 2.12851278e+10,
1.75661779e+10, 1.92219842e+10, 9.91200756e+09, 3.34535475e+09,
3.95317971e+09, 2.59333412e+09, 1.23423769e+09])and will produce Jupiter's surface spectrum:
The plotting can be suppressed setting the logical p.spec(iplot=False).
This function writes a vts file that can be used to produce 3D visualizations of field lines with Paraview/VisIt. Usage:
p.writeVtsFile(potExtra=True, ratio_out=2, nrout=32)where,
potExtra: bool, whether to use potential extrapolationratio_out: float, radius till which the field would be extrapolated in terms of the surface radiusnrout: radial resolution for extrapolation
Example of a 3D image produced using Paraview for Neptune's field, extrapolated till 5 times the surface radius.
The planet class also provides a function for producing a filtered view of the radial magnetic field using the function plot_filt.
This function can take in either an arbitrary array of spherical harmonic degrees and orders or cut-off values. This is illustrated
below with examples, assuming the user is in the repository directory.
from planetmagfields import *
p = planet(name='saturn')
p.plot_filt(r=0.75,lCutMin=4,proj='Mollweide')
Compare this with Fig. 20 B from Cao et al. 2020
from planetmagfields import *
p = planet(name='jupiter')
p.plot_filt(r=1,larr=[1,2,3],marr=[3],proj='Mollweide')
from planetmagfields import *
p = planet(name='earth')
p.plot_filt(r=1,lCutMin=5,mmin=4,proj='Mollweide')
The potextra module provides a method for potential extrapolation of a planet's magnetic field.
This uses the SHTns library for spherical harmonic transforms.
Usage example:
import numpy as np
from planetmagfields import *
p = planet('saturn')
ratio_out = 5 # Ratio (> 1) in terms of surface radius to which to extrapolate
nrout = 32 # Number of grid points in radius between 1 and ratio_out
rout = np.linspace(1,ratio_out,nrout)
brout, btout, bpout = potextra.extrapot(p.lmax,1.,p.Br,rout)$ ./magField.py <planet> <radius> <projection>
This will plot the radial magnetic field of a planet (any of the names from the list below, case insensitive) at a radius given in terms of the surface radius. The default is the surface field. For example,
$ ./magField.py earth Mollweide
displays the same information as above about Earth's field and produces the surface field of Earth while
$ ./magField.py jupiter 0.85 Mollweide
produces the same plot of Jupiter's field as shown before.
$ ./magField.py all <radius> <projection>
would produce a table of information about dipole tilt for each planet and magnetic field maps of all different planets at the given radius in a single figure.
For example:
$ ./magField all 0.9 Mollweide
would give
|=========|======|=======|
|Planet | Theta| Phi |
|=========|======|=======|
|Mercury | 0.0 | 0.0 |
|Earth | -9.4 | -72.7 |
|Jupiter | 10.3 | -16.6 |
|Saturn | 0.0 | 0.0 |
|Uranus | 58.6 | -53.6 |
|Neptune | 46.9 | -72.0 |
|Ganymede | -4.2 | 25.5 |
|---------|------|-------|
followed by the following plot
All the Gauss coefficients in the collected data are Schmidt semi-normalized. Only the data for Earth uses a Condon-Shortley phase, the others do not.
By default, the plot function tries to use the Mollweide projection. However, using the power of the cartopy library, any projection from this list is supported. In the absence of the cartopy library, the 2D plots fall back to the internally written Hammer projection. Examples of Jupiter's radial magnetic field at r=0.85 with different projections are shown below:
In [1]: from planetmagfields import *
In [2]: p = planet(name='jupiter')
Planet: Jupiter
l_max = 10
Dipole tilt (degrees) = 10.307870
In [3]: p.plot(r=0.85)
In [4]: projlist=['Mercator','Robinson','Stereographic','AzimuthalEquidistant']
In [5]: for k,proj in enumerate(projlist):
...: p.plot(r=0.85,proj=proj)
...: savefig('jup_r0_85'+proj+'.png',dpi=150,bbox_inches='tight')
...: close()In[3] produces the plot of Jupiter's field already shown above. In[5] produces the following figures with the Mercator, Robinson, Stereographic and azimuthal equidistant projections, respectively.
This also works with the magField.py script for quick plotting. Examples:
./magField.py earth 0.9 Robinsonor even with plots of all planets together
./magField.py all 0.9 RobinsonNote that I have chosen to keep the projection information out of the plot titles to prevent too much text.
| ❗ | The Orthographic projection often does not create correct plots, be cautious while using it |
|---|
Mercury : Anderson et al. 2012
Earth : IGRF 13, Alken et al. 2021
Jupiter : JRM09, Connerney et al. 2018
Saturn : Cassini 11+, Cao et al. 2020
Uranus : Connerny et al. 1987
Neptune : Connerny et al. 1991
Ganymede: Kivelson et al. 2002
I would like to thank Regupathi Angappan for motivating me to create this package. I thank Jon Aurnou for testing it out and pointing out runtime errors. Lastly, I would like to thank Thomas Gastine for comparing the plots with real data and pointing out a normalization error which has been fixed.