{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Two Clusters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the SOXS Python interface, this example shows how to generate photons from two thermal spectra and two $\\beta$-model spatial distributions, as an approximation of two galaxy clusters. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "\n",
    "matplotlib.rc(\"font\", size=18)\n",
    "import soxs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Spectral Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We want to generate thermal spectra, so we first create a spectral generator using the ``ApecGenerator`` class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "emin = 0.05  # keV\n",
    "emax = 20.0  # keV\n",
    "nbins = 20000\n",
    "agen = soxs.ApecGenerator(emin, emax, nbins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we'll generate the two thermal spectra. We'll set the APEC norm for each to 1, and renormalize them later:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kT1 = 6.0\n",
    "abund1 = 0.3\n",
    "redshift1 = 0.05\n",
    "norm1 = 1.0\n",
    "spec1 = agen.get_spectrum(kT1, abund1, redshift1, norm1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kT2 = 4.0\n",
    "abund2 = 0.4\n",
    "redshift2 = 0.1\n",
    "norm2 = 1.0\n",
    "spec2 = agen.get_spectrum(kT2, abund2, redshift2, norm2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, re-normalize the spectra using energy fluxes between 0.5-2.0 keV:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux1 = 1.0e-13  # erg/s/cm**2\n",
    "flux2 = 5.0e-14  # erg/s/cm**2\n",
    "emin = 0.5  # keV\n",
    "emax = 2.0  # keV\n",
    "spec1.rescale_flux(flux1, emin=0.5, emax=2.0, flux_type=\"energy\")\n",
    "spec2.rescale_flux(flux2, emin=0.5, emax=2.0, flux_type=\"energy\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also apply foreground galactic absorption to each spectrum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_H = 0.04  # 10^20 atoms/cm^2\n",
    "spec1.apply_foreground_absorption(n_H)\n",
    "spec2.apply_foreground_absorption(n_H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "``spec1`` and ``spec2`` are ``Spectrum`` objects. Let's have a look at the spectra:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = spec1.plot(\n",
    "    xmin=0.1,\n",
    "    xmax=20.0,\n",
    "    ymin=1.0e-8,\n",
    "    ymax=1.0e-3,\n",
    "    label=\"$\\mathrm{kT\\ =\\ 6\\ keV,\\ Z\\ =\\ 0.3\\ Z_\\odot,\\ z\\ =\\ 0.05}$\",\n",
    ")\n",
    "spec2.plot(\n",
    "    label=\"$\\mathrm{kT\\ =\\ 4\\ keV,\\ Z\\ =\\ 0.4\\ Z_\\odot,\\ z\\ =\\ 0.1}$\", fig=fig, ax=ax\n",
    ")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate Spatial Models"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now what we want to do is associate spatial distributions with these spectra. Each cluster will be represented using a $\\beta$-model. For that, we use the ``BetaModel`` class. For fun, we'll give the second ``BetaModel`` an ellipticity and tilt it by 45 degrees (a bit extreme, but it demonstrates the functionality nicely):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters for the clusters\n",
    "r_c1 = 30.0  # in arcsec\n",
    "r_c2 = 20.0  # in arcsec\n",
    "beta1 = 2.0 / 3.0\n",
    "beta2 = 1.0\n",
    "\n",
    "# Center of the field of view\n",
    "ra0 = 30.0  # degrees\n",
    "dec0 = 45.0  # degrees\n",
    "\n",
    "# Space the clusters roughly a few arcminutes apart on the sky.\n",
    "# They're at different redshifts, so one is behind the other.\n",
    "dx = 3.0 / 60.0  # degrees\n",
    "ra1 = ra0 - 0.5 * dx\n",
    "dec1 = dec0 - 0.5 * dx\n",
    "ra2 = ra0 + 0.5 * dx\n",
    "dec2 = dec0 + 0.5 * dx\n",
    "\n",
    "# Now actually create the models\n",
    "pos1 = soxs.BetaModel(ra1, dec1, r_c1, beta1, ellipticity=0.5, theta=45.0)\n",
    "pos2 = soxs.BetaModel(ra2, dec2, r_c2, beta2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Create SIMPUT files"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now, what we want to do is generate energies and positions from these models. We want to create a large sample that we'll draw from when we run the instrument simulator, so we choose a large exposure time and a large collecting area (should be bigger than the maximum of the ARF). To do this, we use the `from_models()` method of the `SimputPhotonList` class to make instances of the latter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_exp = (500.0, \"ks\")\n",
    "area = (3.0, \"m**2\")\n",
    "cluster_phlist1 = soxs.SimputPhotonList.from_models(\n",
    "    \"cluster1\", spec1, pos1, t_exp, area\n",
    ")\n",
    "cluster_phlist2 = soxs.SimputPhotonList.from_models(\n",
    "    \"cluster2\", spec2, pos2, t_exp, area\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can quickly show the positions using the `plot()` method of the `SimputPhotonList` instances. For simplicity, we'll only show every 100th event using the ``stride`` argument, and restrict ourselves to a roughly $20'\\times~20'$ field of view."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = cluster_phlist1.plot(\n",
    "    [30.0, 45.0], 6.0, marker=\".\", stride=100, label=\"Cluster 1\"\n",
    ")\n",
    "cluster_phlist2.plot(\n",
    "    [30.0, 45.0], 6.0, marker=\".\", stride=100, fig=fig, ax=ax, label=\"Cluster 2\"\n",
    ")\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the positions and the energies of the photons in the `SimputPhotonList`s, we can write them to a SIMPUT catalog, using the `SimputCatalog` class. Each cluster will have its own photon list, but be part of the same SIMPUT catalog:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the SIMPUT catalog \"sim_cat\" from the photon lists \"cluster1\" and \"cluster2\"\n",
    "sim_cat = soxs.SimputCatalog.from_source(\n",
    "    \"clusters_simput.fits\", cluster_phlist1, overwrite=True\n",
    ")\n",
    "sim_cat.append(cluster_phlist2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate an Observation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can use the instrument simulator to simulate the two clusters by ingesting the SIMPUT file, setting an output file `\"evt.fits\"`, setting an exposure time of 50 ks (less than the one we used to generate the source), the `\"lynx_hdxi\"` instrument, and the pointing direction of (RA, Dec) = (30.,45.) degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.instrument_simulator(\n",
    "    \"clusters_simput.fits\",\n",
    "    \"evt.fits\",\n",
    "    (50.0, \"ks\"),\n",
    "    \"lynx_hdxi\",\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 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": [
    "Now we show the resulting image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fig, ax = soxs.plot_image(\n",
    "    \"img.fits\", stretch=\"log\", cmap=\"viridis\", vmin=0.1, vmax=10.0, width=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Alternative Way to Generate the SIMPUT Catalog"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the above example, we generated the SIMPUT catalog for the observation of the two clusters using `SimputPhotonList`s, which in previous versions was the only option available in SOXS. It is also possible to use two `SimputSpectrum` objects, which is another type of SIMPUT source that consists of a spectrum and (optionally) an image. The image is used by SOXS to serve as a model to generate photon positions on the sky. If no image is included, then the source is simply a point source. \n",
    "\n",
    "In this case of course, the clusters are two extended sources, so we can use the `from_models` method in a similar way as we did above, but in this case we have to supply the `width` and the resolution (`nx`) of the image that we want to associate with the spectrum:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "width = 10.0  # arcmin by default\n",
    "nx = 1024  # resolution of image\n",
    "cluster_spec1 = soxs.SimputSpectrum.from_models(\"cluster1\", spec1, pos1, width, nx)\n",
    "cluster_spec2 = soxs.SimputSpectrum.from_models(\"cluster2\", spec2, pos2, width, nx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we create the SIMPUT catalog in essentially the same way as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create the SIMPUT catalog \"sim_cat\" from the spectra \"cluster1\" and \"cluster2\" in the same way\n",
    "sim_cat2 = soxs.SimputCatalog.from_source(\n",
    "    \"clusters2_simput.fits\", cluster_spec1, overwrite=True\n",
    ")\n",
    "sim_cat2.append(cluster_spec2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the `instrument_simulator`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.instrument_simulator(\n",
    "    \"clusters2_simput.fits\",\n",
    "    \"evt2.fits\",\n",
    "    (50.0, \"ks\"),\n",
    "    \"lynx_hdxi\",\n",
    "    [30.0, 45.0],\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and make an image:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.write_image(\"evt2.fits\", \"img2.fits\", emin=0.5, emax=2.0, overwrite=True)\n",
    "fig, ax = soxs.plot_image(\n",
    "    \"img2.fits\", stretch=\"log\", cmap=\"viridis\", vmin=0.1, vmax=10.0, width=0.1\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We used the same models, so the resulting images are the same except that different random numbers were used."
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}