{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# pyXSIM Example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To show how to make a set of photons from a 3D dataset using pyXSIM and yt for \n",
    "reading into SOXS, we'll look at is that of thermal emission from a galaxy \n",
    "cluster. In this case, the gas in the core of the cluster is \"sloshing\" in \n",
    "the center, producing spiral-shaped cold fronts. The dataset we want to use \n",
    "for this example is available for download from the [yt Project](http://yt-project.org) \n",
    "at [this link](http://yt-project.org/data/GasSloshing.tar.gz). \n",
    "\n",
    "First, import our necessary modules:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import yt\n",
    "import pyxsim\n",
    "import soxs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we `load` the dataset with yt (this dataset does not have species fields, \n",
    "so we specify that the gas is fully ionized in this case so that the emission measure\n",
    "field can be computed correctly):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = yt.load(\n",
    "    \"GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150\", default_species_fields=\"ionized\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use yt to take a slice of density and temperature through the center of \n",
    "the dataset so we can see what we're looking at: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "slc = yt.SlicePlot(\n",
    "    ds, \"z\", [(\"gas\", \"density\"), (\"gas\", \"temperature\")], width=(1.0, \"Mpc\")\n",
    ")\n",
    "slc.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, sloshing gas as advertised. Next, we'll create a sphere object to serve as \n",
    "a source for the photons. Place it at the center of the domain with `\"c\"`, and \n",
    "use a radius of 500 kpc:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sp = ds.sphere(\"c\", (0.5, \"Mpc\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, we need to set up the emission model for our source. We said we were going \n",
    "to look at the thermal emission from the hot plasma, so we'll do that here by \n",
    "using `ThermalSourceModel`. The first four arguments are the name of the underlying \n",
    "spectral model, the maximum and minimum energies, and the number of bins in the \n",
    "spectrum. We've chosen these numbers so that the spectrum has an energy resolution \n",
    "of about 1 eV. Setting `thermal_broad=True` turns on thermal broadening. This \n",
    "simulation does not include metallicity, so we'll do something simple and say \n",
    "that it uses the above spectral model and the metallicity is a constant $Z = 0.3~Z_\\odot$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "source_model = pyxsim.CIESourceModel(\n",
    "    \"apec\", 0.5, 9.0, 9000, thermal_broad=True, Zmet=0.3\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We're almost ready to go to generate the photons from this source, but first \n",
    "we should decide what our redshift, collecting area, and exposure time should \n",
    "be. Let's pick big numbers, because remember the point of this first step is \n",
    "to create a Monte-Carlo sample from which to draw smaller sub-samples for mock \n",
    "observations. Note these are all (value, unit) tuples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "exp_time = (300.0, \"ks\")  # exposure time\n",
    "area = (3.0, \"m**2\")  # collecting area\n",
    "redshift = 0.2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, that's everything--let's create the photons!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_photons, n_cells = pyxsim.make_photons(\n",
    "    \"my_photons\", sp, redshift, area, exp_time, source_model\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, that was easy. Now we have a photon list that we can use to create events, \n",
    "using the `project_photons()` function. To be realistic, we're going to want to assume \n",
    "foreground Galactic absorption, using the \"TBabs\" absorption model and assuming a \n",
    "foreground absorption column of $N_H = 4 \\times 10^{20}~{\\rm cm}^{-2}$. Here we'll \n",
    "just do a simple projection along the z-axis, reducing the exposure time, and \n",
    "centering the photons at RA, Dec = (30, 45) degrees:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_events = pyxsim.project_photons(\n",
    "    \"my_photons\", \"my_events\", \"z\", (30.0, 45.0), absorb_model=\"tbabs\", nH=0.04\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can then use this event list that we wrote as an input to the instrument simulator in SOXS. We'll use a smaller exposure time (100 ks instead of 500 ks), and observe it with the Lynx calorimeter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.instrument_simulator(\n",
    "    \"my_events.h5\",\n",
    "    \"evt.fits\",\n",
    "    (100.0, \"ks\"),\n",
    "    \"lynx_lxm\",\n",
    "    [30.0, 45.0],\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use the `write_image()` function in SOXS to bin the events into an image \n",
    "and write them to a file, restricting the energies between 0.5 and 2.0 keV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.write_image(\"evt.fits\", \"img.fits\", emin=0.5, emax=2.0, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can show the resulting image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = soxs.plot_image(\n",
    "    \"img.fits\", stretch=\"sqrt\", cmap=\"afmhot\", vmax=1000.0, width=0.05\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also bin the events into a spectrum using `write_spectrum()` and write \n",
    "the spectrum to disk:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.write_spectrum(\"evt.fits\", \"evt.pha\", overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and plot the spectrum using `plot_spectrum()`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = soxs.plot_spectrum(\"evt.pha\", xmin=0.5, xmax=7.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's zoom into the region of the spectrum around the iron line to look at \n",
    "the detailed structure afforded by the resolution of the calorimeter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ax.set_xlim(5.4, 5.7)\n",
    "fig"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}