ACTNotebook_Section_4_visualize_objects.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Section 4: The DR5 Cluster Catalog and Visualizing Objects in the Maps\n",
    "\n",
    "This notebook will demonstrate how to use a catalog to visualize objects in the maps as well as how to stack on the locations of the objects in that catalog.  In order to look at these objects we will start by reading in a cluster catalog from [Hilton et al (2020)](https://arxiv.org/abs/2009.11043) that has been released as part of ACT DR5.  The catalog contains 4195 optically confirmed galaxy clusters that were obtained using the coadded maps discussed in [Section 2](Section_2_coadded_maps.ipynb).\n",
    "\n",
    "We ask that when using the catalog you reference the original work: [Hilton et al (2020)](https://arxiv.org/abs/2009.11043)\n",
    "\n",
    "\n",
    "Notebook Info | Value\n",
    "---|---\n",
    "Dependencies | `numpy`, `scipy`, `healpy`, `matplotlib`, `pixell`, `astropy`, `plotly`\n",
    "Data products | \"act_planck_dr5.01_s08s18_AA_f150_night_map_dg_I.fits\"\n",
    " | \"act_planck_dr5.01_s08s18_AA_f150_night_map_d56_I.fits\"\n",
    " | \"DR5_cluster-catalog_v1.0b3.fits\"\n",
    " | \"messier_objs.csv\"\n",
    "Optional Data Products | \"act_planck_dr5.01_s08s18_AA_f150_night_map.fits\"\n",
    "Memory usage | Minimal unless using the optional full resolution and area map\n",
    "Contributors | Kyra Fichman, Maya Mallaby-Kay, Zach Atkins\n",
    "\n",
    "To return to the main notebook or to view a different section use the links in the table of contents. See the main ACT notebook repo (https://github.com/ACTCollaboration/DR4_DR5_Notebooks) for the full set of things to do.\n",
    "\n",
    "## Table of contents\n",
    "\n",
    " --- \n",
    "> [Section 1](ACTNotebook_Section_1_Introduction.ipynb): Overview and Introduction\n",
    "\n",
    "> [Section 4](ACTNotebookSection_4_visualize_objects.ipynb):  The DR5 Cluster Catalog from [Hilton et al (2020)](https://arxiv.org/abs/2009.11043) and Visualizing Objects in the Maps\n",
    "\n",
    "> [Section 5](ACTNotebookSection_5_Lensing_maps.ipynb):  Introduction to the Lensing maps from [Darwish et al (2020)](https://arxiv.org/pdf/2004.01139.pdf)\n",
    "\n",
    " ---"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Initialize notebook\n",
    "\n",
    "Again, we rerun the first notebook in order to import packages and set the data path.  From there we will read in the coadded map.  Users can choose to import the lower resolution, full-size map along with the full-resolution cutout, or the full resolution, full-size map depending on their computer's capabilities.  To that end make the appropriate choice for your computer and edit the code if need be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%run ACTNotebook_Section_1_Introduction.ipynb\n",
    "\n",
    "from scipy.interpolate import interp1d\n",
    "from astropy.io import ascii, fits\n",
    "import astropy.coordinates as coord\n",
    "import astropy.table as atpy\n",
    "import astropy.units as u\n",
    "import random \n",
    "import os, sys\n",
    "import plotly\n",
    "import plotly.graph_objs as go\n",
    "\n",
    "\n",
    "# Redefine our plotting function from section 2\n",
    "def eshow(x,**kwargs): \n",
    "    ''' Define a function to help us plot the maps neatly '''\n",
    "    plots = enplot.get_plots(x, **kwargs)\n",
    "    enplot.show(plots, method = \"ipython\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As in notebook 2 in the original ACT Notebook repo, due to the large size of the coadded map, we give you two options when working with these maps:\n",
    "\n",
    ">1) You can read in the full sky-area, downgraded intensity map, along with a full resolution cutout of the sky.  Between the two maps you'll be able to see the scope of our data as well as the depth.\n",
    "\n",
    ">2) Alternatively, if you have space and resources for the full file, but don't want to read in the whole map, we demonstrate how to read in just a small area of the full, higher-resolution map.\n",
    "\n",
    "Here we default to the first option but feel free to change to the second by switching which cell is commented out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Comment this block and run the next one if you are using the full maps\n",
    "\n",
    "#Option 1: Read in the downgraded map and the full-res cutout\n",
    "\n",
    "imap_dg = enmap.read_map(path + \"act_planck_dr5.01_s08s18_AA_f150_night_map_dg_I.fits\")\n",
    "\n",
    "imap_d56 = enmap.read_map(path + \"act_planck_dr5.01_s08s18_AA_f150_night_map_d56_I.fits\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This second box is commented out currently but if you have the full maps downloaded then feel free to read that in instead of the maps in the cell above.\n",
    "\n",
    "One thing worth looking at here is that the `[0]` in `imap = enmap.read_map(path + \"act_planck_dr5.01_s08s18_AA_f150_night_map.fits\",   box = box)[0]` is there because the map contains the temperature map and then the Q and U maps and (for now) we only want to read in the first component which is the temperature map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# # Uncomment this block to use the full map\n",
    "\n",
    "# # Option 2: Read in a subset of the full map\n",
    "\n",
    "# # Set the size of the box in degrees and convert to radians\n",
    "# dec_from, dec_to = np.deg2rad([-30,5])\n",
    "# ra_from, ra_to = np.deg2rad([-10, 50])\n",
    "\n",
    "# # Create the box\n",
    "# box = [[dec_from,ra_from],[dec_to,ra_to]]\n",
    "\n",
    "# # read in the map\n",
    "# imap = enmap.read_map(path + \"act_planck_dr5.01_s08s18_AA_f150_night_map.fits\", box = box)[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Purely to make it easier to check which maps we should be working with in the next steps we run a couple of try except statements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# We check which version of the maps we have read in so that later on we can use the right ones\n",
    "try:\n",
    "    imap_dg\n",
    "except NameError:\n",
    "    imap_dg = None\n",
    "\n",
    "try:\n",
    "    imap_d56\n",
    "except NameError:\n",
    "    imap_d56 = None"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Begin Analysis"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "PF92UVTxtOh6"
   },
   "source": [
    "In the previous sections you were probably able to see a few bright objects, which were galaxy clusters, just by visually examining the map. However it would be great to be able to look at these more closely. One option is to use a catalog of clusters and compare the catalog to the map. In order to do so we need to get the RA and DEC of the relevant sources and then cut a box surrounding that location.\n",
    "\n",
    "Here we use a catalog that is available through LAMBDA from ACT DR5:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "kCj6mwWQP3px",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "# Read in ras and decs from a cluster catalog\n",
    "tab = atpy.Table().read(path + 'DR5_cluster-catalog_v1.0b3.fits')\n",
    "\n",
    "# We also pull the Ra and Dec information so we have that easily available for later on\n",
    "ras = tab['RADeg']\n",
    "decs = tab['decDeg']\n",
    "\n",
    "# We're going to plot the objects in the catalog according to their mass\n",
    "# but feel free to choose a different attribute to highlight\n",
    "mass = tab['M500']\n",
    "mass[np.where(mass<0)] = 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualizing the objects in the DR5 Catalog\n",
    "\n",
    "To begin with lets take a look at a 3D representation of the clusters in this catalog.  \n",
    "\n",
    "In order to do so we first set our fiducial cosmology using astropy.  Setting the cosmology allows us to easily calculate the comoving distance to objects in the catalog using their redshifts.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.cosmology import WMAP9 as cosmo\n",
    "\n",
    "# This astropy function will give us a distance in Mpc\n",
    "r = cosmo.comoving_distance(tab['redshift'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have our comoving distances we need to get the angles to the objects and then convert these quantities to x, y, z coordinates.  We do that in the following cell:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the theta, and phi angles from RA and Dec\n",
    "theta=np.radians(tab['RADeg'])\n",
    "phi=np.radians(90-tab['decDeg'])\n",
    "\n",
    "# Use the comoving distance, theta and phi to get the x, y, z coordinates for the objects in Mpc\n",
    "x=r*np.cos(theta)*np.sin(phi)\n",
    "y=r*np.sin(theta)*np.sin(phi)\n",
    "z=r*np.cos(phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's plot our objects in 3D to look at the distribution of clusters in space!  \n",
    "\n",
    "We're using plotly here which makes really awesome interactive plots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Plotly to be rendered inline in the notebook.\n",
    "plotly.offline.init_notebook_mode()\n",
    "\n",
    "# Configure the trace.\n",
    "trace = go.Scatter3d(\n",
    "    x=x, y=y, z=z,\n",
    "    mode='markers',\n",
    "    marker={'size':2*mass, 'opacity': 0.8,}, \n",
    "    hovertext=  [\"Mass: {} Trillion Solar Masses <br>Redshift: {} \"\n",
    "                 .format(round(tab['M500'][x], 2)*1e2, \n",
    "                 round(tab[\"redshift\"][x], 3)) \n",
    "                 for x in range(len(tab['M500'])) ] )\n",
    "\n",
    "# Configure the layout.\n",
    "layout = go.Layout(margin={'l': 0, 'r': 0, 'b': 0, 't': 0})\n",
    "data = [trace]\n",
    "fig = go.Figure(data=data)\n",
    "\n",
    "\n",
    "fig.update_layout(scene = dict(\n",
    "                    xaxis_title='Mpc',\n",
    "                    yaxis_title='Mpc',\n",
    "                    zaxis_title='Mpc'))\n",
    "\n",
    "# Render the plot.\n",
    "plotly.offline.iplot(fig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Clusters in our CMB maps\n",
    "Now that we've looked at the clusters let's try looking at cutouts of these objects in our CMB maps."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "mG5bAZF4P-H1"
   },
   "source": [
    "As noted previously, one way to look at a source is to build a \"bounding box\" and to call the `submap` method.  Let's make the box slightly smaller this time in order to hone in on the objects we are interested in.  \n",
    "\n",
    "This time we will use matplotlib to plot.  This makes substantially less pretty plots but is much faster - if you'd like to use the eshow function feel free to switch back to that!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 336
    },
    "colab_type": "code",
    "id": "amE6DZdmP-aZ",
    "outputId": "1d7db7bc-556b-4dc4-e8c7-71c39606a7f5"
   },
   "outputs": [],
   "source": [
    "# Choose an object from the catalog\n",
    "# Choose an object from the catalog\n",
    "# --------  Try Changing this number to look at different clusters -------- \n",
    "n = 128   # <<<<<<<<\n",
    "\n",
    "ra, dec = [ras[n],decs[n]]\n",
    "width = 30/60.\n",
    "\n",
    "# Build the box around the object and cut out a submap from the larger map\n",
    "box = np.deg2rad([[dec-width/2.,ra-width/2.],[dec+width/2.,ra+width/2.]])\n",
    "\n",
    "if imap_d56 is not None:\n",
    "    stamp = imap_d56.submap(box)\n",
    "else:\n",
    "    stamp = imap.submap(box)\n",
    "\n",
    "# Visualize the submap\n",
    "plt.imshow(stamp, cmap = cm.coolwarm)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "7wuon2PcQDH3"
   },
   "source": [
    "## Stacking clusters and using reproject\n",
    "\n",
    "If we want to stack objects in the map it's relatively easy to do so, but we need to use a slightly different approach. Instead of using a bounding box and cutting a submap from the map we can use pixell's `reproject.thumbnails` function to cut out \"stamps\" that are reprojected at the center of each object. \n",
    "\n",
    "Reprojecting means that we avoid stacking clusters that are distorted due to the rectangular projection of the overall map, as the reprojection will automatically center the \"stamp\" basis on each individual point source. Users can specify the type of projection using the \"proj\" keyword. The default is a \"gnomonic\" projection which projects to a plane that's tangential to the center source so that any distortions are radially symmetric."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Reproject the above source and plot - again we need to check if we're using the d56 or the full map \n",
    "\n",
    "# Set the size of the stamp:\n",
    "r_arcmin = 10.  #Define the radius of the image in arcminutes\n",
    "r = r_arcmin*np.pi/(60 * 180)  # convert arcminutes to radians \n",
    "\n",
    "if imap_d56 is not None:\n",
    "    stamp = reproject.thumbnails(imap_d56, coords = np.deg2rad([decs[n], ras[n]]), res =  0.0000454440, r = r)\n",
    "else:\n",
    "    stamp = reproject.thumbnails(imap, coords = np.deg2rad([decs[n], ras[n]]), res =  0.0000454440, r = r)\n",
    "\n",
    "\n",
    "extent = [-r_arcmin/60., r_arcmin/60., -r_arcmin/60., r_arcmin/60.]\n",
    "plt.imshow(stamp, extent=extent,cmap = cm.coolwarm)\n",
    "plt.colorbar()\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Degrees\", fontsize = 14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see the utility of `reproject.thumbnails` in action: the source is now centered in the new coordinates! Now we stack all of the reprojected sources.  Note this may take a few minutes to run due to the number of sources !"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 265
    },
    "colab_type": "code",
    "id": "zNeonpiJQBQf",
    "outputId": "1d58d7bc-e62b-4f18-ac5b-7d0c1cca36a6"
   },
   "outputs": [],
   "source": [
    "N = len(ras)\n",
    "stack = 0\n",
    "c = 0 # counter to check how many thumbnails we are stacking\n",
    "\n",
    "if imap_d56 is not None:\n",
    "    stamps = reproject.thumbnails(imap_d56, coords=np.deg2rad([decs,ras]).T, res=0.0000454440, r=r)\n",
    "else:\n",
    "    stamps = reproject.thumbnails(imap, coords=np.deg2rad([decs,ras]).T, res=0.0000454440, r=r)\n",
    "stack  = np.mean(stamps,0)\n",
    "\n",
    "# Plot the Stack\n",
    "plt.figure(figsize = (7,7))\n",
    "plt.imshow(stack, extent = extent)\n",
    "plt.grid(True)\n",
    "plt.colorbar()\n",
    "plt.xlabel(\"Degrees\", fontsize = 14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It looks like the average source in our catalog is a well-defined cold spot relative to the surrounding CMB. This is qualitatively consistent with our understanding of a typical thermal-SZ cluster observed at 150 GHz (see e.g. [Carlstrom et al 2002](https://doi.org/10.1146/annurev.astro.40.060401.093803)). We examine tSZ clusters in  more detail in [Section 6](Section_6_Component_separated_maps.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='blue'> Exercise <font>\n",
    "\n",
    "The stack above uses the entire catalog but we can also stack in bins according to redshift, mass or whatever other parameters we're interested in. \n",
    "\n",
    "If you call `tab.columns` you'll be able to examine some of these different attributes (these attributes are also described in more detail [here](https://heasarc.gsfc.nasa.gov/W3Browse/ground-based/actszclust.html). Try selecting objects from the catalog based on different attributes and creating stacks of point sources that correspond to those attributes.  E.g. stacking on objects by redshift bins or mass bins. \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Applications: Aperture Photometry and Radial Binning\n",
    "Stamps are especially useful for additional analysis techniques that rely on the source being centered in the frame. One such technique is Aperture Photometry (AP), which can be used to calculate the brightness of objects. \n",
    "\n",
    "This code calculates the average temperature of an aperture: a circle of radius $x$ centered at the origin. It then finds the temperature of a ring of equal area around the circle (outer radius $\\sqrt{2}x$), and subtracts the temperature of this ring from the temperature of the aperture.  We repeat this for multiple radii and can then plot the AP result against radii.\n",
    "\n",
    "This code gives a good introduction to Aperture Photometry, however, there are a few drawbacks.  For example we don't account for distortions in the area of each pixel due to projection effects.  In order to better deal with these effects, and if you're interested in running more complex routines on these stamps we encourage you to consult the package developed by Emmanuel Schaan called [Thumbstack](https://github.com/EmmanuelSchaan/ThumbStack)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def aper_photo(x, radii_map, plot):\n",
    "    '''calculates the flux-background for a single aperature'''\n",
    "    zeros = np.zeros(np.shape(radii_map))\n",
    "    \n",
    "    circle_sum = np.sum(np.where(radii_map < x, plot, zeros))\n",
    "    circle_num = np.sum(np.where(radii_map < x, 1, zeros))\n",
    "    \n",
    "    ring_sum = np.sum(np.where(radii_map < np.sqrt(2)*x, plot, zeros)) - circle_sum\n",
    "    ring_num = np.sum(np.where(radii_map < np.sqrt(2)*x, 1, zeros)) - circle_num\n",
    "    return np.nan_to_num(circle_sum/circle_num - ring_sum/ring_num)\n",
    "\n",
    "def aper_plot(r, radii, stamp):\n",
    "    '''Creates the aperature photometry plot for a given stamp'''\n",
    "    aper_plot = np.zeros(np.shape(r))\n",
    "    for i in range(np.size(r)):\n",
    "        aper_plot[i] = aper_photo(r[i], radii, stamp)\n",
    "    return aper_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Call the function modrmap that stores the radii in degrees of each pixel from the center\n",
    "radii = enmap.modrmap(stamp.shape, stamp.wcs)*180/np.pi\n",
    "\n",
    "# Determine the bounds for the bins\n",
    "r_max = np.max(radii)\n",
    "r_min = np.min(radii)\n",
    "\n",
    "# define the radii we want for the aperature plot, note the maximum aperature radius is chosen to avoid edge effects\n",
    "r = np.linspace(r_min, r_max/2, 20)\n",
    "\n",
    "\n",
    "plt.figure(figsize = (8,7))\n",
    "plt.title(\"Aperture Photometry Plot\")\n",
    "plt.xlabel(\"Radius [degrees]\")\n",
    "plt.ylabel(\"Average Temperature [$ \\mu $K]\")\n",
    "plt.grid(True)\n",
    "plt.plot(r, aper_plot(r, radii, stamp))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The idea with aperture photometry is that the circle is a guess for the size of the object itself and the ring is the \"background\". In the plot above, the curve is about zero initially, this is because the outer ring falls entirely within the object, so the background and the circle are close to the same temperature. For greater radii the curve tilts away from zero as the temperature in the two increasingly differ and the ring \"ascends\" out of the source. Finally, we see the curve flatten off; this indicates that we've fully enclosed the source.\n",
    "\n",
    "The stamp corresponding to this AP profile is shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize = (7,7))\n",
    "plt.imshow(stamp, extent=extent,cmap = cm.coolwarm)\n",
    "plt.colorbar()\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Degrees\", fontsize = 14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='blue'> Exercise <font>\n",
    "\n",
    "Looking at the profile for a single source is always going to give you a reasonably noisy result. Try applying the AP filter technique to the full stack instead of a single point source.\n",
    "    \n",
    "You can also approximate some simple error bars for this measurement by looking at the standard deviation in AP profiles for the sources. Since selecting stamps, applying filters, etc is a bit of a slow process you might want to try this for a subset of the data instead of the full data set.\n",
    "   \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Radial Binning\n",
    "We will also take a moment to demonstrate how to radially bin a source to look at the profile. \n",
    "\n",
    "Below is code that bins data radially, fits the binned data with an exponential curve, and then creates a 2D profile from the fit to approximate the 2D profile of the source. This code can of course be adapted for other uses but here we will look at the stack of clusters, since single sources tend to be less radially symmetric. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def bin_stamp(data2d, modrmap, bin_edges):\n",
    "    '''gives the centers and averages for radial bins of a 2d plot'''\n",
    "    #get centers of each bin\n",
    "    centers = (bin_edges[1:]+bin_edges[:-1])/2.\n",
    "    \n",
    "    #sort indices by corresponding radial bin\n",
    "    digitized = np.digitize(modrmap.reshape(-1), bin_edges,right=True)\n",
    "\n",
    "    #get average of each radial bin\n",
    "    count = np.bincount(digitized)[1:-1]\n",
    "    res = np.bincount(digitized,(data2d).reshape(-1))[1:-1]/count\n",
    "    return centers,res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define the fit function\n",
    "def func(x, A, B, C):\n",
    "    return A*np.exp(B*x) + C\n",
    "\n",
    "# Set the number of bins\n",
    "N_bins = 50\n",
    "\n",
    "# r_edge is the radius of the largest circle that can fit entirely within the box\n",
    "r_edge = r_max/np.sqrt(2)\n",
    "\n",
    "step = (r_edge - r_min)/N_bins\n",
    "\n",
    "# Create the bin edges\n",
    "bin_edges = np.arange(r_min, r_edge, step)\n",
    "\n",
    "# Use the binner function to get the average temp for each radius\n",
    "cents, res = bin_stamp(stack, radii, bin_edges)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Fit the results from the binner with the exponential\n",
    "popt, pcov = curve_fit(func, cents, res, p0 = [-1, -1, 80])\n",
    "res_fit = func(cents, popt[0], popt[1], popt[2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot the radial temp curve and its fit\n",
    "fig, (ax1, ax2) = plt.subplots(2, figsize=(10,10), gridspec_kw={'height_ratios': [2, 1]})\n",
    "fig.figsize = (10, 9)\n",
    "fig.suptitle('Radial Binning with Fit', fontsize = 15)\n",
    "\n",
    "ax1.plot(cents, res, 'b.')\n",
    "ax1.plot(cents, res_fit, 'r')\n",
    "ax1.set(ylabel = r\"Temperature [$ \\mu $K]\", title = \"Average Radial Temperature\")\n",
    "ax2.plot(cents, res - res_fit, 'r.')\n",
    "\n",
    "ax2.plot(cents, np.zeros(np.size(cents)), 'k')\n",
    "ax2.set(xlabel = \"Radius [degrees]\",\n",
    "        ylabel = r\"Temperature [$ \\mu $K]\", title = \"Residuals\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As shown in the plot above the radial binning curve is nearly exponential, with relatively small residuals. This code uses a fit instead of the raw binning data to \"smooth\" the values a bit, and so that data can be extrapolated past r-edge (where the corners of the stamp would create errors in the radial binning).  As is evident in the plots above the point sources can't be fit perfectly with this curve, feel free to try other fits by adjusting the function in the previous cells!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font color='blue'> Exercise <font>\n",
    "    \n",
    "Here we've looked at the radial profile of an average SZ cluster in the CMB maps but we could do the same thing for the radial profile of an average cluster of a specific mass or at a specific redshift. \n",
    "    \n",
    "Try recreating the plots above but for the stacks you created in the first exercise where you stacked in redshift/mass/something else bins. How does the radial profile vary with the parameter you chose to look at?\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use the radial temp fit to get a 2d map\n",
    "res_fit_2D = interp1d(cents,res_fit,kind='linear',fill_value='extrapolate')(radii)\n",
    "\n",
    "# Plot the 2d fit map and the original stack with the fit subtracted\n",
    "print(\"Radially binned 2d map\")\n",
    "plt.imshow(res_fit_2D,extent=extent,cmap = cm.coolwarm)\n",
    "plt.colorbar()\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Degrees\", fontsize = 14)\n",
    "plt.show()\n",
    "\n",
    "print(\"Stack with radial binning subtracted\")\n",
    "plt.imshow(stack - res_fit_2D,extent=extent,cmap = cm.coolwarm)\n",
    "plt.colorbar()\n",
    "plt.grid(True)\n",
    "plt.xlabel(\"Degrees\", fontsize = 14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the plot above we can see that the stack is slightly asymetric, this is most likely due to random variation which the stack is not large enough to average out, or due to the inadequacy of the exponential model.  It's also worth noting we've adjusted the scale above in order to make the residuals visible."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "UqVls-28Fe6d"
   },
   "source": [
    "## Plot Messier objects\n",
    "\n",
    "While users are free to search for objects in the maps using catalogs we've also provided a tool to plot random Messier objects from our maps below.  You can enter the Messier number of an object you'd like to look at or feel free to let the function generate an object at random."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {},
    "colab_type": "code",
    "id": "DtqkPlJxGT2E"
   },
   "outputs": [],
   "source": [
    "def plot_messier(M = None):\n",
    "    # M should be an integer \n",
    "\n",
    "    # read in a table of messier objects that overlap with our feild of view\n",
    "    tbl = ascii.read(path + \"messier_objs.csv\", encoding = \"utf-8-sig\")\n",
    "\n",
    "    # Check if user has given an object to look for\n",
    "    if M != None:\n",
    "        M = str('M') + str(M)\n",
    "        A = list(tbl['M'])\n",
    "        if M in A:\n",
    "            obj = tbl[A.index(M)]\n",
    "        else:\n",
    "            print(\"Sorry that object isn't visible in this map\")\n",
    "\n",
    "    # If no object is given select a random one\n",
    "    if M == None:\n",
    "        obj = tbl[random.randrange(0, len(tbl))] \n",
    "\n",
    "    # Get the RA and DEC of the object\n",
    "    ra = coord.Angle(obj[\"RA\"], unit=u.hour) \n",
    "    dec = coord.Angle(obj[\"Declination\"], unit = u.degree) \n",
    "    \n",
    "    r = 20.*np.pi/(60 * 180)\n",
    "    \n",
    "    # Get a plot of the object\n",
    "    if imap_d56 is not None:\n",
    "        stamp = reproject.thumbnails(imap_d56 , coords = [dec.rad, ra.rad], r=r, res =  0.0000454440)\n",
    "    else:\n",
    "        stamp = reproject.thumbnails(imap , coords = [dec.rad, ra.rad],  r=r,res =  0.0000454440)\n",
    "    \n",
    "    if np.sum(stamp) == 0:\n",
    "        print(\"This object wasn't in the D56 full resolution map so we'll switch to a downgraded view.\")\n",
    "        print(\"Note this may mean it's difficult to see this particular object, consider looking at it with the full depth map instead!\")\n",
    "        stamp = reproject.thumbnails(imap_dg , coords = [dec.rad, ra.rad],  r=r,res =  0.0000454440)\n",
    "        \n",
    "\n",
    "    # Establish name for the object, if the common name is given use that, if not \n",
    "    # we just use the NGC and M indentifiers\n",
    "    if type(obj[\"Common Name\"] ) == np.str_:\n",
    "        name = str(obj[\"Common Name\"])  + str(' NGC: ') + str(obj[\"NGC\"]) + \\\n",
    "        str(' ') + str(obj[\"M\"])\n",
    "    else:\n",
    "        name = str('NGC: ') + str(obj[\"NGC\"])+ str(' ') + str(obj[\"M\"])\n",
    "\n",
    "    # Plot the object \n",
    "    plt.imshow((stamp))\n",
    "    plt.title(name)\n",
    "    plt.show()\n",
    "\n",
    "    return "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 280
    },
    "colab_type": "code",
    "id": "fi4KFCUqGVDy",
    "outputId": "4253ba53-4c71-4c7e-dd6c-7a9a3935e8b6"
   },
   "outputs": [],
   "source": [
    "# Try this a few times - some of the objects may be clearer than others\n",
    "plot_messier()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Continue to the next section\n",
    "\n",
    "[Section 5](ACTNotebook_Section_5_Lensing_maps.ipynb):  Introduction to the Lensing maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}