{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6dc1075b",
   "metadata": {},
   "source": [
    "# Example 5 - Reproducing ArchKalmag14k"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "751de93f",
   "metadata": {},
   "source": [
    "In this notebook we show how to reproduce the model ArchKalmag14k.r and import it in pymagglobal. We also distribute ArchKalmag14k.r this way. Additional access is provided via a [dedicated website](https://ionocovar.agnld.uni-potsdam.de/Kalmag/Archeo/)."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e55228b",
   "metadata": {},
   "source": [
    "## Set up the data\n",
    "\n",
    "We first set up the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40431586",
   "metadata": {},
   "outputs": [],
   "source": [
    "import io\n",
    "import requests\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "# load the rejection list\n",
    "rejection_response = requests.get(\n",
    "    'https://nextcloud.gfz-potsdam.de/s/'\n",
    "    'WLxDTddq663zFLP/download/rejection_list.npz',\n",
    ")\n",
    "rejection_response.raise_for_status()\n",
    "\n",
    "with np.load(io.BytesIO(rejection_response.content), allow_pickle=True) as fh:\n",
    "    to_reject = fh['to_reject']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b651f720",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paleokalmag import ChunkedData\n",
    "from paleokalmag.data_handling import read_data\n",
    "\n",
    "data_response = requests.get(\n",
    "    'https://nextcloud.gfz-potsdam.de/s/'\n",
    "    'r6YxrrABRJjideS/download/archkalmag_data.csv',\n",
    ")\n",
    "data_response.raise_for_status()\n",
    "# Since we load the data from the web, we read it manually\n",
    "# and pass the dataframe to ChunkedData\n",
    "data = read_data(\n",
    "    io.BytesIO(data_response.content),\n",
    "    rejection_lists=to_reject,\n",
    ")\n",
    "\n",
    "cdat = ChunkedData(\n",
    "    data,\n",
    "    delta_t=10,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ebacd2b7",
   "metadata": {},
   "source": [
    "## Set up and run the model\n",
    "\n",
    "We next set up a model with the estimated hyperparameters. They have been estimated using the procedure illustrated in [example 3](https://sec23.git-pages.gfz-potsdam.de/korte/paleokalmag/example_3.html)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d0e4e88d",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paleokalmag.corefieldmodel import CoreFieldModel\n",
    "\n",
    "a14kModel = CoreFieldModel(\n",
    "    lmax=20,\n",
    "    R=2800,\n",
    "    gamma=-38.00,\n",
    "    alpha_dip=39.00,\n",
    "    tau_dip=171.34,\n",
    "    alpha_wodip=118.22,\n",
    "    tau_wodip=379.59,\n",
    "    rho=2.25,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67e51c62",
   "metadata": {},
   "source": [
    "Next we run the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "cd22259e",
   "metadata": {},
   "outputs": [],
   "source": [
    "a14kModel.forward(cdat, output_every=5, quiet=False)\n",
    "a14kModel.save('./ArchKalmag14k_forward.npz')\n",
    "a14kModel.backward(quiet=False)\n",
    "a14kModel.save('./ArchKalmag14k_smoothed.npz')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a16f9caa",
   "metadata": {},
   "source": [
    "## Use the model with pymagglobal\n",
    "\n",
    "Finally, we load the model and compare it to other models by using [pymagglobal](https://sec23.git-pages.gfz-potsdam.de/korte/pymagglobal)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "742212ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "from paleokalmag import PaleoKalmag\n",
    "\n",
    "archkalmag14k = PaleoKalmag(\n",
    "    './ArchKalmag14k_smoothed.npz',\n",
    "    name='ArchKalmag14k',\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6f64528a",
   "metadata": {},
   "outputs": [],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "from pymagglobal import Model, dipole_series\n",
    "\n",
    "plt.rcParams['font.size'] = '14'\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots(1, 1, figsize=(10, 6))\n",
    "\n",
    "a14ktimes = np.linspace(-12000, 1900, 1001)   \n",
    "dip, a16, a84 = dipole_series(a14ktimes, archkalmag14k)\n",
    "ax.plot(a14ktimes, dip, color='C0', zorder=2, label=archkalmag14k.name)\n",
    "ax.fill_between(\n",
    "    a14ktimes,\n",
    "    a16,\n",
    "    a84,\n",
    "    color='C0', \n",
    "    alpha=0.4,\n",
    "    zorder=0,\n",
    ")\n",
    "\n",
    "arch10k = Model('ARCH10k.1')\n",
    "# different times due to different valid interval\n",
    "a10ktimes = np.linspace(-8000, 1900, 1001)\n",
    "dip = dipole_series(a10ktimes, arch10k)\n",
    "ax.plot(a10ktimes, dip, color='C1', zorder=2, label=arch10k.name)\n",
    "\n",
    "ax.legend();\n",
    "ax.set_ylabel(r'Dipole moment [$10^{22}$Am$^2$]');\n",
    "ax.set_xlabel(r'time [yrs.]');"
   ]
  }
 ],
 "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.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}