{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "681fb4de",
   "metadata": {},
   "source": [
    "# Generate synthetic data"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c938d6b",
   "metadata": {},
   "source": [
    "Here we briefly show how `pymagglobal` can be used to generate synthetic data. We first set up the model we want to use to generate the synthetic data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f52a366b",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymagglobal import Model\n",
    "\n",
    "myModel = Model('CALS10k.2')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13988168",
   "metadata": {},
   "source": [
    "Next we have to generate a data distribution. We use `pymagglobal`s `get_grid` routine, to generate `n_at` random locations, that are uniformly distributed on the sphere. Times are drawn uniformly as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f57b72a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymagglobal.utils import get_grid\n",
    "\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "from cartopy import crs as ccrs\n",
    "\n",
    "# the number of artificial records\n",
    "n_at = 400\n",
    "grid = get_grid(n_at, random=True)\n",
    "# set the times to something arbitrary\n",
    "grid[3] = np.random.uniform(-1000, 1900, size=n_at)\n",
    "\n",
    "# plot the records, convert colat to lat\n",
    "fig, ax = plt.subplots(1, 1, subplot_kw={'projection': ccrs.Mollweide()})\n",
    "ax.scatter(grid[1], 90-grid[0], transform=ccrs.PlateCarree())\n",
    "ax.set_global()\n",
    "ax.coastlines();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "806081dc",
   "metadata": {},
   "source": [
    "Finally we can evaluate the model at the inputs, to get synthetic ''observations'' of the field"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5842b198",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymagglobal.core import field\n",
    "obs = field(grid, myModel)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "14502d90",
   "metadata": {},
   "source": [
    "Let's have a look at the records:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "67cd6ed7",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},\n",
    "                        figsize=(15, 4))\n",
    "titles = [r'$B_N$', r'$B_E$', r'$B_Z$']\n",
    "for it in range(3):\n",
    "    axs[it].set_title(titles[it])\n",
    "    axs[it].scatter(\n",
    "        grid[1], \n",
    "        90 - grid[0], \n",
    "        c=obs[it], transform=ccrs.PlateCarree(),\n",
    "    )\n",
    "    axs[it].set_global()\n",
    "    axs[it].coastlines();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e96b4ec7",
   "metadata": {},
   "source": [
    "Declination, inclination and intensity records are also easily generated, by using the `field` kwargs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57b38352",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs_dif = field(grid, myModel, field_type='dif')\n",
    "fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},\n",
    "                        figsize=(15, 4))\n",
    "titles = [r'$D$', r'$I$', r'$F$']\n",
    "for it in range(3):\n",
    "    axs[it].set_title(titles[it])\n",
    "    axs[it].scatter(grid[1], 90-grid[0], c=obs_dif[it], transform=ccrs.PlateCarree())\n",
    "    axs[it].set_global()\n",
    "    axs[it].coastlines();"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a16058b9",
   "metadata": {},
   "source": [
    "##  Data from an arbitray set of coefficients"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "003f37e3",
   "metadata": {},
   "source": [
    "Sometimes it is helpful to generate data not from a model that's included in pymagglobal, but from an arbitrary set of coefficients. Provided they are in the ''standard order'', i.e. $g_1^0, g_1^1, g_1^{-1}, g_2^0, g_2^1, g_2^{-1}, ...$, this is also straight-forward using utility routines from pymagglobal:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fdbdfd1c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymagglobal.core import coefficients\n",
    "\n",
    "epoch = -3000\n",
    "# First generate coefficients at epoch\n",
    "_, _, coeffs = coefficients(epoch, myModel)\n",
    "\n",
    "# Generate new input locations at the given epoch\n",
    "grid = get_grid(n_at, random=True, t=epoch)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e81a742d",
   "metadata": {},
   "source": [
    "We use the pymagglobal routine `dsh_basis`, which evaluates the derivatives of the spherical harmonics up to a given degree at given inputs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5924dde3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymagglobal.utils import dsh_basis\n",
    "\n",
    "base = dsh_basis(myModel.l_max, grid)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d875959",
   "metadata": {},
   "source": [
    "Field values are then generated via a dot product:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a0a43e60",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = coeffs @ base"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c804ae2",
   "metadata": {},
   "source": [
    "The outputs are now given by obs, every third value is N, E, Z. To have a more convenient form, we reshape the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6eb5bc3",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = obs.reshape(n_at, 3).T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d3ffc91a",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(1, 3, subplot_kw={'projection': ccrs.Mollweide()},\n",
    "                        figsize=(15, 4))\n",
    "titles = [r'$B_N$', r'$B_E$', r'$B_Z$']\n",
    "for it in range(3):\n",
    "    axs[it].set_title(titles[it])\n",
    "    axs[it].scatter(grid[1], 90-grid[0], c=obs[it], transform=ccrs.PlateCarree())\n",
    "    axs[it].set_global()\n",
    "    axs[it].coastlines();"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}