{
"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",
"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": [
"## 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": [
"## 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": [
"## 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.<br>"
]
},
{
"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: <br>\n",
"D = declination <br>\n",
"H = horizontal intensity <br>\n",
"I = inclination <br>\n",
"X = north component <br>\n",
"Y = east component <br>\n",
"Z = vertically (downwards) component <br>\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.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}