{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Spherical harmonic models 1" ] }, { "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", "from src import sha, mag, IGRF13_FILE" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Spherical harmonics and representing the geomagnetic field" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The north (X), east (Y) and vertical (Z) (downwards) components of the internally-generated geomagnetic field at colatitude $\\theta$, longitude $\\phi$ and radial distance $r$ (in geocentric coordinates with reference radius $a=6371.2$ km for the IGRF) are written as follows:\n", "\n", "\\begin{align}\n", "&X= \\sum_{n=1}^N\\left(\\frac{a}{r}\\right)^{n+2}\\left[ g_n^0X_n^0+\\sum_{m=1}^n\\left( g_n^m\\cos m\\phi+h_n^m\\sin m\\phi \\right)X_n^m\\right]\\\\[6pt]\n", "&Y= \\sum_{n=1}^N\\left(\\frac{a}{r}\\right)^{n+2} \\sum_{m=1}^n \\left(g_n^m\\sin m\\phi-h_n^m\\cos m\\phi \\right)Y_n^m \\\\[6pt]\n", "&Z= \\sum_{n=1}^N\\left(\\frac{a}{r}\\right)^{n+2} \\left[g_n^0Z_n^0+\\sum_{m=1}^n\\left( g_n^m\\cos m\\phi+h_n^m\\sin m\\phi \\right)Z_n^m\\right]\\\\[6pt]\n", "\\text{with}&\\\\[6pt]\n", "&X_n^m=\\frac{dP_n^n}{d\\theta}\\\\[6pt]\n", "&Y_n^m=\\frac{m}{\\sin \\theta}P_n^m \\kern{10ex} \\text{(Except at the poles where Y_n^m=X_n^m\\cos \\theta.)}\\\\[6pt]\n", "&Z_n^m=-(n+1)P_n^m\n", "\\end{align}\n", "\n", "\n", "where $n$ and $m$ are spherical harmonic degree and order, respectively, and the ($g_n^m, h_n^m$) are the Gauss coefficients for a particular model (e.g. the IGRF) of maximum degree $N$.\n", "\n", "The Associated Legendre functions of degree $n$ and order $m$ are defined, in Schmidt semi-normalised form by\n", "\n", "$$P^m_n(x) = \\frac{1}{2^n n!}\\left[ \\frac{(2-\\delta_{0m})(n-m)!\\left(1-x^2\\right)^m}{(n+m)!} \\right]^{1/2}\\frac{d^{n+m}}{dx^{n+m}}\\left(1-x^2\\right)^{n},$$\n", "\n", "\n", "where $x = \\cos(\\theta)$. \n", "\n", "Referring to Malin and Barraclough (1981), the recurrence relations\n", "\n", "\\begin{align}\n", "P_n^n&=\\left(1-\\frac{1}{2n}\\right)^{1/2}\\sin \\theta \\thinspace P_{n-1}^{n-1} \\\\[6pt]\n", "P_n^m&=\\left[\\left(2n-1\\right) \\cos \\theta \\thinspace P_{n-1}^m-\\left[ \\left(n-1\\right)^2-m^2\\right]^{1/2}P_{n-2}^m\\right]\\left(n^2-m^2\\right)^{-1/2},\\\\[6pt]\n", "\\end{align}\n", "\n", "and\n", "\n", "\\begin{align}\n", "X_n^n&=\\left(1-\\frac{1}{2n}\\right)^{1/2}\\left( \\sin \\theta \\thinspace X_{n-1}^{n-1}+ \\cos \\theta \\thinspace P_{n-1}^{n-1} \\right)\\\\[6pt]\n", "X_n^m&=\\left[\\left(2n-1\\right) \\cos \\theta \\thinspace X_{n-1}^m- \\sin \\theta \\thinspace P_{n-1}^m\\right] - \\left[ \\left(n-1\\right)^2-m^2\\right]^{1/2}X_{n-2}^m\\left(n^2-m^2\\right)^{-1/2}.\n", "\\end{align}\n", "\n", "may be used to calculate the $X^m_n$, $Y^m_n$ and $Z^m_n$, given $P_0^0=1$, $P_1^1=\\sin(\\theta)$, $X_0^0=0$ and $X_1^1=\\cos(\\theta)$.\n", "

" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting $P_n^m$ and $X_n^m$\n", "The $P_n^m(\\theta)$ and $X_n^m(\\theta)$ are building blocks for computing geomagnetic field models given a spherical harmonic model. It's instructive to visualise these functions and below you can experiment by setting different values of spherical harmonic degree ($n$) and order ($m \\le n$). Note how the choice of $n$ and $m$ affects the number of zeroes of the functions. \n", "\n", "The functions are plotted on a semi-circle representing the surface of the Earth, with the inner core added for cosmetic purposes only! Again, purely for cosmetic purposes, the functions are scaled to fit within $\\pm$10% of the Earth's surface." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "**>> USER INPUT HERE: Set the spherical harmonic degree and order for the plot**" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "degree = 13\n", "order = 7" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate Pnm and Xmn values every 0.5 degrees\n", "colat = np.linspace(0,180,361)\n", "pnmvals = np.zeros(len(colat))\n", "xnmvals = np.zeros(len(colat))\n", "\n", "idx = sha.pnmindex(degree,order)\n", "for i, cl in enumerate(colat):\n", " p,x = sha.pxyznm_calc(degree, cl)[0:2]\n", " pnmvals[i] = p[idx]\n", " xnmvals[i] = x[idx]\n", " \n", "theta = np.deg2rad(colat)\n", "ct = np.cos(theta)\n", "st = np.sin(theta)\n", "\n", "# Numbers mimicking the Earth's surface and outer core radii\n", "e_rad = 6.371\n", "c_rad = 3.485\n", "\n", "# Scale values to fit within 10% of \"Earth's surface\". Firstly the P(n,m),\n", "shell = 0.1*e_rad\n", "pmax = np.abs(pnmvals).max()\n", "pnmvals = pnmvals*shell/pmax + e_rad\n", "xp = pnmvals*st\n", "yp = pnmvals*ct\n", "\n", "# and now the X(n,m)\n", "xmax = np.abs(xnmvals).max()\n", "xnmvals = xnmvals*shell/xmax + e_rad\n", "xx = xnmvals*st\n", "yx = xnmvals*ct\n", "\n", "# Values to draw the Earth's and outer core surfaces as semi-circles\n", "e_xvals = e_rad*st\n", "e_yvals = e_rad*ct\n", "c_xvals = e_xvals*c_rad/e_rad\n", "c_yvals = e_yvals*c_rad/e_rad\n", "\n", "# Earth-like background framework for plots\n", "def eplot(ax):\n", " ax.set_aspect('equal')\n", " ax.set_axis_off()\n", " ax.plot(e_xvals,e_yvals, color='blue')\n", " ax.plot(c_xvals,c_yvals, color='black')\n", " ax.fill_between(c_xvals, c_yvals, y2=0, color='lightgrey')\n", " ax.plot((0, 0), (-e_rad, e_rad), color='black')\n", "\n", "# Plot the P(n,m) and X(n,m)\n", "fig, axes = plt.subplots(nrows=1, ncols=2, figsize=(18, 8))\n", "fig.suptitle('Degree (n) = '+str(degree)+', order (m) = '+str(order), fontsize=20)\n", " \n", "axes.plot(xp,yp, color='red')\n", "axes.set_title('P('+ str(degree)+',' + str(order)+')', fontsize=16)\n", "eplot(axes)\n", "\n", "axes.plot(xx,yx, color='red')\n", "axes.set_title('X('+ str(degree)+',' + str(order)+')', fontsize=16)\n", "eplot(axes)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Try again!**" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## The International Geomagnetic Reference Field\n", "The latest version of the IGRF is IGRF13 which consists of a main-field model every five years from 1900.0 to 2020.0 and a secular variation model for 2020-2025. The main field models have (maximum) spherical harmonic degree (n) and order (m) 10 up to 1995 and n=m=13 from 2000 onwards. The secular variation model has n=m=8.\n", "\n", "The coefficients are first loaded into a pandas database: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "igrf13 = pd.read_csv(IGRF13_FILE, delim_whitespace=True, header=3)\n", "igrf13.head() # Check the values have loaded correctly" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### a) Calculating geomagnetic field values using the IGRF\n", "The function below calculates geomagnetic field values at a point defined by its colatitude, longitude and altitude, using a spherical harmonic model of maximum degree _nmax_ supplied as an array _gh_. The parameter _coord_ is a string specifying whether the input position is in geocentric coordinates (when _altitude_ should be the geocentric distance in km) or geodetic coordinates (when altitude is distance above mean sea level in km). \n", "\n", "(It's unconventional, but I've chosen to include a monopole term, set to zero, at index zero in the _gh_ array.)