CMB_School_Part_11_comp_sep_SPT_Planck.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Component Separation With SPT Data\n",
    "\n",
    "### Tom Crawford"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What will we be doing in this notebook?\n",
    "\n",
    "In the previous notebook, you created independent CMB and tSZ component maps from simulated three-band data. The simulated noise was uncorrelated between pixels and uniform across the maps, so you could use a single covariance matrix for the entire map, and that matrix was diagonal. Also, each band had the same angular resolution (beam).\n",
    "\n",
    "Both ACT (Madhavacheril et al. 2020) and SPT (Bleem et al. 2021) have recently published and released the results of component-separation analyses, in which one of the main results were independent CMB and tSZ maps. In both of these analyses, the high-resolution ground-based data (ACT or SPT) was combined with Planck data. In this notebook, we will do a scaled-down version of the SPT+Planck component separation. We will still only compute a single $N_\\mathrm{band}$-by-$N_\\mathrm{band}$ covariance matrix. This is very much not the full description of the noise (including astrophysical signals we are not interested in) in the SPT+Planck maps. Luckily, this only makes the component maps less optimal; it does not make them biased.* We do have to deal with the fact that the different frequencies have different beams, though, otherwise the component maps will look very weird. We will do the easiest possible thing and smooth all the maps to roughly the same resolution.\n",
    "\n",
    "(* This is actually a fairly subtle feature of linear-least squares solutions and applications like matched filtering and component separation. Think about why this might be true...)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## How to get the SPT data?\n",
    "\n",
    "The single-frequency SPT maps that went into the Bleem et al. analysis are available at the NASA/LAMBDA site (https://lambda.gsfc.nasa.gov/) and from Lindsey Bleem's box server at Argonne (https://anl.app.box.com/s/5yr7z0xrgkebpzukcjvw7ve8hr4dwony). But we are going to cheat and use some maps that have already been combined with Planck data at the nearest frequency to each SPT band. We will work with one small subfield (centered at RA = 6h30m, declination = -55d). The files for this are on my personal KICP page: https://kicp.uchicago.edu/~tcrawfor/summer_school_data/.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import the usual stuff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib\n",
    "import sys\n",
    "import matplotlib.cm as cm\n",
    "import matplotlib.mlab as mlab\n",
    "import matplotlib.pyplot as plt\n",
    "import astropy.io.fits as fits\n",
    "from astropy.coordinates import SkyCoord\n",
    "from astropy.wcs import WCS\n",
    "\n",
    "%matplotlib inline\n",
    "import constants as cs # the constants module\n",
    "\n",
    "from cmb_modules import * # the module of functions\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Get SPT (+Planck) maps"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "ename": "IOError",
     "evalue": "[Errno 2] No such file or directory: 'map_sptsz90_planck100_ra6h30dec-55.fits'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mIOError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-3-f3132b9b73fb>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m#spt_files = ['sptsz_ymaprun_map_ra6h30dec-55_' + bn + 'GHz.fits' for bn in spt_band_names]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      6\u001b[0m \u001b[0mspt_files\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0;34m'map_'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mfile_stub\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;34m'_ra6h30dec-55.fits'\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mfile_stub\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mfile_stubs\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 7\u001b[0;31m \u001b[0mhdulist_temp\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfits\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mspt_files\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      8\u001b[0m \u001b[0mheader\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhdulist_temp\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mheader\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      9\u001b[0m \u001b[0mnx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mheader\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'NAXIS1'\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/cvmfs/spt.opensciencegrid.org/py2-v1/RHEL_7_x86_64/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc\u001b[0m in \u001b[0;36mfitsopen\u001b[0;34m(name, mode, memmap, save_backup, cache, **kwargs)\u001b[0m\n\u001b[1;32m    136\u001b[0m         \u001b[0;32mraise\u001b[0m \u001b[0mValueError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Empty filename: %s'\u001b[0m \u001b[0;34m%\u001b[0m \u001b[0mrepr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    137\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 138\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mHDUList\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfromfile\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmemmap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msave_backup\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    139\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    140\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/cvmfs/spt.opensciencegrid.org/py2-v1/RHEL_7_x86_64/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc\u001b[0m in \u001b[0;36mfromfile\u001b[0;34m(cls, fileobj, mode, memmap, save_backup, cache, **kwargs)\u001b[0m\n\u001b[1;32m    278\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    279\u001b[0m         return cls._readfrom(fileobj=fileobj, mode=mode, memmap=memmap,\n\u001b[0;32m--> 280\u001b[0;31m                              save_backup=save_backup, cache=cache, **kwargs)\n\u001b[0m\u001b[1;32m    281\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    282\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mclassmethod\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/cvmfs/spt.opensciencegrid.org/py2-v1/RHEL_7_x86_64/lib/python2.7/site-packages/astropy/io/fits/hdu/hdulist.pyc\u001b[0m in \u001b[0;36m_readfrom\u001b[0;34m(cls, fileobj, data, mode, memmap, save_backup, cache, **kwargs)\u001b[0m\n\u001b[1;32m    799\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_File\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    800\u001b[0m                 \u001b[0;31m# instantiate a FITS file object (ffo)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 801\u001b[0;31m                 \u001b[0mffo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_File\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmemmap\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mmemmap\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mcache\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mcache\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    802\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    803\u001b[0m                 \u001b[0mffo\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfileobj\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/cvmfs/spt.opensciencegrid.org/py2-v1/RHEL_7_x86_64/lib/python2.7/site-packages/astropy/io/fits/file.pyc\u001b[0m in \u001b[0;36m__init__\u001b[0;34m(self, fileobj, mode, memmap, clobber, cache)\u001b[0m\n\u001b[1;32m    139\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_open_fileobj\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mclobber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    140\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstring_types\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 141\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_open_filename\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mclobber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    142\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    143\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_open_filelike\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfileobj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mclobber\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/cvmfs/spt.opensciencegrid.org/py2-v1/RHEL_7_x86_64/lib/python2.7/site-packages/astropy/io/fits/file.pyc\u001b[0m in \u001b[0;36m_open_filename\u001b[0;34m(self, filename, mode, clobber)\u001b[0m\n\u001b[1;32m    491\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_file\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbz2\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mBZ2File\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mbzip2_mode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    492\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 493\u001b[0;31m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_file\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfileobj_open\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mPYFITS_MODES\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mmode\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    494\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    495\u001b[0m         \u001b[0;31m# Make certain we're back at the beginning of the file\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m/cvmfs/spt.opensciencegrid.org/py2-v1/RHEL_7_x86_64/lib/python2.7/site-packages/astropy/io/fits/util.pyc\u001b[0m in \u001b[0;36mfileobj_open\u001b[0;34m(filename, mode)\u001b[0m\n\u001b[1;32m    413\u001b[0m         \"\"\"\n\u001b[1;32m    414\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 415\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfilename\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmode\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    416\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    417\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mIOError\u001b[0m: [Errno 2] No such file or directory: 'map_sptsz90_planck100_ra6h30dec-55.fits'"
     ]
    }
   ],
   "source": [
    "spt_band_names = ['90','150','220']\n",
    "file_stubs = ['sptsz90_planck100','sptsz150_planck143','sptsz220_planck217']\n",
    "spt_band_centers = np.array([95.,152.,219.])\n",
    "nb = len(spt_band_names)\n",
    "#spt_files = ['sptsz_ymaprun_map_ra6h30dec-55_' + bn + 'GHz.fits' for bn in spt_band_names]\n",
    "spt_files = ['map_' + file_stub + '_ra6h30dec-55.fits' for file_stub in file_stubs]\n",
    "hdulist_temp = fits.open(spt_files[0])\n",
    "header = hdulist_temp[0].header\n",
    "nx = header['NAXIS1']\n",
    "ny = header['NAXIS2']\n",
    "w = WCS(header,naxis=2)\n",
    "spt_maps = np.zeros([nb,ny,nx])\n",
    "spt_noise_maps = np.zeros([nb,ny,nx])\n",
    "i = 0\n",
    "c_min = -6e-4\n",
    "c_max = 6e-4\n",
    "X_width = 3000\n",
    "Y_width = 3000\n",
    "for spt_file in spt_files:\n",
    "    hdulist_temp = fits.open(spt_file)\n",
    "    maps_temp = hdulist_temp[0].data\n",
    "    spt_maps[i,:,:] = maps_temp\n",
    "    hdulist_temp = fits.open(spt_file.replace('map','dmap'))\n",
    "    maps_temp = hdulist_temp[0].data\n",
    "    spt_noise_maps[i,:,:] = maps_temp\n",
    "    # plot the result\n",
    "    plt.title(spt_band_names[i]+' GHz, signal map',fontsize=18)\n",
    "    p=Plot_CMB_Map(spt_maps[i,:,:],c_min,c_max,X_width,Y_width)\n",
    "    plt.title(spt_band_names[i]+' GHz, noise map',fontsize=18)\n",
    "    p=Plot_CMB_Map(spt_noise_maps[i,:,:],c_min,c_max,X_width,Y_width)\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Smooth SPT maps to a common resolution\n",
    "\n",
    "SPT beams are shaped roughly like Gaussians, and the full width at half maximum (FWHM) of those Gaussians are roughly 1.7, 1.2, and 1.1 at 90, 150, and 220 GHz, respectively. Let's smooth the maps so that they all have a final resolution of roughly 2 arcminutes. How would you do that?\n",
    "\n",
    "#### Three hints\n",
    "\n",
    "1. A Gaussian convolved with another Gaussian has a FWHM equal to the quadrature sum (square root of the sum of the squares) of the two individual Gaussians.\n",
    "2. The \"sigma\" of a Gaussian is equal to the FWHM divided by $\\sqrt{8 \\times \\ln{2}}$, or about 2.35.\n",
    "3. The scipy function ndimage.gaussian_filter(map,sigma) smooths an array by a Gaussian with the input sigma (IN UNITS OF ARRAY INDICES OR MAP PIXELS)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# smooth to a common resolution---but only the signal maps\n",
    "from scipy import ndimage\n",
    "destination_fwhm = 2.0 # arcmin\n",
    "spt_fwhm = np.array([1.7,1.2,1.1]) # arcmin\n",
    "# your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## noise levels in the SPT maps\n",
    "\n",
    "There are two types of maps in the SPT files: \"signal\" maps and \"noise\" maps. The noise maps are created by making a map from half the data, making another map from half the data, and differencing the two. Recall that the white noise levels in the sim notebook were defined in units of microKelvin-arcminute (or uK-arcmin). Using the SPT noise maps (knowing that the SPT maps are in units of Kelvin), how would you measure the noise in units of uK-arcmin?\n",
    "\n",
    "#### Warning: \n",
    "\n",
    "You can see from the plots above that the maps don't fill the entire array. Maybe only use a cutout around the center of the noise map to measure the noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "# measure SPT noise levels\n",
    "spt_white_noise_levels = np.zeros(nb)\n",
    "# your code goes here\n",
    "\n",
    "print(spt_white_noise_levels)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Constructing the band weights\n",
    "\n",
    "Do just what you did for the sims, but with the SPT data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now multiply the SPT frequency maps by the band weights."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'ny' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mNameError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-6-eea544e9340c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mspt_component_maps\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mzeros\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mny\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mnx\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# your code goes here\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m \u001b[0;31m# plot the result\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'ny' is not defined"
     ]
    }
   ],
   "source": [
    "spt_component_maps = np.zeros([2,ny,nx])\n",
    "\n",
    "# your code goes here\n",
    "\n",
    "# plot the result\n",
    "compnames = ['CMB','tSZ']\n",
    "c_mins = [-6e-4,-2e-4]\n",
    "c_maxs = [6e-4,2e-4]\n",
    "for i in np.arange(2):\n",
    "    plt.title('reconstructed ' + compnames[i],fontsize=18)\n",
    "    p=Plot_CMB_Map(spt_component_maps[i,:,:]-np.mean(spt_component_maps[i,:,:]),c_mins[i],c_maxs[i],X_width,Y_width)#,Title=str(bands[i])+' GHz')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise: Find a cluster!\n",
    "\n",
    "The tSZ component map is in units of \"SZ strength\" or---depending on exactly how you defined the tSZ SED, maybe even units of Compton-y parameter---so massive galaxy clusters will show up as bright spots. Do you see a very bright spot along the left edge of the map, a bit less than halfway up? Figure out what the x,y location of that bright spot is (either by zooming in and eyeballing or by finding a local maximum of the tSZ map). Use the astropy WCS utilities and the WCS information in the FITS header of the SPT map file to convert that x,y location to an RA,dec location, then search one of the astronomical databases (like NED: http://ned.ipac.caltech.edu/ for that location. Is there a famous galaxy cluster nearby?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'ra' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m\u001b[0m",
      "\u001b[0;31mNameError\u001b[0mTraceback (most recent call last)",
      "\u001b[0;32m<ipython-input-7-91a3f227888c>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;31m# your code goes here\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      4\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 5\u001b[0;31m \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'RA: '\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mra\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m' deg'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      6\u001b[0m \u001b[0;32mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'dec: '\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdec\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;34m' deg'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mNameError\u001b[0m: name 'ra' is not defined"
     ]
    }
   ],
   "source": [
    "from astropy import wcs\n",
    "\n",
    "# your code goes here\n",
    "\n",
    "print('RA: '+str(ra)+' deg')\n",
    "print('dec: '+str(dec)+' deg')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Advanced question:\n",
    "\n",
    "What's that wispy, filamentary stuff in the tSZ map?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}