{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Magnetic apex coordinates and MLT\n",
    "\n",
    "I will show how to calculate magnetic apex (QD) coordinates, QD vector compontents, and magnetic local time. I will use apexpy for the coordinates and pyamps for MLT, since the MLT code in pyamps is faster than in apexpy.\n",
    "\n",
    "Start by defining some dummy coordinates and measurements:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "glat, glon = np.array([70, 30, -60, 80]), np.zeros(4) # geodetic latitude and longitude\n",
    "\n",
    "B = np.random.random((len(glat), 10000, 2)) # 10000 measurements, 2 components, from each station"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "set up Apex object - epoch 2010, and calculate coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import apexpy\n",
    "\n",
    "a = apexpy.Apex(2010)\n",
    "qdlat, qdlon = a.geo2qd(glat, glon, 0) # last parameter is height\n",
    "\n",
    "# get QD base vectors:\n",
    "f1, f2 = a.basevectors_qd(glat, glon, 0, coords = 'geo')\n",
    "\n",
    "# calculate F parameters:\n",
    "F = f1[0]*f2[1] - f1[1]*f2[0]\n",
    "\n",
    "# calcualte QD components:\n",
    "Eqd = f1[0].reshape((-1, 1)) * B[:, :, 0] / F.reshape((-1, 1)) + f1[1].reshape((-1, 1)) * B[:, :, 1] / F.reshape((-1, 1))\n",
    "Nqd = f2[0].reshape((-1, 1)) * B[:, :, 0] / F.reshape((-1, 1)) + f2[1].reshape((-1, 1)) * B[:, :, 1] / F.reshape((-1, 1))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calcualte magnetic local time using pyamps:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyamps.mlt_utils import mlon_to_mlt\n",
    "import pandas as pd # need to set up some dummy times:\n",
    "\n",
    "t = pd.date_range(start = '2010-01-01 00:00', end = '2012-01-01 00:00', periods = B.shape[1])\n",
    "\n",
    "# loop through each location and calculate the MLT for each time (same length as )\n",
    "mlt = np.array([mlon_to_mlt(lon, t, epoch = 2010) for lon in qdlon])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "done!"
   ]
  },
  {
   "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}