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

IAGA Summer School 2019

\n", "\n", "

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", "sys.path.append('..')\n", "from src import sha_lib as sha, mag_lib as mag\n", "\n", "IGRF12_FILE = os.path.abspath('../data/external/igrf12coeffs.txt')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 1. 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).\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)$. The following 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", "(e.g. Malin and Barraclough, 1981) may be used to calculate the $X^m_n$, $Y^m_n$ and $Z^m_n$, for example\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}$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 2. 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": [ "### >> 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 = 5" ] }, { "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[0].plot(xp,yp, color='red')\n", "axes[0].set_title('P('+ str(degree)+',' + str(order)+')', fontsize=16)\n", "eplot(axes[0])\n", "\n", "axes[1].plot(xx,yx, color='red')\n", "axes[1].set_title('X('+ str(degree)+',' + str(order)+')', fontsize=16)\n", "eplot(axes[1])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# 3. The International Geomagnetic Reference Field\n", "The latest version of the IGRF is IGRF12 which consists of a main-field model every five years from 1900.0 to 2015.0 and a secular variation model for 2015-2020. The main field models have 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": [ "igrf12 = pd.read_csv(IGRF12_FILE, delim_whitespace=True, header=3)\n", "igrf12.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.
" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def shm_calculator(gh, nmax, altitude, colat, long, coord):\n", " RREF = 6371.2 #The reference radius assumed by the IGRF\n", " degree = nmax\n", " phi = long\n", "\n", " if (coord == 'Geodetic'):\n", " # Geodetic to geocentric conversion using the WGS84 spheroid\n", " rad, theta, sd, cd = sha.gd2gc(altitude, colat)\n", " else:\n", " rad = altitude\n", " theta = colat\n", "\n", " # Function 'rad_powers' to create an array with values of (a/r)^(n+2) for n = 0,1, 2 ..., degree\n", " rpow = sha.rad_powers(degree, RREF, rad)\n", "\n", " # Function 'csmphi' to create arrays with cos(m*phi), sin(m*phi) for m = 0, 1, 2 ..., degree\n", " cmphi, smphi = sha.csmphi(degree,phi)\n", "\n", " # Function 'gh_phi_rad' to create arrays with terms such as [g(3,2)*cos(2*phi) + h(3,2)*sin(2*phi)]*(a/r)**5 \n", " ghxz, ghy = sha.gh_phi_rad(gh, degree, cmphi, smphi, rpow)\n", "\n", " # Function 'pnm_calc' to calculate arrays of the Associated Legendre Polynomials for n (&m) = 0,1, 2 ..., degree\n", " pnm, xnm, ynm, znm = sha.pxyznm_calc(degree, theta)\n", "\n", " # Geomagnetic field components are calculated as a dot product\n", " X = np.dot(ghxz, xnm)\n", " Y = np.dot(ghy, ynm)\n", " Z = -np.dot(ghxz, znm)\n", "\n", " # Convert back to geodetic (X, Y, Z) if required\n", " if (coord == 'Geodetic'):\n", " t = X\n", " X = X*cd + Z*sd\n", " Z = Z*cd - t*sd\n", "\n", " return((X, Y, Z))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### >> >> USER INPUT HERE: Set the input parameters" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "location = 'Erehwon'\n", "ctype = 'Geodetic' # coordinate type\n", "altitude = 0 # in km above the spheroid if ctype = 'Geodetic', radial distance if ctype = 'Geocentric'\n", "colat = 35 # NB colatitude, not latitude\n", "long = -3 # longitude\n", "NMAX = 13 # Maxiimum spherical harmonic degree of the model\n", "date = 2015.0 # Date for the field estimates" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now calculate the IGRF geomagnetic field estimates." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Calculate the gh values for the supplied date\n", "if date == 2015.0:\n", " gh = igrf12['2015.0']\n", "elif date < 2015.0:\n", " date_1 = (date//5)*5\n", " date_2 = date_1 + 5\n", " w1 = date-date_1\n", " w2 = date_2-date\n", " gh = np.array((w2*igrf12[str(date_1)] + w1*igrf12[str(date_2)])/(w1+w2))\n", "elif date > 2015.0:\n", " gh =np.array(igrf12['2015.0'] + (date-2015.0)*igrf12['2015-20'])\n", "\n", "gh = np.append(0., gh) # Add a zero monopole term corresponding to g(0,0)\n", "\n", "bxyz = shm_calculator(gh, NMAX, altitude, colat, long, ctype)\n", "dec, hoz ,inc , eff = mag.xyz2dhif(bxyz[0], bxyz[1], bxyz[2])\n", "\n", "print('\\nGeomagnetic field values at: ', location+', '+ str(date), '\\n')\n", "print('Declination (D):', '{: .1f}'.format(dec), 'degrees')\n", "print('Inclination (I):', '{: .1f}'.format(inc), 'degrees')\n", "print('Horizontal intensity (H):', '{: .1f}'.format(hoz), 'nT')\n", "print('Total intensity (F) :', '{: .1f}'.format(eff), 'nT')\n", "print('North component (X) :', '{: .1f}'.format(bxyz[0]), 'nT')\n", "print('East component (Y) :', '{: .1f}'.format(bxyz[1]), 'nT')\n", "print('Vertical component (Z) :', '{: .1f}'.format(bxyz[2]), 'nT')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# b) Maps of the IGRF\n", "Now draw maps of the IGRF at the date selected above. The latitude range is set at -85 degrees to +85 degrees and the longitude range -180 degrees to +180 degrees and IGRF values for (X, Y, Z) are calculated on a 5 degree grid (this may take a few seconds to complete).\n", "\n", "## >> >> USER INPUT HERE: Set the element to plot\n", "Enter the geomagnetic element to plot below:
\n", "D = declination
\n", "H = horizontal intensity
\n", "I = inclination
\n", "X = north component
\n", "Y = east component
\n", "Z = vertically (downwards) component
\n", "F = total intensity.)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "el2plot = 'H'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "def IGRF_plotter(el_name, vals, date):\n", " if el_name=='D':\n", " cvals = np.arange(-25,30,5)\n", " else:\n", " cvals = 15\n", " fig, ax = plt.subplots(figsize=(16, 8))\n", " cplt = ax.contour(longs, lats, vals, levels=cvals)\n", " ax.clabel(cplt, cplt.levels, inline=True, fmt='%d', fontsize=10)\n", " ax.set_title('IGRF: '+ el_name + ' (' + str(date) + ')', fontsize=20)\n", " ax.set_xlabel('Longitude', fontsize=16)\n", " ax.set_ylabel('Latitude', fontsize=16)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "longs = np.linspace(-180, 180, 73)\n", "lats = np.linspace(-85, 85, 35)\n", "Bx, By, Bz = zip(*[sha.shm_calculator(gh,13,6371.2,90-lat,lon,'Geocentric') \\\n", " for lat in lats for lon in longs])\n", "X = np.asarray(Bx).reshape(35,73)\n", "Y = np.asarray(By).reshape(35,73)\n", "Z = np.asarray(Bz).reshape(35,73)\n", "D, H, I, F = [mag.xyz2dhif(X, Y, Z)[el] for el in range(4)]\n", "\n", "el_dict={'X':X, 'Y':Y, 'Z':Z, 'D':D, 'H':H, 'I':I, 'F':F}\n", "IGRF_plotter(el2plot, el_dict[el2plot], date)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "Produce a similar plot for the secular variation according to the IGRF." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### References\n", "\n", "Malin, S. R. . and Barraclough, D., (1981). An algorithm for synthesizing the geomagnetic field, Computers & Geosciences. Pergamon, 7(4), pp. 401–405. doi: 10.1016/0098-3004(81)90082-0." ] } ], "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 }