{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "Studying Ptolemy's sun declination using visualization tools\n",
    "----"
   ],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "from kanon.units import Sexagesimal\n",
    "from kanon.tables import HTable\n",
    "import math\n",
    "from astropy.units import arcsecond\n",
    "from matplotlib import pyplot as plt\n",
    "import numpy as np"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# This function redefines what is in the declination notebook\n",
    "\n",
    "def build_declination(sin_table, obliquity):\n",
    "    arcsin_table = sin_table.copy(set_index=\"Val\")\n",
    "    obl = sin_table.get(obliquity)\n",
    "    obl_table = sin_table.apply(\"Val\", lambda x: x * obl)\n",
    "    return obl_table.apply(\"Val\", lambda x: round(arcsin_table.get(x), 2))"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "OBLIQUITY = \"23;51,20\""
   ],
   "outputs": [],
   "metadata": {
    "tags": [
     "parameters"
    ]
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "sin = lambda x: math.sin(math.radians(x))\n",
    "asin = lambda x: math.degrees(math.asin(x))\n",
    "\n",
    "# This table holds computed sine values\n",
    "x = list(Sexagesimal.range(91))\n",
    "y = [round(Sexagesimal.from_float(sin(n), 3)) for n in x]\n",
    "sin_table_true = HTable([x, y], names=(\"Arg\", \"Val\"), index=\"Arg\")\n",
    "\n",
    "# This table holds computed declination values\n",
    "sin_obl = sin(Sexagesimal(OBLIQUITY))\n",
    "x = list(Sexagesimal.range(1, 91))\n",
    "y = [round(Sexagesimal.from_float(asin(sin(n) * sin_obl), 2)) for n in x]\n",
    "decl_table_true = HTable(\n",
    "    [x, y],\n",
    "    names=(\"Arg\", \"Val\"),\n",
    "    index=\"Arg\",\n",
    ")\n"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# We import Ptolemy's sun declination table from DISHAS\n",
    "ptolemy_decl = HTable.read(214)\n",
    "\n",
    "ptolemy_decl.plot2d()\n",
    "plt.show()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# We rebuild a declination table from a sine table without odd arguments\n",
    "\n",
    "sin_table_grid2 = sin_table_true[::2]\n",
    "decl_table_grid2 = build_declination(sin_table_grid2, Sexagesimal(OBLIQUITY))[1:]"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "# Let's compare Ptolemy's declination with our reconstructed tables\n",
    "\n",
    "residue_nogrid = ptolemy_decl.apply(\n",
    "    \"Entries\", lambda x: x - decl_table_true[\"Val\"], \"Declination residue\"\n",
    ")\n",
    "residue_nogrid[\"Declination residue\"] = residue_nogrid[\"Declination residue\"].to(\n",
    "    arcsecond\n",
    ")\n",
    "# For the computed declination table residue, we plot in blue\n",
    "residue_nogrid.plot2d(\"bx\", linestyle=\"dashed\", lw=0.4)\n",
    "\n",
    "residue_grid2 = ptolemy_decl[1::2].apply(\n",
    "    \"Entries\", lambda x: x - decl_table_grid2[\"Val\"], \"Declination residue\"\n",
    ")\n",
    "residue_grid2[\"Declination residue\"] = residue_grid2[\"Declination residue\"].to(\n",
    "    arcsecond\n",
    ")\n",
    "# For the built declination with grid 2 sine residue, we plot in red\n",
    "residue_grid2.plot2d(\"r+\", linestyle=\"dashed\", lw=0.4)\n",
    "\n",
    "plt.show()"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [
    "arr_nogrid = residue_nogrid[\"Declination residue\"].astype(float)\n",
    "arr_grid2 = residue_grid2[\"Declination residue\"].astype(float)\n",
    "\n",
    "print(\"Residue with no interpolation\")\n",
    "print(f\"mean : {np.mean(arr_nogrid):.4f}, std : {np.std(arr_nogrid):.4f}\")\n",
    "print(f\"quartiles : {np.quantile(arr_nogrid, [0.25, 0.5, 0.75])} \\n\")\n",
    "print(\"Residue with grid 2 interpolation\")\n",
    "print(f\"mean : {np.mean(arr_grid2):.4f}, std : {np.std(arr_grid2):.4f}\")\n",
    "print(f\"quartiles : {np.quantile(arr_grid2, [0.25, 0.5, 0.75])}\")"
   ],
   "outputs": [],
   "metadata": {}
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "source": [],
   "outputs": [],
   "metadata": {}
  }
 ],
 "metadata": {
  "kernelspec": {
   "name": "python3",
   "display_name": "Python 3.8.6 64-bit"
  },
  "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.8.6"
  },
  "metadata": {
   "interpreter": {
    "hash": "4cd7ab41f5fca4b9b44701077e38c5ffd31fe66a6cab21e0214b68d958d0e462"
   }
  },
  "interpreter": {
   "hash": "4cd7ab41f5fca4b9b44701077e38c5ffd31fe66a6cab21e0214b68d958d0e462"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}