COBE + FIRAS CMB Spectrum.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "In this section of the lab, we'll use some data from NASA to determine the temperature of the cosmic microwave background (CMB).\n",
    "\n",
    "The CMB is the radiation left over from the big bang, and it exhibits a nearly perfect black body distribution. Deviation for a perfectly uniform distibution as a function of direction (there are very small temperature fluctuations as you look in different directions) are of intense interest in cosmology.\n",
    "\n",
    "We'll be fitting the monopole/uniform background data, which we have downloaded here\n",
    "\n",
    "Data source: https://lambda.gsfc.nasa.gov/product/cobe/firas_products.cfm\n",
    "\n",
    "If you'd like to read more about COBE and FIRAS, here are some great starting points:\n",
    "\n",
    "https://en.wikipedia.org/wiki/Cosmic_Background_Explorer\n",
    "\n",
    "https://lambda.gsfc.nasa.gov/product/cobe/firas_overview.cfm\n",
    "\n",
    "Execute the cells below to import some required libraries"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "import numpy as np\n",
    "from scipy.optimize import curve_fit\n",
    "from matplotlib import pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "In class we discussed the intensity distribution of a black body radiator. To use the NASA data, we'll need to instead work with spectral irradiance per solid angle $B$.\n",
    "\n",
    "$$ B(f, T) = \\frac{2 h f^3}{c^2}\\frac{1}{exp\\left(\\frac{hf}{k_BT}\\right) -1}$$\n",
    "\n",
    "It describes the intensity per frequency per solid angle at a given temperature. \n",
    "\n",
    "\n",
    "Solid angle wikipedia article: https://en.wikipedia.org/wiki/Solid_angle\n",
    "\n",
    "**Question 1** \n",
    "What are the SI units of $B$?\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Next, define a function to evaluate $B(f, t)$, along with some constants we need"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "outputs": [],
   "source": [
    "h = 6.62606e-34 # m^2 kg / s      Planck's constant\n",
    "c = 2.998e8 # m / s               Speed of light\n",
    "kB = 1.38065e-23 # J / k\n",
    "\n",
    "def bb_spectrum_freq(freqs, T):\n",
    "    return 2 * h * np.power(freqs, 3) / ( c * c * (np.exp(h*freqs/(kB * T)) - 1))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Next, we'll load the data. It contains frequencies $f$ (`freqs`) at which $B$ (`bs`) was measured, along with uncertainties in $B$, $\\sigma$ (`sigmas`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "outputs": [],
   "source": [
    "def load_firas_data():\n",
    "    freqs = list()\n",
    "    bs = list()\n",
    "    sigmas = list()\n",
    "    Jy_conversion = 1e-20 # 1MJy to W/m^2\n",
    "\n",
    "    with open('data/firas_monopole_spec_v1.txt', 'rt') as data_file:\n",
    "        for line in data_file:\n",
    "            if line.startswith('#'):\n",
    "                continue\n",
    "            freq, intensity, resid, sigma,  __ = map(float, line.split())\n",
    "            freqs.append(freq)\n",
    "            bs.append(intensity + resid * 1e-3)\n",
    "            sigmas.append(sigma)\n",
    "    sigmas = np.array(sigmas) * Jy_conversion * 1e-3\n",
    "    freqs = np.array(freqs) *  c * 100 # convert to hertz (from cm-1)\n",
    "    bs = np.array(bs) * (Jy_conversion)\n",
    "    \n",
    "    return freqs, bs, sigmas\n",
    "\n",
    "freqs, bs, sigmas = load_firas_data()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Now let's plot it!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.errorbar(freqs, bs, sigmas, fmt='.')\n",
    "plt.ylabel('Irradiance ($W m^{-2} Hz^{-1}str^{-1}$)')\n",
    "plt.xlabel('$f (Hz)$')\n",
    "plt.title('COBE FIRAS Cosmic Microwave Background')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Wien's displacement Law for frequencies states that\n",
    "\n",
    "$$f_{peak}/T = 5.879 \\times 10^{10} Hz / K $$\n",
    "\n",
    "\n",
    "**Question 2**\n",
    "\n",
    "What is value of f_peak (find from the graph above)? What temperature does this predict for the CMB?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f_peak = # ENTER YOUR VALUE\n",
    "T_peak =  f_peak / 5.879e10 # compute T\n",
    "\n",
    "print('T = {}'.format(T_peak))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Next, we'll plot this distribution versus our data. The function `bb_spectrum_freq(frequencies, T)` takes an array of frequencies (in Hz) as its first argument and a temperature (in Kelvin). It returns the values of the black body irradiance spectrum evaluated at the array of frequencies.\n",
    "\n",
    "\n",
    "We'll use a few functions below to generate our plot.\n",
    "\n",
    "`np.linspace(min_value, max_value, number_samples)` will return an array of evenly spaced values between min_value and max_value.\n",
    "\n",
    "`plt.plot(x, y, label='label')` will plot x vs y in the current figure.\n",
    "\n",
    "`plt.errorbar(x, y, yerr, fmt='.')` will plot x vs y and add error bars of width `yerr` to the plot\n",
    "\n",
    "Plot the theoretical data at the temperature you computed above against the FIRE data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    },
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "plt.figure(figsize=(8, 9))\n",
    "f_theory = np.linspace(5e9, 8e11, 1000)\n",
    "\n",
    "bb_theory = bb_spectrum_freq(f_theory, T=T_peak)\n",
    "\n",
    "plt.plot(f_theory, bb_theory, '-k', label='Theory {:0.3f}K'.format(T_peak))\n",
    "plt.errorbar(freqs, bs, 6*sigmas, fmt='.r', label='COBE FIRAS')#  plot 95% confidence error bars (+/- 3 sigma)\n",
    "plt.legend()\n",
    "\n",
    "plt.ylabel('B ($W m^{-2} Hz^{-1}str^{-1}$)')\n",
    "plt.xlabel('$f (Hz)$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "**Question 3**\n",
    "\n",
    "How well does this fit predict the irradiance distribution? How could you get a better estimate?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Next, we'll fit this using a method called non-linear least squares. You can read about it [here](https://en.wikipedia.org/wiki/Non-linear_least_squares).\n",
    "\n",
    "The specific implementation we'll use is from scipy (\"Scientific Python\"). This is an open source library that contains an incredibly powerful set of tools for scientific analysis and it is absolutely free! The scipy documentation for the `curve_fit` function we will use is [here](https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html)\n",
    "\n",
    "The way non-linear least squares (as well a most regression methods work) is that you provide a set of data points of the dependent variable (`ydata`) at a set of independent variables (`xdata`). You also specify the uncertainty in your measurements of the dependent variable `sigma`.\n",
    "\n",
    "You supply a function that takes the independent variable as it first variable and any additional parameters as additional arguments. This function is one that you expect to fit your data. In our case we will use the black body spectrum `bb_spectrum_freq`\n",
    "\n",
    "**The output of the fitting process is the parameters of this function that best fit your data, and some error bound on those parameters**\n",
    "\n",
    "\n",
    "We can also seed the fitting process with a good guess of our parameters. In our case we can use the temperature $T_{peak}$ that we computed above that we computed above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Run this cell\n",
    "\n",
    "parameters, covariance = curve_fit(\n",
    "    f=bb_spectrum_freq,  # this is a function that computes the theoretical irradiance at any frequency for a given temp\n",
    "    xdata=freqs,   # this is the list of frequencies at which irradiance was measures\n",
    "    ydata=bs,      # this is the measure irradiance\n",
    "    sigma=sigmas,  # this is the standard error of each data point\n",
    "    p0=[T_peak],   # this is the initial guess we give the fitting algorithm\n",
    ")\n",
    "T_fit = parameters[0]\n",
    "sigma_T_fit = np.sqrt(covariance[0][0])\n",
    "print('T_fit = {:0.4f}K\\nsigma_T = {:0.4}K'.format(T_fit, sigma_T_fit))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "**Question 4**\n",
    "\n",
    "What are the values of T and your uncertainty (sigma)? Does your initial estimate fall within this range?\n",
    "\n",
    "Now plot the theoretical curve at `T=T_fit` against the data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "outputs": [],
   "source": [
    "fig = plt.figure(figsize=(10, 10))\n",
    "ax1 = fig.add_subplot(111)\n",
    "theory_f = np.linspace(5e10, 8e11, 1000)\n",
    "ax1.plot(\n",
    "    theory_f,\n",
    "    bb_spectrum_freq(theory_f, T_fit),\n",
    "    '-k', label='theory, T={:0.3f}'.format(T_fit)\n",
    ")\n",
    "ax1.errorbar(freqs, bs, 2*3*sigmas, fmt='.r', label=\"COBE FIRAS\") # plot 95% confidence error bars (+/- 3 sigma)\n",
    "plt.legend()\n",
    "plt.ylabel('B ($W m^{-2} Hz^{-1}str^{-1}$)')\n",
    "plt.xlabel('$f (Hz)$')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "**Question 5**\n",
    "\n",
    "How does this theoretical curve compare to the the one you initially found? Zoom in to a few points - does the theoretical curve fall within the bounds of the error bars? This experiment was very precise, you may have to zoom in far!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Answer:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Question 6**\n",
    "\n",
    "Why do you think the temperature we determined with this data is a better match for the data? How could you quantify how closely a theoretical curve matches experimentally observed data?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "metadata": {
     "deletable": false,
     "editable": false,
     "trusted": true
    }
   },
   "source": [
    "Answer:"
   ]
  },
  {
   "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}