{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Thermal Emission"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The first example we'll look at is that of thermal emission from a galaxy cluster. In this case, the gas in the core of the cluster is \"sloshing\" in the center, producing spiral-shaped cold fronts. The dataset we want to use for this example is available for download from the [yt Project](http://yt-project.org) 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": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "import yt\n",
    "import pyxsim\n",
    "import soxs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Next, we `load` the dataset with yt. Note that this dataset does not have \n",
    "species fields in it, so we'll set `default_species_fields=\"ionized\"` to \n",
    "assume full ionization (as appropriate for galaxy clusters):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "ds = yt.load(\n",
    "    \"GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150\", default_species_fields=\"ionized\"\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Let's use yt to take a slice of density and temperature through the center of the dataset so we can see what we're looking at: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "slc = yt.SlicePlot(\n",
    "    ds, \"z\", [(\"gas\", \"density\"), (\"gas\", \"temperature\")], width=(1.0, \"Mpc\")\n",
    ")\n",
    "slc.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Ok, sloshing gas as advertised. Next, we'll create a sphere object to serve as a source for our investigations. Place it at the center of the domain with `\"c\"`, and use a radius of 500 kpc:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "sp = ds.sphere(\"c\", (500.0, \"kpc\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now, we need to set up a source model. We said we were going to look at the thermal emission from the hot plasma, which in a galaxy cluster is in collisional ionization equilibrium (CIE), so to do that we can set up a `CIESourceModel`. The first argument specifies which model we want to use. Let's use `\"spex\"`. The next three arguments are the maximum and minimum energies, and the number of bins in the spectrum. We've chosen these numbers so that the spectrum has an energy resolution of about 1 eV. \n",
    "\n",
    "`CIESourceModel` takes a lot of optional arguments, which you can investigate in the docs, but here we'll do something simple and say that the metallicity is a constant $Z = 0.3~Z_\\odot$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "source_model = pyxsim.CIESourceModel(\"spex\", 0.05, 11.0, 1000, 0.3, binscale=\"log\")"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "We can use this `source_model` object to compute X-ray fields for use in yt calculations. For example, we can create fields for emissivity, luminosity, and photon emissivity within a particular band, in this case 0.5-7 keV:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "xray_fields = source_model.make_source_fields(ds, 0.5, 7.0)\n",
    "print(xray_fields)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "From this we can compute the total luminosity in our sphere object:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "print(sp[\"gas\", \"xray_luminosity_0.5_7.0_keV\"])\n",
    "print(sp.sum((\"gas\", \"xray_luminosity_0.5_7.0_keV\")))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "and we can make a projection of the emissivity:"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%% md\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "prj = yt.ProjectionPlot(\n",
    "    ds, \"z\", (\"gas\", \"xray_emissivity_0.5_7.0_keV\"), width=(1.0, \"Mpc\")\n",
    ")\n",
    "prj.show()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We're almost ready to go to generate the photons from this source, but first we should decide what our redshift, collecting area, and exposure time should be. Let's pick big numbers, because remember the point of this first step is to create a Monte-Carlo sample from which to draw smaller sub-samples for mock observations. Note these are all (value, unit) tuples:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [
    "exp_time = (300.0, \"ks\")  # exposure time\n",
    "area = (1000.0, \"cm**2\")  # collecting area\n",
    "redshift = 0.05"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "So, that's everything--let's create the photons! We use the `make_photons` function for this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "n_photons, n_cells = pyxsim.make_photons(\n",
    "    \"sloshing_photons\", sp, redshift, area, exp_time, source_model\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Ok, that was easy. Now we have a photon list that we can use to create events using the `project_photons` function. Here, we'll just do a simple projection along the z-axis, and center the photons at RA, Dec = (45, 30) degrees. Since we want to be realistic, we'll want to apply foreground galactic absorption using the `\"tbabs\"` model, assuming a neutral hydrogen column of $N_H = 4 \\times 10^{20}~{\\rm cm}^{-2}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "n_events = pyxsim.project_photons(\n",
    "    \"sloshing_photons\",\n",
    "    \"sloshing_events\",\n",
    "    \"z\",\n",
    "    (45.0, 30.0),\n",
    "    absorb_model=\"tbabs\",\n",
    "    nH=0.04,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now that we have a set of \"events\" on the sky, we can use them as an input to the instrument simulator in SOXS. We'll use a small exposure time (100 ks instead of 300 ks), and observe it with the as-launched ACIS-I model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "soxs.instrument_simulator(\n",
    "    \"sloshing_events.h5\",\n",
    "    \"evt.fits\",\n",
    "    (100.0, \"ks\"),\n",
    "    \"chandra_acisi_cy0\",\n",
    "    [45.0, 30.0],\n",
    "    overwrite=True,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "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": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "soxs.write_image(\"evt.fits\", \"img.fits\", emin=0.5, emax=2.0, overwrite=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we can take a quick look at the image: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "pycharm": {
     "name": "#%%\n",
     "is_executing": true
    }
   },
   "outputs": [],
   "source": [
    "soxs.plot_image(\n",
    "    \"img.fits\", stretch=\"sqrt\", cmap=\"arbre\", vmin=0.0, vmax=10.0, width=0.2\n",
    ")"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "interpreter": {
   "hash": "4d0667920515a25babb24d7efbecfb85048fc348313e91f189b3ebd8f3c9393c"
  },
  "kernelspec": {
   "display_name": "Python 3.10.4 ('base')",
   "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}