{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import chandra_limits as cl\n",
    "from kadi.commands import read_backstop, states\n",
    "import matplotlib as mpl\n",
    "import matplotlib.pyplot as plt\n",
    "import xija\n",
    "from xija.get_model_spec import get_xija_model_spec\n",
    "from cxotime import CxoTime\n",
    "import numpy as np\n",
    "from pprint import pprint\n",
    "import logging\n",
    "from ska_matplotlib import pointpair\n",
    "logging.getLogger(\"xija\").setLevel(logging.ERROR)\n",
    "mpl.rc(\"font\", size=18)\n",
    "mpl.rc(\"axes\", linewidth=2)\n",
    "mpl.rc(\"xtick.major\", width=2)\n",
    "mpl.rc(\"ytick.major\", width=2)\n",
    "mpl.rc(\"xtick.minor\", width=2)\n",
    "mpl.rc(\"ytick.minor\", width=2)\n",
    "np.set_printoptions(threshold=20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic Usage"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For this example, we will get states from a backstop file, but in theory\n",
    "they could come from the ``kadi`` states database, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs_cmds = read_backstop(\"/data/acis/LoadReviews/2022/MAR2822/ofls/CR086_2107.backstop\")\n",
    "sts = states.get_states(cmds=bs_cmds, merge_identical=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The first example is a simple and boring one, the 1DEAMZT upper limit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the DEALimit object\n",
    "dea_limit = cl.DEALimit()\n",
    "# Get the \"high\" limit, with states as input (though for this case they are not used)\n",
    "dea_limit_line = dea_limit.get_limit_line(sts, which='high')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``dea_limit_line`` is a ``LimitLine`` object with various useful\n",
    "attributes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(dea_limit_line.times) # the times of the limit line\n",
    "print(dea_limit_line.values) # the values of the limit line\n",
    "print(dea_limit_line.reasons) # the reason for the limit at particular times"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the line by itself (we’ll make this more useful below):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = dea_limit_line.plot(lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A more interesting limit is the zero-FEPs 1DPAMZT limit, which is a lower limit and only applies if 0 FEPs are on:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpa_limit = cl.DPALimit()\n",
    "dpa_limit_line = dpa_limit.get_limit_line(sts, which='low')\n",
    "print(dpa_limit_line.values)\n",
    "print(dpa_limit_line.reasons)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that in this case there is no limit at certain times. We can also see this if we plot it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = dpa_limit_line.plot(lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The ACIS Focal Plane Limit"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The ACIS FP limit is more complex, as it depends on a number of conditions pertaining to the schedule of observations. ACIS observations have different limits depending on which array of chips is used, whether or not there is a grating inserted, and the number of expected counts in the observation. The information about the limit values can be determined by querying the obscat, or they can be read directly from an OR list. \n",
    "\n",
    "This limit is handled by the `ACISFPLimit` class. Before it can be used to create a limit line, one needs to run the `ACISFPLimit.set_obs_info` method, which takes as a required input a dictionary of four lists and/or NumPy arrays: `obsid`, `start_science`, `stop_science`, and `simode`. These are the observation IDs, start times, stop times, and SIMODEs of the observations in the schedule. Note that the start and stop times are defined by the times of the ACIS `startScience` and `stop_science` commands, respectively. Note also that the SIMODE names must reflect whether or not the bias should be calculated (e.g., `TE_0057AB` with bias, and `TE_0057A` without bias). This dictionary can be constructed using any method (so long as the OBSIDs are valid). One way to make one easily is to use the `determine_obsid_info` helper function, which takes a states table as input (from ``kadi``, for example) and determines the OBSIDs, their start and stop times, and the SIMODEs from the states: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bs_cmds2 = read_backstop(\"/data/acis/LoadReviews/2024/JUL1524/ofls/CR196_1802.backstop\")\n",
    "sts2 = states.get_states(cmds=bs_cmds2, merge_identical=True)\n",
    "obs_list = cl.determine_obsid_info(sts2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting structure is very simple (showing here only a portion of \n",
    "the output):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pprint(obs_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this `obs_list` in hand, one can determine the ACIS FP limit line in a very similar manner to the other limits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acisfp_limit = cl.ACISFPLimit()\n",
    "acisfp_limit.set_obs_info(obs_list)\n",
    "fp_limit_line = acisfp_limit.get_limit_line(sts2)\n",
    "fig, ax = fp_limit_line.plot(lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the limit values for the different observations are determined by an obscat query. However, if an OR list is available, the limit values can be read directly from it like so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acisfp_limit = cl.ACISFPLimit()\n",
    "acisfp_limit.set_obs_info(\n",
    "    obs_list,\n",
    "    orlist=\"/data/acis/LoadReviews/2024/JUL1524/ofls/mps/or/JUL1524_A.or\"\n",
    ")\n",
    "fp_limit_line = acisfp_limit.get_limit_line(sts2)\n",
    "fig, ax = fp_limit_line.plot(lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the two limit lines obtained by either method are identical, as\n",
    "they should be.\n",
    "\n",
    "If you want to plot the limit so that the different limit reasons have\n",
    "different colors, you can set `use_colors=True`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = fp_limit_line.plot(lw=2, use_colors=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running with Models and Checking for Violations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we show an example where we run a thermal model, check the limit against it, and see if there are violations. We’ll use the 1DPAMZT model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "dpa_json = get_xija_model_spec(\"dpa\")[0]\n",
    "tm_dpa = xija.XijaModel(\"dpa\", sts[\"datestart\"][0], sts[\"datestop\"][-1], \n",
    "                        model_spec=dpa_json, cmd_states=sts)\n",
    "tm_dpa.comp[\"dpa0\"].set_data(10.0)\n",
    "tm_dpa.make()\n",
    "tm_dpa.calc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We feed the states and the times from the model, and then check the violations using the `check_violations` method of the `LimitLine` object. In this case, we’ll check the low-temperature limit:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "dpa_limit_line = dpa_limit.get_limit_line(sts, which='low')\n",
    "viols = dpa_limit_line.check_violations(tm_dpa)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Two violations were found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pprint(viols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can double-check the first one by eye by plotting the model, the FEP count, and the limit line. We can also use the helper function `plot_viols` to plot bands where the violations occur:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "ax2 = ax.twinx()\n",
    "x = pointpair(CxoTime(sts[\"datestart\"]).plot_date, \n",
    "              CxoTime(sts[\"datestop\"]).plot_date)\n",
    "y = pointpair(sts[\"fep_count\"])\n",
    "ax2.plot(x, y, lw=2, color='magenta', drawstyle='steps')\n",
    "ax.plot(CxoTime(tm_dpa.times).plot_date, tm_dpa.comp[\"1dpamzt\"].mvals, lw=2)\n",
    "dpa_limit_line.plot(lw=5, fig_ax=(fig, ax), color=\"C1\")\n",
    "ax.set_xlim(CxoTime(\"2022:086:21:00:00\").plot_date, \n",
    "            CxoTime(\"2022:087:04:00:00\").plot_date)\n",
    "ax.set_ylim(5, 20)\n",
    "ax.set_ylabel(\"1DPAMZT ($^\\circ$C)\")\n",
    "cl.plot_viols(ax, viols, alpha=0.25)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can try the same thing for the focal plane model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "acisfp_json = get_xija_model_spec(\"acisfp\")[0]\n",
    "tm = xija.XijaModel(\"acisfp\", sts2[\"datestart\"][0], sts2[\"datestop\"][-1], \n",
    "                    model_spec=acisfp_json, cmd_states=sts2)\n",
    "tm.make()\n",
    "tm.calc()\n",
    "\n",
    "obs_list = cl.determine_obsid_info(sts2)\n",
    "acisfp_limit = cl.ACISFPLimit(margin=0.0)\n",
    "acisfp_limit.set_obs_info(obs_list)\n",
    "fp_limit_line = acisfp_limit.get_limit_line(sts, which='high')\n",
    "viols = fp_limit_line.check_violations(tm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case, many violations were found:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pprint(viols)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot some of these also:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = plt.subplots(figsize=(10,10))\n",
    "ax.plot(CxoTime(tm.times).plot_date, tm.comp[\"fptemp\"].mvals, lw=2)\n",
    "fp_limit_line.plot(lw=2, use_colors=True, fig_ax=(fig, ax))\n",
    "ax.set_xlim(CxoTime(\"2024:200:00:00:00\").plot_date, \n",
    "            CxoTime(\"2024:204:00:00:00\").plot_date)\n",
    "ax.set_ylim(-121, -97)\n",
    "ax.set_ylabel(\"FPTEMP_11 ($^\\circ$C)\")\n",
    "cl.plot_viols(ax, viols, alpha=0.25)"
   ]
  }
 ],
 "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.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}