CMB_School_Part_03.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a href=\"https://colab.research.google.com/github/jeffmcm1977/CMBAnalysis_SummerSchool/blob/master/CMB_School_Part_03.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Instrument Beam, Instrument Noise, and Filtering\n",
    "### Jeff McMahon and Renée Hložek "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will now add noise and filter the map - so that we have a map with all astrophysical and instrumental effects."
   ]
  },
  {
   "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",
    "N = int(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",
    "beam_size_fwhp = cs.beam_size_fwhp\n",
    "\n",
    "pix_size = cs.pix_size\n",
    "\n",
    "Number_of_Sources  = cs.Number_of_Sources\n",
    "Amplitude_of_Sources = cs.Amplitude_of_Sources\n",
    "Number_of_Sources_EX = cs.Number_of_Sources_EX\n",
    "Amplitude_of_Sources_EX = cs.Amplitude_of_Sources_EX\n",
    "\n",
    "Number_of_SZ_Clusters  = cs.Number_of_SZ_Clusters\n",
    "Mean_Amplitude_of_SZ_Clusters = cs.Mean_Amplitude_of_SZ_Clusters\n",
    "SZ_beta = cs.SZ_beta\n",
    "SZ_Theta_core = cs.SZ_Theta_core"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Recalculating the results from the previous stages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## Make a CMB map\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",
    "## make a point source map\n",
    "PSMap = cmb_modules.Poisson_source_component(N,pix_size,Number_of_Sources,Amplitude_of_Sources) \n",
    "PSMap +=  cmb_modules.Exponential_source_component(N,pix_size,Number_of_Sources_EX,Amplitude_of_Sources_EX)\n",
    "\n",
    "## make an SZ map\n",
    "SZMap,SZCat = cmb_modules.SZ_source_component(N,pix_size,Number_of_SZ_Clusters,\\\n",
    "                                  Mean_Amplitude_of_SZ_Clusters,SZ_beta,SZ_Theta_core,False)\n",
    "\n",
    "## add them all together to get the sky map at a single freuqency\n",
    "total_map = CMB_T + PSMap + SZMap\n",
    "\n",
    "## plot the result\n",
    "p=cmb_modules.Plot_CMB_Map(total_map,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Sky Map Convolved with a Beam\n",
    "\n",
    "Telescopes suffer from diffraction which leads to finite resolution effects.  To account for this we generate a Gaussian beam pattern and convolve the map with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "beam_size_fwhp = 1.25 # arcminutes\n",
    "\n",
    "def convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,Map):\n",
    "    \"convolves a map with a Gaussian beam pattern.  NOTE: pix_size and beam_size_fwhp need to be in the same units\" \n",
    "    # make a 2d gaussian \n",
    "    gaussian = make_2d_gaussian_beam(N,pix_size,beam_size_fwhp)\n",
    "  \n",
    "    # do the convolution\n",
    "    FT_gaussian = np.fft.fft2(np.fft.fftshift(gaussian)) # first add the shift so that it is central\n",
    "    FT_Map = np.fft.fft2(np.fft.fftshift(Map)) #shift the map too\n",
    "    convolved_map = np.fft.fftshift(np.real(np.fft.ifft2(FT_gaussian*FT_Map))) \n",
    "    \n",
    "    # return the convolved map\n",
    "    return(convolved_map)\n",
    "  ###############################   \n",
    "\n",
    "def make_2d_gaussian_beam(N,pix_size,beam_size_fwhp):\n",
    "     # make a 2d coordinate system\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",
    "    R = np.sqrt(X**2. + Y**2.)\n",
    "    # plt.title('Radial coordinates')\n",
    "    # plt.imshow(R)\n",
    "    # plt.colorbar()\n",
    "    # plt.show()\n",
    "  \n",
    "    # make a 2d gaussian \n",
    "    beam_sigma = beam_size_fwhp / np.sqrt(8.*np.log(2))\n",
    "    gaussian = np.exp(-.5 *(R/beam_sigma)**2.)\n",
    "    gaussian = gaussian / np.sum(gaussian)\n",
    "    # return the gaussian\n",
    "    # plt.title('Gaussian beam')\n",
    "    # plt.imshow(gaussian)\n",
    "    # plt.colorbar()\n",
    "    # plt.show()\n",
    "    return(gaussian)\n",
    "  ###############################  \n",
    "    \n",
    "    \n",
    "## convolve the signal part of the map\n",
    "\n",
    "CMB_T_convolved = convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,total_map)\n",
    "\n",
    "p=cmb_modules.Plot_CMB_Map(CMB_T_convolved,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows the simulated map that includes the CMB, point source, and SZ signals and is convolved with an instrumental beam.  The agreement with the real sky maps is now reasonable, but does not yet include instrumental noise. \n",
    "\n",
    "If you want to see a bit more of the inner-workings of what happened, uncomment the lines in the make_2d_gaussian_beam function that make plots of the R coordinates and beam itself."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  There are three typical beam sizes for CMB experiments: 1) large angular scale B-mode observatories which have ~30 arcminute beams; 2) medium scale observatories such as the Planck satellite that has a ~3 arcminute beam, and 3) high resolution observatories such as ACT and SPT that have ~1 arcminute beams.  Convolve your map with each of these cases and compare the differences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "comment on the differences here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noise Map\n",
    "\n",
    "Ground based CMB instruments suffer from several types of noise.  This includes: (1) White noise drawn from a gaussian distribution; (2) atmospheric noise which grows larger on large angular scales, this can be understood in terms of Kolmogorov turbulence; and (3) $1/f$ noise in the detectors.  Both atmospheric and detector $1/f$ manifest as correlated noise in map space.  However, like the CMB these noise components are very nearly uncorrelated in Fourier space.  We will build these noise terms in Fourier space and display them in map space. \n",
    "\n",
    "We choose a white noise level of 10 $\\mu\\mathrm{K}$-arcmin to approximate the current deep maps coming from ACTPol and SPTPol.  The $1/f$ noise is similar to what would be seen in an experiment at the south pole.  When turned on, this noise leads to the effect of \"striping\" in the maps.  In later notebooks, we default to turning the atmospheric and 1/f noise components off since if they are included we may want to filter the maps before estimating the power spectrum.  (This is easy, but is not yet implemented).  The atmospheric noise (Kolmogorov turbulence) has a 2D spectrum similar to CMB leading to similar patterns in the maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "white_noise_level = 10.\n",
    "atmospheric_noise_level = 0.5  # multiply by zero to turn this off\n",
    "one_over_f_noise_level = 0.2 # multiply by zero to turn this off\n",
    "\n",
    "def make_noise_map(N,pix_size,white_noise_level,atmospheric_noise_level,one_over_f_noise_level):\n",
    "    \"makes a realization of instrument noise, atmosphere and 1/f noise level set at 1 degrees\"\n",
    "    ## make a white noise map\n",
    "    white_noise = np.random.normal(0,1,(N,N)) * white_noise_level/pix_size\n",
    " \n",
    "    ## make an atnosperhic noise map\n",
    "    atmospheric_noise = 0.\n",
    "    if (atmospheric_noise_level != 0):\n",
    "        ones = np.ones(N)\n",
    "        inds  = (np.arange(N)+.5 - N/2.) \n",
    "        X = np.outer(ones,inds)\n",
    "        Y = np.transpose(X)\n",
    "        R = np.sqrt(X**2. + Y**2.) * pix_size /60. ## angles relative to 1 degrees  \n",
    "        mag_k = 2 * np.pi/(R+.01)  ## 0.01 is a regularization factor\n",
    "        atmospheric_noise = np.fft.fft2(np.random.normal(0,1,(N,N)))\n",
    "        atmospheric_noise  = np.fft.ifft2(atmospheric_noise * np.fft.fftshift(mag_k**(5/3.)))\n",
    "        atmospheric_noise = atmospheric_noise * atmospheric_noise_level/pix_size\n",
    "\n",
    "    ## make a 1/f map, along a single direction to illustrate striping \n",
    "    oneoverf_noise = 0.\n",
    "    if (one_over_f_noise_level != 0): \n",
    "        ones = np.ones(N)\n",
    "        inds  = (np.arange(N)+.5 - N/2.) \n",
    "        X = np.outer(ones,inds) * pix_size /60. ## angles relative to 1 degrees \n",
    "        kx = 2 * np.pi/(X+.01) ## 0.01 is a regularization factor\n",
    "        oneoverf_noise = np.fft.fft2(np.random.normal(0,1,(N,N)))\n",
    "        oneoverf_noise = np.fft.ifft2(oneoverf_noise * np.fft.fftshift(kx))\n",
    "        oneoverf_noise = oneoverf_noise * one_over_f_noise_level/pix_size\n",
    "\n",
    "    ## return the noise map\n",
    "    noise_map = np.real(white_noise + atmospheric_noise + oneoverf_noise)\n",
    "    return(noise_map)\n",
    "  ###############################\n",
    "\n",
    "## make an instrument noise map\n",
    "Noise = make_noise_map(N,pix_size,white_noise_level,atmospheric_noise_level,one_over_f_noise_level)\n",
    "\n",
    "p=cmb_modules.Plot_CMB_Map(Noise,c_min/2.,c_max/2.,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This plot shows a realization of instrumental + atmospheric noise.  The atmospheric noise looks similar to the CMB since it has a similar 2d spectrum.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  The above simulation has both atmospheric and detector $1/f$ noise.   To understand these components seperately, turn each off and make plots of each component seperately.  Comment on how these look.  \n",
    "\n",
    "Optional Part 2: Plot the absolute value of the 2-d FFT of these maps and note how the non-white component of the noise is localized.   \n",
    "\n",
    "NOTE: This map assumes no cross-linking, the $1/f$ would look different if there were two or more scan directions in the map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code and plots go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Your comments go here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Filtering\n",
    "\n",
    "In maps with frequency dependent noise it is often advantageous to filter the maps to 'whiten' the noise spectrum.   There are a wide variety of choices of what to do here.  For both atmospheric and $1/f$ noise, a high pass filter would be an effective choice.  This could either be a square filter (where the modes below some cutoff are set to zero while shorter wavelength modes are kept unchanged) or a filter that weights the fourier modes by $1/A$ where $A$ is the amplitude of the noise distribution.    Below is code that masks out $N_{mask}$ fourier modes along the scan direction.  See the exercise below and play with it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "N_mask =4\n",
    "\n",
    "def Filter_Map(Map,N,N_mask):\n",
    "    ## set up a x, y, and r coordinates for mask generation\n",
    "    ones = np.ones(N) # N is the size of the map\n",
    "    inds  = (np.arange(N)+.5 - N/2.) \n",
    "\n",
    "    X = np.outer(ones,inds)\n",
    "    Y = np.transpose(X)\n",
    "    R = np.sqrt(X**2. + Y**2.)  ## angles realative to 1 degrees  \n",
    "    \n",
    "    ## make a mask\n",
    "    mask  = np.ones([N,N])\n",
    "    mask[np.where(np.abs(X) < N_mask)]  = 0\n",
    "    \n",
    "    ## apply the filter in fourier space\n",
    "    FMap = np.fft.fftshift(np.fft.ifft2(np.fft.fftshift(Map)))\n",
    "    FMap_filtered = FMap * mask\n",
    "    Map_filtered = np.real(np.fft.fftshift(np.fft.fft2(np.fft.fftshift(FMap_filtered))))\n",
    "    \n",
    "    ## return the output\n",
    "    return(Map_filtered)\n",
    "\n",
    "\n",
    "filtered = Filter_Map(Noise,N,N_mask)\n",
    "p=cmb_modules.Plot_CMB_Map(filtered,c_min/4.,c_max/4.,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<font color='red'>EXERCISE: </font>  Find the minimum $N_{mask}$ that effectively whitens the noise spectrum with the bad $1/f$ turned on.  See how this filter works for atmospheric noise.  Comment on whether this filter works and whether it is a good choice for atmospheric noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "## your code goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "your notes go here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Complete Simulated CMB Map\n",
    "\n",
    "We complete our simulated CMB map by adding the simulated sky map convolved with the beam to our simulated noise map.   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "CMB_T_convolved =convolve_map_with_gaussian_beam(N,pix_size,beam_size_fwhp,total_map)\n",
    "\n",
    "total_map_plus_noise = CMB_T_convolved + Noise\n",
    "\n",
    "p=cmb_modules.Plot_CMB_Map(total_map_plus_noise,c_min,c_max,X_width,Y_width)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "source": [
    "This plot shows a complete simulated CMB map including the key astrophysical and instrumental effects. Compare it to the CMB map without noise at the top of this notebook and note the change!"
   ]
  }
 ],
 "metadata": {
  "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
}