planck_namaster_with_ff.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "4db40b20",
   "metadata": {},
   "source": [
    "# Prepare PySM 'f2' free-free model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c57500c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import healpy as hp\n",
    "import pysm3\n",
    "from pathlib import Path\n",
    "import os\n",
    "import pysm3.units as u\n",
    "import pymaster as nmt\n",
    "from scipy.optimize import curve_fit\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84dfbcec",
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "import os\n",
    "\n",
    "os.environ[\n",
    "    \"OMP_NUM_THREADS\"\n",
    "] = \"64\"  # for jupyter.nersc.gov otherwise the notebook only uses 2 cores\n",
    "\n",
    "output_files = []\n",
    "\n",
    "from datetime import date\n",
    "today = date.today()\n",
    "version = today.strftime(\"%Y.%m.%d\")\n",
    "version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2d921e57",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "ffbd7ea8",
   "metadata": {},
   "source": [
    "## Reading the HE free-free EM map, Processing"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f9bad52f",
   "metadata": {},
   "outputs": [],
   "source": [
    "imapfile = Path(\"HE_EM_mean_std.fits\")\n",
    "if not imapfile.exists():\n",
    "    os.system(f\"wget -O {imapfile} https://zenodo.org/records/10523170/files/EM_mean_std.fits?download=1\")\n",
    "else:\n",
    "    print(f\"Using existing file {imapfile}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8cb74ffd",
   "metadata": {},
   "outputs": [],
   "source": [
    "HE_EM=hp.read_map(imapfile) #EM is the first column, cm^-6 pc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c0c08398",
   "metadata": {},
   "outputs": [],
   "source": [
    "hp.mollview(HE_EM, title='HE_EM', unit=r'${\\rm cm}^{-6} {\\rm pc}$', cmap='turbo',norm='hist')\n",
    "hp.graticule(dpar=30)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8f541f5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "nside=hp.get_nside(HE_EM)\n",
    "print(nside)\n",
    "lmax=3*nside-1\n",
    "print(lmax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f07e280",
   "metadata": {},
   "source": [
    "#### Converting the Emission Measure map to Intensity in brightness temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bf211327",
   "metadata": {},
   "outputs": [],
   "source": [
    "#this function converts the EM map to a T_uRJ map\n",
    "\n",
    "def convert_EM_to_TRJ(EMmap,freq,TeMap=7000):   #we use 7,000 K as the average value of the Electron temperature in the ISM (Dickinson et al 2003, Planck 2015 X)\n",
    "\n",
    "    gff=np.log((np.exp(5.960-((np.sqrt(3)/np.pi)*np.log(freq*((TeMap*1e-4)**(-3/2))))))+np.exp(1))\n",
    "\n",
    "    tau=0.05468*(TeMap**(-3/2))*EMmap*(freq**(-2))*gff\n",
    "\n",
    "    TRJmap=1e6*TeMap*(1-np.exp(-tau))\n",
    "\n",
    "    return TRJmap\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ef30cfbd",
   "metadata": {},
   "outputs": [],
   "source": [
    "HE_ff_T=convert_EM_to_TRJ(HE_EM,freq=30.0)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b0cad0f5",
   "metadata": {},
   "source": [
    "# HE_ff_T Analysis with Planck Mask and NaMaster\n",
    "This notebook loads a HEALPix temperature-only map (`HE_ff_T`), applies a Planck galactic mask, computes the angular power spectrum using NaMaster, and compares the original and filtered maps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec412789",
   "metadata": {},
   "outputs": [],
   "source": [
    "import healpy as hp\n",
    "import numpy as np\n",
    "import pymaster as nmt\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Load HE_ff_T temperature-only HEALPix map\n",
    "# Replace 'path_to_HE_ff_T.fits' with your file path\n",
    "HE_ff_T = hp.read_map('path_to_HE_ff_T.fits')\n",
    "nside = hp.get_nside(HE_ff_T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ada45dd7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot original HE_ff_T map\n",
    "hp.mollview(HE_ff_T, title='Original HE_ff_T map', unit='K')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a053a003",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Load Planck galactic mask (NSIDE=2048)\n",
    "mask_2048 = hp.read_map('COM_Mask_GalPlane-apo2_2048_R2.00.fits')\n",
    "\n",
    "# Downgrade to match HE_ff_T resolution\n",
    "mask = hp.ud_grade(mask_2048, nside_out=nside)\n",
    "mask = np.where(mask > 0.9, 1.0, 0.0)  # Binary mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e08b1241",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apodize mask (cosine apodization of 2 degrees)\n",
    "mask_apo = nmt.mask_apodization(mask, aposize_deg=2.0, apotype='C1')\n",
    "\n",
    "# Plot apodized mask\n",
    "hp.mollview(mask_apo, title='Apodized Mask')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7ed7bc58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define NaMaster field\n",
    "field = nmt.NmtField(mask_apo, [HE_ff_T])\n",
    "\n",
    "# Define binning scheme (e.g., 0-128 in steps of 10)\n",
    "binner = nmt.NmtBin(nside, bpws=np.arange(0, 130, 10))\n",
    "\n",
    "# Compute coupled Cl\n",
    "cl_coupled = nmt.compute_coupled_cell(field, field)\n",
    "\n",
    "# Compute and apply coupling matrix to decouple\n",
    "workspace = nmt.NmtWorkspace()\n",
    "workspace.compute_coupling_matrix(field, field, binner)\n",
    "cl_decoupled = workspace.decouple_cell(cl_coupled)\n",
    "\n",
    "# Plot power spectrum\n",
    "ells = binner.get_effective_ells()\n",
    "plt.figure()\n",
    "plt.plot(ells, cl_decoupled[0], label='Temperature Cl')\n",
    "plt.xlabel(r'$\\ell$')\n",
    "plt.ylabel(r'$C_\\ell$')\n",
    "plt.title('Angular Power Spectrum (Decoupled)')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7b84f58a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Reconstruct smoothed map (low-pass filter example)\n",
    "# Apply mask before transform\n",
    "alm = hp.map2alm(HE_ff_T * mask_apo)\n",
    "# Zero out modes above ell=50\n",
    "lmax_alm = hp.Alm.getlmax(len(alm))\n",
    "filter_array = np.array([1 if ell < 50 else 0 for ell in range(lmax_alm + 1)])\n",
    "hp.almxfl(alm, filter_array)\n",
    "filtered_map = hp.alm2map(alm, nside)\n",
    "\n",
    "# Plot original vs filtered maps side by side\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "plt.subplot(121)\n",
    "hp.mollview(HE_ff_T, title='Original Map', unit='K', sub=(1,2,1))\n",
    "plt.subplot(122)\n",
    "hp.mollview(filtered_map, title='Filtered Map (low-l)', unit='K', sub=(1,2,2))\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}