{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Magnetic Poles\n", "## drift between 2000-01-01 and 2030-01-01\n", "\n", "The calucation uses the Fortran libraty of Emmert et al.[1] and wrapped by the [apexpy](https://apexpy.readthedocs.io/en/latest/readme.html) Python package.\n", "\n", "This version also calculates the azimut angle of the prime magnetic meridian at the poles locations.\n", "\n", "[1] Emmert, J. T., A. D. Richmond, and D. P. Drob (2010),\n", " A computationally compact representation of Magnetic-Apex\n", " and Quasi-Dipole coordinates with smooth base vectors,\n", " J. Geophys. Res., 115(A8), A08322,\n", " [doi:10.1029/2010JA015326](http://dx.doi.org/10.1029/2010JA015326)" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from apexpy import Apex\n", "\n", "def qd2geo(time, qdlat, qdlon, height=0, precision=1e-5):\n", " \"\"\" Get geographic coordinates for the given QD coordinates. \"\"\"\n", " apex = Apex(time)\n", " lat, lon, error = apex.qd2geo(qdlat, qdlon, height, precision)\n", " if error > precision:\n", " raise RuntimeError(\n", " \"The result precision %gdeg is not within the required limit %gdeg!\"\n", " \"\" % (error, precision)\n", " )\n", " return lat, lon" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "from datetime import datetime\n", "from itertools import product\n", "\n", "start_year = 2000\n", "end_year = 2030\n", "\n", "times = [datetime(year, 1, 1) for year in range(start_year, end_year + 1)]\n", "\n", "eps = 1 # deg\n", "\n", "pole_north = [qd2geo(time, 90, 0) for time in times]\n", "pole_south = [qd2geo(time, -90, 0) for time in times]\n", "pole_north_ngh = [qd2geo(time, 90 - eps, 180) for time in times]\n", "pole_south_ngh = [qd2geo(time, -90 + eps, 0) for time in times]\n", "\n", "prime_meridian = [\n", " qd2geo(times[len(times)//2], lat , 0) \n", " for lat in range(-89, 90)\n", "]\n", "\n", "anti_meridian = [\n", " qd2geo(times[len(times)//2], lat , 180) \n", " for lat in range(-89, 90)\n", "]" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "from math import pi\n", "from numpy import array, dot, arccos, sqrt, tan, cross, arctan2\n", "from eoxmagmod import vnorm, convert,GEOCENTRIC_SPHERICAL, GEOCENTRIC_CARTESIAN\n", "\n", "def get_tanget_vector(cart0, cart1):\n", " \"\"\" Get tangent vector at cart0 pointing to cart1.\n", " cart0 and cart1 are Cartesian coordinates of two distinct points\n", " on shell of a unit Sphere.\n", " \"\"\"\n", " angle = arccos(dot(cart0, cart1))\n", " cart1_tan = cart1 * sqrt(1.0 + tan(angle)**2)\n", " d_cart = cart1_tan - cart0\n", " return d_cart / vnorm(d_cart)\n", "\n", "\n", "def calculate_azimuth(coord0, coord1):\n", " \"\"\" Calculate local azimuth at coord0 pointing to coord1.\n", " coords0 and coord1 are Spherical (lat/lon) coordinates of two\n", " distinct points on a Sphere.\n", " \"\"\"\n", " lat0, lon0 = coord0\n", " lat1, lon1 = coord1\n", " lat2, lon2 = lat0 + 0.1, lon0 # north\n", " if lat2 > 90:\n", " lat2 = 180 - lat2, lon2 + 180\n", " print(lat2, lon2)\n", " cart0, cart1, cart2 = convert(\n", " [(lat0, lon0, 1.0), (lat1, lon1, 1.0), (lat0 + 0.1, lon0, 1.0)],\n", " GEOCENTRIC_SPHERICAL, GEOCENTRIC_CARTESIAN\n", " )\n", " cart_dir = get_tanget_vector(cart0, cart1)\n", " cart_north = get_tanget_vector(cart0, cart2)\n", " cart_east = cross(cart_north, cart0)\n", " return arctan2(dot(cart_dir, cart_east), dot(cart_dir, cart_north)) * 180 / pi\n", "\n", "\n", "azimuth_north = [\n", " calculate_azimuth(np, npm) for np, npm in zip(pole_north, pole_north_ngh)\n", "]\n", "\n", "azimuth_south = [\n", " calculate_azimuth(sp, spm) for sp, spm in zip(pole_south, pole_south_ngh)\n", "]" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/html": [ "
\n", " | date | \n", "northMagPoleLat | \n", "northMagPoleLon | \n", "northPrimeMeridAzimuth | \n", "southMagPoleLat | \n", "southPoleLon | \n", "southPrimeMeridAzimuth | \n", "
---|---|---|---|---|---|---|---|
0 | \n", "2000-01-01 | \n", "81.866020 | \n", "-82.987709 | \n", "-8.001296 | \n", "-74.128334 | \n", "126.434753 | \n", "160.460766 | \n", "
1 | \n", "2001-01-01 | \n", "81.938301 | \n", "-83.054321 | \n", "-8.021357 | \n", "-74.138107 | \n", "126.396301 | \n", "160.447508 | \n", "
2 | \n", "2002-01-01 | \n", "82.010452 | \n", "-83.121162 | \n", "-8.041747 | \n", "-74.147903 | \n", "126.357841 | \n", "160.433509 | \n", "
3 | \n", "2003-01-01 | \n", "82.082474 | \n", "-83.188232 | \n", "-8.062271 | \n", "-74.157692 | \n", "126.319382 | \n", "160.420156 | \n", "
4 | \n", "2004-01-01 | \n", "82.154373 | \n", "-83.255516 | \n", "-8.083352 | \n", "-74.167488 | \n", "126.280930 | \n", "160.406628 | \n", "
5 | \n", "2005-01-01 | \n", "82.226151 | \n", "-83.323059 | \n", "-8.104661 | \n", "-74.177299 | \n", "126.242462 | \n", "160.392930 | \n", "
6 | \n", "2006-01-01 | \n", "82.306572 | \n", "-83.440826 | \n", "-8.147176 | \n", "-74.198517 | \n", "126.196213 | \n", "160.364156 | \n", "
7 | \n", "2007-01-01 | \n", "82.386887 | \n", "-83.559349 | \n", "-8.190690 | \n", "-74.219749 | \n", "126.150002 | \n", "160.335197 | \n", "
8 | \n", "2008-01-01 | \n", "82.467102 | \n", "-83.678673 | \n", "-8.235127 | \n", "-74.241013 | \n", "126.103851 | \n", "160.306180 | \n", "
9 | \n", "2009-01-01 | \n", "82.547211 | \n", "-83.798813 | \n", "-8.280491 | \n", "-74.262299 | \n", "126.057755 | \n", "160.277201 | \n", "
10 | \n", "2010-01-01 | \n", "82.627205 | \n", "-83.919785 | \n", "-8.326688 | \n", "-74.283615 | \n", "126.011719 | \n", "160.248113 | \n", "
11 | \n", "2011-01-01 | \n", "82.716835 | \n", "-84.043762 | \n", "-8.375749 | \n", "-74.303459 | \n", "125.956863 | \n", "160.224074 | \n", "
12 | \n", "2012-01-01 | \n", "82.806381 | \n", "-84.168839 | \n", "-8.426153 | \n", "-74.323334 | \n", "125.902046 | \n", "160.200052 | \n", "
13 | \n", "2013-01-01 | \n", "82.895844 | \n", "-84.295052 | \n", "-8.477822 | \n", "-74.343239 | \n", "125.847267 | \n", "160.175624 | \n", "
14 | \n", "2014-01-01 | \n", "82.985222 | \n", "-84.422462 | \n", "-8.530882 | \n", "-74.363159 | \n", "125.792534 | \n", "160.151249 | \n", "
15 | \n", "2015-01-01 | \n", "83.074516 | \n", "-84.551102 | \n", "-8.585155 | \n", "-74.383118 | \n", "125.737839 | \n", "160.127008 | \n", "
16 | \n", "2016-01-01 | \n", "83.158089 | \n", "-84.696404 | \n", "-8.633016 | \n", "-74.400307 | \n", "125.661491 | \n", "160.100657 | \n", "
17 | \n", "2017-01-01 | \n", "83.241623 | \n", "-84.842804 | \n", "-8.681990 | \n", "-74.417557 | \n", "125.585182 | \n", "160.074246 | \n", "
18 | \n", "2018-01-01 | \n", "83.325127 | \n", "-84.990356 | \n", "-8.732282 | \n", "-74.434868 | \n", "125.508926 | \n", "160.047718 | \n", "
19 | \n", "2019-01-01 | \n", "83.408600 | \n", "-85.139091 | \n", "-8.783854 | \n", "-74.452225 | \n", "125.432732 | \n", "160.021112 | \n", "
20 | \n", "2020-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
21 | \n", "2021-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
22 | \n", "2022-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
23 | \n", "2023-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
24 | \n", "2024-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
25 | \n", "2025-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
26 | \n", "2026-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
27 | \n", "2027-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
28 | \n", "2028-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
29 | \n", "2029-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "
30 | \n", "2030-01-01 | \n", "83.492035 | \n", "-85.289047 | \n", "-8.836948 | \n", "-74.469643 | \n", "125.356590 | \n", "159.994392 | \n", "