{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Radial Profile"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example shows how to create a radial profile from a SOXS event file, including using an exposure map to get flux-based quantities. We'll simulate a simple isothermal cluster."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib\n",
    "\n",
    "matplotlib.rc(\"font\", size=18)\n",
    "import matplotlib.pyplot as plt\n",
    "import soxs\n",
    "import astropy.io.fits as pyfits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, create the spectrum for the cluster using an absorbed thermal APEC model:"
   ]
  },
  {
   "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": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "kT = 6.0\n",
    "abund = 0.3\n",
    "redshift = 0.05\n",
    "norm = 1.0\n",
    "spec = agen.get_spectrum(kT, abund, redshift, norm)\n",
    "spec.rescale_flux(1.0e-13, emin=0.5, emax=2.0, flux_type=\"energy\")\n",
    "spec.apply_foreground_absorption(0.02)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And a spatial distribution based on a $\\beta$-model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "pos = soxs.BetaModel(30.0, 45.0, 50.0, 0.67)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Generate a SIMPUT catalog from these two models, and write it to a file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "width = 10.0  # arcmin by default\n",
    "nx = 1024  # resolution of image\n",
    "cluster = soxs.SimputSpectrum.from_models(\"beta_model\", spec, pos, width, nx)\n",
    "cluster_cat = soxs.SimputCatalog.from_source(\n",
    "    \"beta_model.simput\", cluster, overwrite=True\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and run the instrument simulation (for simplicity we'll turn off the point-source background):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.instrument_simulator(\n",
    "    \"beta_model.simput\",\n",
    "    \"evt.fits\",\n",
    "    (100.0, \"ks\"),\n",
    "    \"lynx_hdxi\",\n",
    "    [30.0, 45.0],\n",
    "    overwrite=True,\n",
    "    ptsrc_bkgnd=False,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make an exposure map so that we can obtain flux-based quantities:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.make_exposure_map(\"evt.fits\", \"expmap.fits\", 2.3, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Make the radial profile, using energies between 0.5 and 5.0 keV, between radii of 0 and 200 arcseconds, with 50 bins:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "soxs.write_radial_profile(\n",
    "    \"evt.fits\",\n",
    "    \"profile.fits\",\n",
    "    [30.0, 45.0],\n",
    "    0,\n",
    "    200,\n",
    "    50,\n",
    "    emin=0.5,\n",
    "    emax=5.0,\n",
    "    expmap_file=\"expmap.fits\",\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can use AstroPy's FITS reader to open the profile and have a look at the columns that are inside:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f = pyfits.open(\"profile.fits\")\n",
    "f[\"PROFILE\"].columns"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and use Matplotlib to plot some quantities. We can plot the surface brightness:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "plt.errorbar(\n",
    "    f[\"profile\"].data[\"rmid\"],\n",
    "    f[\"profile\"].data[\"sur_bri\"],\n",
    "    lw=2,\n",
    "    yerr=f[\"profile\"].data[\"sur_bri_err\"],\n",
    ")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"r (arcsec)\")\n",
    "plt.ylabel(\"S (cts/s/arcsec**2)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and, since we used an exposure map, the surface flux:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 8))\n",
    "plt.errorbar(\n",
    "    f[\"profile\"].data[\"rmid\"],\n",
    "    f[\"profile\"].data[\"sur_flux\"],\n",
    "    lw=2,\n",
    "    yerr=f[\"profile\"].data[\"sur_flux_err\"],\n",
    ")\n",
    "plt.xscale(\"log\")\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"r (arcsec)\")\n",
    "plt.ylabel(\"S (cts/s/cm**2/arcsec**2)\")"
   ]
  }
 ],
 "metadata": {
  "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
}