{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "

IAGA Summer School 2019

\n", "\n", "

Tracing magnetic field lines

" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Import notebook dependencies\n", "\n", "import os\n", "import sys\n", "import numpy as np\n", "import pandas as pd\n", "import matplotlib.pyplot as plt\n", "%matplotlib inline\n", "\n", "sys.path.append('..')\n", "from src import sha_lib as sha\n", "\n", "# Set path to IGRF coefficients file\n", "IGRF12_FILE = os.path.abspath('../data/external/igrf12coeffs.txt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. Tracing magnetic field lines\n", "\n", "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\n", "\n", "$$\n", "\\begin{align}\n", "\\delta \\mathbf{r} &= k \\mathbf{B}\\\\\n", "&=\\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)\n", "\\end{align}\n", "$$\n", "\n", "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'.\n", "\n", "## Axisymmetric multipole field lines\n", "\n", "In this case there is no variation with longitude ($\\phi$), so that\n", "\n", "$$\\delta \\mathbf{r}=\\left(\\delta r, r \\delta \\theta\\right) = k(B_r, B_\\theta)$$\n", "\n", "using spherical coordinates, which gives the equation of the field line as\n", "\n", "$$\\frac{1}{r}\\frac{dr}{d_\\theta}= \\frac{B_r}{B_\\theta}.$$\n", "\n", "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)\n", "\n", "$$r^n=\\sin{\\theta} \\thinspace P_n^1(\\cos{\\theta}).$$\n", "\n", "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$),\n", "\n", "$$r_0^n = k P_n^1(\\cos(\\theta_0)) \\sin(\\theta_0),$$\n", "\n", "from which an expression for the scaling factor in the first equation\n", "\n", "$$k = \\frac{r_0^n}{P_n^1(\\cos(\\theta_0)) \\sin(\\theta_0)}$$\n", "\n", "is derived.\n", "\n", "The (un-normalised) forms of $P_n^1$ for $n=1, \\dots 4$ are\n", "\n", "$$\n", "\\begin{align}\n", "P_1^1&=\\sin^2{\\theta}\\\\\n", "P_2^1&=\\cos{\\theta}\\sin^2{\\theta}\\\\\n", "P_3^1&=(5\\cos^2{\\theta}-1)\\sin^2{\\theta}\\\\\n", "P_4^1&=(7\\cos^3{\\theta}-3\\cos{\\theta})\\sin^2{\\theta}\n", "\\end{align}\n", "$$\n", "\n", "These, along with the expression for $k$, are used in the functions below to compute values of $r^n$ for $n=1, \\dots 4$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def dipole(r0, theta_0, th): # P(n=1, m=0)\n", " k = r0/(np.sin(theta_0)**2)\n", " return (k*np.sin(th)**2)\n", " \n", "def quadpole(r0, theta_0, th): # P(n=2, m=0)\n", " k = r0**2/(np.sin(theta_0)**2*np.cos(theta_0))\n", " P = np.cos(th)*np.sin(th)**2\n", " return (np.sqrt(np.abs(k)*np.abs(P)))\n", "\n", "def octpole(r0, theta_0, th): # P(n=3, m=0)\n", " k = r0**3/(np.sin(theta_0)**2*(5*np.cos(theta_0)**2-1))\n", " P = (5*np.cos(th)**2-1)*np.sin(th)**2\n", " return ((np.abs(k)*np.abs(P))**(1/3))\n", " \n", "def hexdpole(r0, theta_0, th): # P(n=4, m=0)\n", " k = r0**4/((7*np.cos(theta_0)**3-3*np.cos(theta_0))*np.sin(theta_0)**2)\n", " P =(7*np.cos(th)**3-3*np.cos(th))*np.sin(th)**2\n", " return ((np.abs(k)*np.abs(P))**(1/4))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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.\n", "### >> USER INPUT HERE: Set the plotting parameters below (pole type and starting colatitudes)\n", "**pole_type** options:
\n", "1 = dipole
\n", "2 = quadrupole
\n", "3 = octupole
\n", "4 = hexadecapole
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Pole type\n", "pole_type = 1\n", "\n", "# Experiment by changing a list of starting colatitudes\n", "theta = [10, 15, 20, 30, 40]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now create the plot." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "d2r = np.deg2rad\n", "axpoles = {1:dipole, 2:quadpole, 3:octpole, 4:hexdpole}\n", "names = {1:'Dipole', 2:'Quadrupole', 3:'Octupole', 4:'Hexadecapole'}\n", "clines = ['red', 'blue', 'grey', 'purple', 'brown', 'purple', 'pink', 'orange', 'magenta', 'olive', 'cyan']\n", "\n", "# Define the \"Earth\"\n", "r0 = 1\n", "the = d2r(np.linspace(0, 360, 1000))\n", "xe = r0*np.sin(the)\n", "ye = r0*np.cos(the)\n", "\n", "fig, ax = plt.subplots(figsize=(15, 15))\n", "ax.set_aspect('equal')\n", "ax.fill(xe, ye, color='lightgrey') # Plot the Earth\n", "\n", "ic = -1\n", "for i in theta:\n", " ic += 1\n", " theta_0 = d2r(i)\n", " th = np.linspace(theta_0, d2r(90), 1000)\n", " rad = axpoles[pole_type](r0, theta_0, th)\n", " xb = rad*np.sin(th)\n", " yb = rad*np.cos(th)\n", " xb[np.where(rad<1)] = np.nan\n", " yb[np.where(rad<1)] = np.nan\n", " ax.plot(xb, yb, color = clines[ic%10])\n", " # Assume a symmetrical distribution in the four quadrants\n", " ax.plot( xb, -yb, color = clines[ic%10])\n", " ax.plot(-xb, yb, color = clines[ic%10])\n", " ax.plot(-xb, -yb, color = clines[ic%10])\n", "ax.set_xlabel('Earth radii', fontsize=16)\n", "ax.set_ylabel('Earth radii', fontsize=16)\n", "ax.set_title(names[pole_type]+' field lines', fontsize=24);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", "Jeffreys, B. (1988) ‘Derivations of the equation for the field lines of an axisymmetric multipole’, Geophysical Journal International. John Wiley & Sons, Ltd (10.1111), 92(2), pp. 355–356. doi: 10.1111/j.1365-246X.1988.tb01148.x.\n", "\n", "Willis, D. M. and Young, L. R. (1987) ‘Equation for the field lines of an axisymmetric magnetic multipole’, Geophysical Journal International. Oxford University Press, 89(3), pp. 1011–1022. doi: 10.1111/j.1365-246X.1987.tb05206.x." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2. Using the IGRF to trace field lines (and find conjugate points)\n", "\n", "The IGRF gives a fuller 3-D representation of the geomagnetic field, rather than the axisymmetric 2-D representation used above, and so the 'stepping strategy' is needed to follow field lines. The dipole field dominates, as we would expect, but it is interesting to to see how the colatitude and longitude change along the path. The starting and ending points at the Earth's surface are 'connected' by a field line - they are conjugate points. (The IGRF for 2015.0 is used in the code below.)\n", "### >> >> USER INPUT HERE: Set the input parameters (listed below)\n", " 1. Starting colatitude (degrees)\n", " 2. Starting longitude (degrees)\n", " 3. Step size (km)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Starting colatitude and longitude (east) in degrees\n", "theta_0 = 30\n", "phi_0 = 330\n", "\n", "# Step size for the field line tracing in km\n", "step_size = 10" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now do the calculation and plot the results. *** This may take a few seconds - be patient! ***" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "d2r = np.deg2rad\n", "r2d = np.rad2deg\n", "fcalc = lambda x: np.sqrt(np.dot(x,x))\n", "\n", "igrf12 = pd.read_csv(IGRF12_FILE, delim_whitespace=True, header=3)\n", "gh2015 = np.append(0., igrf12['2015.0'])\n", "NMAX = 13\n", "\n", "# Initialise variables at starting point\n", "r0 = 6371.2\n", "thrd = d2r(theta_0)\n", "phrd = d2r(phi_0)\n", "track = [(r0, thrd, phrd)] # Store coordinates of points on the field line\n", "bxyz = sha.shm_calculator(gh2015, NMAX, r0, theta_0, phi_0, 'Geocentric')\n", "eff = fcalc(bxyz)\n", "lamb = step_size/eff\n", "rad = r0\n", "\n", "# Allow a maximum number of steps in the iteration for the field line to\n", "# return to the Earth's surface\n", "maxstep = 10000\n", "newrad = r0+0.001\n", "\n", "step = 0\n", "while step <= maxstep and newrad >= r0:\n", " rad, th, ph = track[step]\n", " lx, ly, lz = tuple(el*lamb for el in bxyz)\n", " newrad = rad+lz\n", " newth = th+lx/rad\n", " newph = ph-ly/(rad*np.sin(th))\n", " track += [(newrad, newth, newph)]\n", " bxyz = sha.shm_calculator(gh2015, NMAX, newrad, r2d(newth), \\\n", " r2d(newph),'Geocentric')\n", " lamb = step_size/fcalc(bxyz)\n", " step += 1\n", "\n", "rads = np.array([r[0] for r in track])\n", "ths = np.array([t[1] for t in track])\n", "phs = np.array([p[2] for p in track])\n", "\n", "fig, (ax0, ax1, ax2) = plt.subplots(nrows=3, ncols=1, figsize=(10, 12))\n", "ax0.plot(rads, r2d(ths))\n", "ax0.invert_yaxis()\n", "ax0.set_xlabel('Radial distance (km)', fontsize=14)\n", "ax0.set_ylabel('Colatitude (degrees)', fontsize=14)\n", "ax0.set_title('Field line tracing using the IGRF', fontsize=20)\n", "ax1.plot(rads, r2d(phs), color='red')\n", "ax1.set_xlabel('Radial distance (km)', fontsize=14)\n", "ax1.set_ylabel('Longitude (degrees)', fontsize=14)\n", "ax2.plot(r2d(phs), r2d(ths), color='green')\n", "ax2.set_xlabel('Longitude (degrees)', fontsize=14)\n", "ax2.set_ylabel('Colatitude (degrees)', fontsize=14)\n", "ax2.invert_yaxis()\n", "fig.tight_layout()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "print('\\nCoordinates of the starting point:')\n", "print('\\tColatitude :', '{0:.1f}'.format(theta_0), 'degrees')\n", "print('\\tLongitude :', '{0:.1f}'.format(phi_0), 'degrees')\n", "print('\\nCoordinates of the conjugate point:')\n", "print('\\tColatitude :', '{0:.1f}'.format(r2d(newth)), 'degrees')\n", "print('\\tLongitude :', '{0:.1f}'.format(r2d(newph)%360), 'degrees')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 2 }