# Tracing magnetic field lines¶

# Import notebook dependencies
import os
import sys
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

from src import sha, IGRF13_FILE


## Tracing magnetic field lines¶

The basic principle in following a field line at a point is to calculate the field vector there, then ‘take a step’ in the direction of the field to a new point, and repeat. The smaller the step, the higher the accuracy. This can be expressed as

\begin{split} \begin{align} \delta \mathbf{r} &= k \mathbf{B}\\ &=\left(\delta r, \thinspace r \delta \theta, \thinspace r\sin{\theta}\delta \phi\right) = k\left(B_r, \thinspace B_\theta, \thinspace B_\phi\right) \end{align} \end{split}

with $$k$$ a constant to scale the step size. We will use this to investigate the path of field lines computed using the IGRF, but first explore the simpler case of axisymmetric fields, where an explicit formula can be derived, avoiding the need for the ‘stepping strategy’.

### Axisymmetric multipole field lines¶

In this case there is no variation with longitude ($$\phi$$), so that

$\delta \mathbf{r}=\left(\delta r, r \delta \theta\right) = k(B_r, B_\theta)$

using spherical coordinates, which gives the equation of the field line as

$\frac{1}{r}\frac{dr}{d_\theta}= \frac{B_r}{B_\theta}.$

This equation can be solved when the field is expressed in terms of spherical harmonics for the axisymmetric terms with $$m=0$$. In a spherical harmonic expansion, the zonal terms $$P_n^o$$ correspond to axial multipoles, ($$P_1^0$$ is a dipole, $$P_2^0$$ is a quadrupole, $$P_3^0$$ is an octupole and $$P_4^0$$ is an hexadecapole.) The field line equation for $$P_n^0$$ is (Willis and Young, 1987; Jeffreys, 1988)

$r^n=\sin{\theta} \thinspace P_n^1(\cos{\theta}).$

Note the power of spherical harmonic degree $$n$$, and that the Associated Legendre polynomial in the equation has $$m=1$$ not $$m=0$$.

For the field line that passes through a given starting position ($$r_0$$, $$\theta_0$$),

$r_0^n = k P_n^1(\cos(\theta_0)) \sin(\theta_0),$

which gives the following expression for the scaling factor,

$k = \frac{r_0^n}{P_n^1(\cos(\theta_0)) \sin(\theta_0)}$

is derived.

The (un-normalised) forms of $$P_n^1$$ for $$n=1, \dots 4$$ are

\begin{split} \begin{align} P_1^1&=\sin^2{\theta}\\ P_2^1&=\cos{\theta}\sin^2{\theta}\\ P_3^1&=(5\cos^2{\theta}-1)\sin^2{\theta}\\ P_4^1&=(7\cos^3{\theta}-3\cos{\theta})\sin^2{\theta} \end{align} \end{split}

These, along with the expression for $$k$$, are used in the functions below to compute values of $$r^n$$ for $$n=1, \dots 4$$.

def dipole(r0, theta_0, th):      # P(n=1, m=0)
k = r0/(np.sin(theta_0)**2)
return (k*np.sin(th)**2)

def quadpole(r0, theta_0, th):    # P(n=2, m=0)
k = r0**2/(np.sin(theta_0)**2*np.cos(theta_0))
P = np.cos(th)*np.sin(th)**2
return (np.sqrt(np.abs(k)*np.abs(P)))

def octpole(r0, theta_0, th):     # P(n=3, m=0)
k = r0**3/(np.sin(theta_0)**2*(5*np.cos(theta_0)**2-1))
P = (5*np.cos(th)**2-1)*np.sin(th)**2
return ((np.abs(k)*np.abs(P))**(1/3))

def hexdpole(r0, theta_0, th):    # P(n=4, m=0)
k = r0**4/((7*np.cos(theta_0)**3-3*np.cos(theta_0))*np.sin(theta_0)**2)
P =(7*np.cos(th)**3-3*np.cos(th))*np.sin(th)**2
return ((np.abs(k)*np.abs(P))**(1/4))


You can plot dipole, quadrupole, octupole and hexadecapole field lines using the code below. Start by choosing the type of multipole field lines to plot by setting the value of pole_type to the desired spherical harmonic degree $$n$$, where $$n=1, \dots 4$$, and also setting some values of starting colatitude $$\theta_0$$ in the list theta. The program will draw the field lines passing through the point on the Earth’s surface at each value of colatitude (represented as radius=1 so the axes are scaled in Earth radii). Different choices of colatitudes may be more illuminating in different cases.

### >> USER INPUT HERE: Set the plotting parameters below (pole type and starting colatitudes)¶

pole_type options:
1 = dipole
3 = octupole

# Pole type
pole_type = 2

# Experiment by changing a list of starting colatitudes
theta   = [5, 10, 15, 20, 30, 40]


Now create the plot.

d2r     = np.deg2rad
axpoles = {1:dipole, 2:quadpole, 3:octpole, 4:hexdpole}
clines  = ['red', 'blue', 'grey', 'purple', 'brown', 'purple', 'pink', 'orange', 'magenta', 'olive', 'cyan']

# Define the "Earth"
r0  = 1
the = d2r(np.linspace(0, 360, 1000))
xe  = r0*np.sin(the)
ye  = r0*np.cos(the)

fig, ax = plt.subplots(figsize=(15, 15))
ax.set_aspect('equal')
ax.fill(xe, ye, color='lightgrey') # Plot the Earth

ic = -1
for i in theta:
ic += 1
theta_0 = d2r(i)
th = np.linspace(theta_0, d2r(90), 1000) 