CMB_School_Part_02.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/jeffmcm1977/CMBAnalysis_SummerSchool/blob/master/CMB_School_Part_02.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Adding Foregrounds to the Simulated CMB\n",
    "### Jeff McMahon and Renée Hložek and Tom Crawford"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The previous notebook simulated a purely CMB realisation of the sky. We now want to add in foregrounds to the CMB map. We will start by reading in the parameters we defined in the previous notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python -c \"import cmb_modules\" || ( \\\n",
    "    wget https://github.com/jeffmcm1977/CMBAnalysis_SummerSchool/raw/master/cmb_school.tar.gz && \\\n",
    "    tar xzvf cmb_school.tar.gz \\\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "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",
    "\n",
    "import constants as cs # the constants module\n",
    "import cmb_modules # the module of functions\n",
    "N = cs.N\n",
    "c_min = cs.c_min\n",
    "c_max = cs.c_max\n",
    "X_width =cs.X_width\n",
    "Y_width = cs.Y_width\n",
    "\n",
    "pix_size = cs.pix_size"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1.3: Point Source Map\n",
    "\n",
    "Point sources in CMB maps arise from a number of astrophysical objects including active galactic nuclei (AGN), dusty star-forming galaxies (DSFGs), and a bright tail of lensed DSFGs.  All of theses objects are interesting in their own right.  We will simulate sources in several bins of brightness from the two main families.  The way we simulate the sources, the number of sources in a given pixel of our simulated map will follow a Poisson distribution---i.e., there will be no angular clustering of the sources. For now we will just make a map that looks like data maps at 150 GHz.  If we were to make multifrequency maps, we would have to take into account the different frequency scalings of these families of sources.  In general DSFGs grow brighter at higher frequencies while AGN have a spectrum that falls with increasing frequencies.\n",
    "\n",
    "Also note that in later notebooks we will use a simplified version of this code, one in which we simulate one dim population and one brighter population, without explicit source count models. The resulting maps will look qualitatively the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### define brightness bins for flux distributions. these are in units of number of sources per square degree per \n",
    "### Kelvin (CMB). normally they are defined in flux units (like per Jansky), but we convert to KCMB at 150 GHz\n",
    "### so that we can add the simulated maps directly to a CMB realization.\n",
    "\n",
    "Number_of_Bins = 100\n",
    "Brightness_Bins = np.zeros([2,Number_of_Bins])\n",
    "Minimum_Brightness = 0\n",
    "Maximum_Brightness = 5e-3 # Kelvin\n",
    "# lower edge of bins\n",
    "Brightness_Bins[0,:] = np.arange(Number_of_Bins)/Number_of_Bins*(Maximum_Brightness - Minimum_Brightness) + Minimum_Brightness\n",
    "# upper edge of bins\n",
    "Brightness_Bins[1,:] = (np.arange(Number_of_Bins) + 1)/Number_of_Bins*(Maximum_Brightness - Minimum_Brightness) + Minimum_Brightness\n",
    "\n",
    "def AGN_Counts(Brightness_Bins):\n",
    "    bins_mid = np.sum(Brightness_Bins,0)/2.\n",
    "    AGN_Counts = 5e4 * (bins_mid/1e-4)**(-1.8)\n",
    "    return(AGN_Counts)\n",
    "\n",
    "def DSFG_Counts(Brightness_Bins):\n",
    "    bins_mid = np.sum(Brightness_Bins,0)/2.\n",
    "    DSFG_Counts = 1e11*(bins_mid/1e-6)**(-2)*np.exp(-bins_mid/1e-4) + 0.5*(bins_mid/1e-2)**(-2.4)\n",
    "    return(DSFG_Counts)\n",
    "\n",
    "def make_source_map(N,pix_size,Brightness_Bins,Counts):\n",
    "    \"makes a realization of a naive Poisson-distributed point source map\"\n",
    "    \"takes in the number of pixels, the size of the pixels, the model source counts, and the brightnesses at which the model is defined\"\n",
    "    PSMap = np.zeros([int(N),int(N)])\n",
    "    area_sqdeg = (N*pix_size/60.)**2\n",
    "    for i in np.arange(len(Counts)):\n",
    "        bin = Brightness_Bins[:,i]\n",
    "        count = Counts[i]\n",
    "        bin_width = bin[1] - bin[0]\n",
    "        nsources_mean = count*bin_width*area_sqdeg\n",
    "        nsources = np.random.poisson(nsources_mean)\n",
    "        brightnesses = np.random.rand(nsources)*(bin[1] - bin[0]) + bin[0]\n",
    "        pix_x = N*np.random.rand(nsources)\n",
    "        pix_y = N*np.random.rand(nsources)\n",
    "        for j in np.arange(nsources):\n",
    "            PSMap[int(pix_x[j]),int(pix_y[j])] += brightnesses[j]\n",
    "    return(PSMap)    \n",
    "  ############################### \n",
    "    \n",
    "## make a point source map\n",
    "\n",
    "Counts1 = AGN_Counts(Brightness_Bins)\n",
    "Counts2 = DSFG_Counts(Brightness_Bins)\n",
    "PSMap1 = make_source_map(N,pix_size,Brightness_Bins,Counts1) * 1e6 # convert from K to micro-K\n",
    "PSMap2 = make_source_map(N,pix_size,Brightness_Bins,Counts2) * 1e6 # convert from K to micro-K\n",
    "PSMap = PSMap1 + PSMap2\n",
    "\n",
    "hist,bin_edges = np.histogram(PSMap,bins = 50,range=[20, 4e3])\n",
    "plt.figure(figsize=(10,10))\n",
    "plt.stairs(hist, bin_edges, fill=True)\n",
    "plt.yscale('log')\n",
    "plt.xlabel('source amplitude [$\\mu$K]')\n",
    "plt.ylabel('number or pixels')\n",
    "plt.show()\n",
    "\n",
    "p=cmb_modules.Plot_CMB_Map(PSMap,c_min,c_max,X_width,Y_width)    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top plot shows a histogram of brightness of the pixels in our source map.   The lower plot shows a map of the point source map we have simulated."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font> Plot a histogram of the AGN brightness values and a histogram of the DSFG brightness values.   How do these two distributions differ?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your discussion goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Section 1.4: SZ Map\n",
    "\n",
    "Clusters of galaxies imprint a subtle distortion into CMB maps that is most apparent on arcminute scales.    While clusters of galaxies are named after the galaxies bound within them, the galaxies represent only a small fraction of the matter contained within a cluster.   Roughly 80% of the baryons are not contained within galaxies, but rather exist as a cloud of gas bound within the gravitational potential well created by a dark matter halo that caries the vast majority of the mass of the cluster.  Within this well, the dilute gas becomes ionized and heated to temperatures of millions of Kelvin.  Occasionally a CMB photon interacts with one of the hot electrons in this ionized gas.  This interaction (inverse Compton scattering) gives the CMB photon a boost in energy.  Detailed calculations show that this effect (the Sunyev-Zeldovich or SZ effect) leads to decrement of power at frequencies below the 'null' at 220 GHz and extra power at higher frequencies.  This result is redshift independent.  Thus the SZ effect provides a clean way to detect clusters of galaxies and the signal which traces the electron density within the cluster.  The SZ signal is a reasonably good tracer of cluster mass.\n",
    "\n",
    "For these simulations we treat each cluster as having a brightness \"beta profile\", and fix each cluster to have an identical angular size.  We draw the distribution of central temperatures from an exponential distribution to simplify the code and reduce the dependance on external libraries.  For more accurate simulations, a range of clusters sizes should be used, a distribution of cluster shapes (with more accurate profiles) should be considered, and the number of clusters as a function of mass and redshift should be chosen to match measurements of the cluster mass function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "### paramaters to set up the SZ point sources\n",
    "Number_of_SZ_Clusters  = 500\n",
    "Mean_Amplitude_of_SZ_Clusters = 50\n",
    "SZ_beta = 0.86\n",
    "SZ_Theta_core = 1.0\n",
    "\n",
    "def SZ_source_component(N,pix_size,Number_of_SZ_Clusters,Mean_Amplitude_of_SZ_Clusters,SZ_beta,SZ_Theta_core,do_plots):\n",
    "    \"makes a realization of a naive SZ map\"\n",
    "    N=int(N)\n",
    "    SZMap = np.zeros([N,N])\n",
    "    SZcat = np.zeros([3,Number_of_SZ_Clusters]) ## catalogue of SZ sources, X, Y, amplitude\n",
    "    # make a distribution of point sources with varying amplitude\n",
    "    i = 0\n",
    "    while (i < Number_of_SZ_Clusters):\n",
    "        pix_x = int(N*np.random.rand())\n",
    "        pix_y = int(N*np.random.rand())\n",
    "        pix_amplitude = np.random.exponential(Mean_Amplitude_of_SZ_Clusters)*(-1.)\n",
    "        SZcat[0,i] = pix_x\n",
    "        SZcat[1,i] = pix_y\n",
    "        SZcat[2,i] = pix_amplitude\n",
    "        SZMap[pix_x,pix_y] += pix_amplitude\n",
    "        i = i + 1\n",
    "    if (do_plots):\n",
    "        hist,bin_edges = np.histogram(SZMap,bins = 50,range=[SZMap.min(),-10])\n",
    "        plt.figure(figsize=(10,10))\n",
    "        plt.semilogy(bin_edges[0:-1],hist)\n",
    "        plt.xlabel('source amplitude [$\\mu$K]')\n",
    "        plt.ylabel('number of pixels')\n",
    "        plt.show()\n",
    "    \n",
    "    # make a beta function\n",
    "    beta = beta_function(int(N),pix_size,SZ_beta,SZ_Theta_core)\n",
    "    \n",
    "    # convolve the beta function with the point source amplitude to get the SZ map\n",
    "    # NOTE: you should go back to the Intro workshop for more practice with convolutions!\n",
    "    FT_beta = np.fft.fft2(np.fft.fftshift(beta))\n",
    "    FT_SZMap = np.fft.fft2(np.fft.fftshift(SZMap))\n",
    "    SZMap = np.fft.fftshift(np.real(np.fft.ifft2(FT_beta*FT_SZMap)))\n",
    "    \n",
    "    # return the SZ map\n",
    "    return(SZMap,SZcat)    \n",
    "  ############################### \n",
    "\n",
    "def beta_function(N,pix_size,SZ_beta,SZ_Theta_core):\n",
    "  # make a beta function\n",
    "\n",
    "    N=int(N)\n",
    "    ones = np.ones(N)\n",
    "    inds  = (np.arange(N)+.5 - N/2.) * pix_size\n",
    "    X = np.outer(ones,inds)\n",
    "    Y = np.transpose(X)\n",
    "    # compute the same real-space R function as before for the PS\n",
    "    R = np.sqrt(X**2. + Y**2.)\n",
    "    \n",
    "    beta = (1 + (R/SZ_Theta_core)**2.)**((1-3.*SZ_beta)/2.)\n",
    "\n",
    "    # return the beta function map\n",
    "    return(beta)\n",
    "  ############################### \n",
    "    \n",
    "## make a point source map\n",
    "\n",
    "SZMap,SZCat = SZ_source_component(N,pix_size,Number_of_SZ_Clusters,Mean_Amplitude_of_SZ_Clusters,SZ_beta,SZ_Theta_core,True)\n",
    "\n",
    "  \n",
    "p=cmb_modules.Plot_CMB_Map(SZMap,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The top plot shows the a histogram of the SZ-decrements from our simulated SZ cluster map.  The bottom plot shows our simulated SZ map. This map is at 150 GHz. How could you tell just from looking at the code above?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  In reality the cluster radii vary from cluster to cluster.  For reference the radius depends on redshift and mass.  The number of clusters as a function of mass and redshift is called the cluster 'mass function' which is a sensitive cosmological proble.  \n",
    "\n",
    "To enhance the realism of our sims, divide the simulated cluster sample into an extremely large radius sample (1 cluster with a 30 arcminute radius, comparable to the Coma cluster; the largest cluster on the sky), a large radius bin (10%) clusters with 5 arcminute radius), a medium bin (30%) with 2 arcminute radius, and a small bin (60% with 0.5 arcminute radius).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# your code goes here\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comment on the differences between this map and the original version.  And can you think of a reason why (given the algorithm we are using) we didn't draw the cluster radius from a continuous distribution?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.5 Full Sky Map\n",
    "\n",
    "The sky map is a combination of the CMB anisotropy, a point source map, and an SZ map.  In an appendix we add the impact of CMB lensing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## remake the CMB map from the pervious excercise\n",
    "ell, DlTT = np.loadtxt(\"CAMB_fiducial_cosmo_scalCls.dat\", usecols=(0, 1), unpack=True) \n",
    "CMB_T = cmb_modules.make_CMB_T_map(N,pix_size,ell,DlTT)\n",
    "\n",
    "## add all the components together\n",
    "total_map = CMB_T + PSMap + SZMap\n",
    "\n",
    "p=cmb_modules.Plot_CMB_Map(total_map,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows our simulated map that includes CMB, point source, and SZ cluster signals.  Note that the sources seem brighter than what we saw in the real observed maps.   This is not a mistake, as will be seen after we fold in the beam (point spread function) of the instrument."
   ]
  }
 ],
 "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}