{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Bandpass integration\n",
    "\n",
    "The `get_emission` method implemented by any subclass of `Model`, including `Sky`, allows a `numpy` array of frequencies instead of a single frequency value to define a tophat bandpass and optionally a `weights` array of the same length to define a custom weight.\n",
    "PySM normalizes the `weights` array to unit integral then generates the emission at each of the specified frequencies, it multiplies it by the corresponding weight, and integrates them with the Trapezoidal rule.\n",
    "Instead of having generating all the maps as PySM 2, PySM 3 keeps in memory just 1 array for the output map and then creates the maps at a specific frequency one at a time and accomulates it properly into the output array."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "nbsphinx": "hidden"
   },
   "outputs": [],
   "source": [
    "import pysm\n",
    "import pysm.units as u\n",
    "import healpy as hp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sky = pysm.Sky(nside=128, preset_strings=[\"d1\", \"s1\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_100GHz_delta = sky.get_emission(100 * u.GHz)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bandpass_frequencies = np.linspace(95, 105, 11) * u.GHz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bandpass_frequencies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_100GHz_tophat = sky.get_emission(bandpass_frequencies)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(map_100GHz_delta[0] - map_100GHz_tophat[0], min=-1, max=1, title=\"I map difference\", unit=map_100GHz_tophat.unit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "bandpass_weights = np.array([2,3,5,9,11,11.5,11,9,5,3,2])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(bandpass_frequencies,bandpass_weights);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "map_100GHz_bandpass = sky.get_emission(bandpass_frequencies, bandpass_weights)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(map_100GHz_bandpass[0] - map_100GHz_tophat[0], min=-1, max=1, title=\"I map difference\", unit=map_100GHz_tophat.unit)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "PySM 3",
   "language": "python",
   "name": "pysm3"
  },
  "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.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
