compare_websky_agora_sz.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OMP_NUM_THREADS set to 256\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "\n",
    "# Automatically set OMP_NUM_THREADS to all available cores\n",
    "n_cores = os.cpu_count()\n",
    "os.environ[\"OMP_NUM_THREADS\"] = str(n_cores)\n",
    "print(f\"OMP_NUM_THREADS set to {n_cores}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pysm3 import Sky\n",
    "from pysm3 import units as u\n",
    "\n",
    "import healpy as hp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from astropy.table import Table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "log = logging.getLogger(\"pysm3\")\n",
    "log.setLevel(logging.DEBUG)\n",
    "handler = logging.StreamHandler()\n",
    "log.addHandler(handler)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def joint_limits(freq_maps, percentile=99.5):\n",
    "    \"\"\"Return symmetric color limits for a dictionary of Quantity maps.\"\"\"\n",
    "    stacked = np.concatenate([freq_maps[model][0].value.ravel() for model in freq_maps])\n",
    "    vmax = np.percentile(np.abs(stacked), percentile)\n",
    "    if vmax <= 0:\n",
    "        vmax = np.max(np.abs(stacked))\n",
    "    if vmax == 0:\n",
    "        vmax = 1.0\n",
    "    vmin = -vmax\n",
    "    return vmin, vmax\n",
    "\n",
    "\n",
    "def remove_monopole(component):\n",
    "    \"\"\"Subtract the mean signal and return the centered Quantity plus the monopole.\"\"\"\n",
    "    data = component.value\n",
    "    mono = float(np.mean(data))\n",
    "    centered = (data - mono) * component.unit\n",
    "    return centered, mono"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Reference frequencies from the literature\n",
    "- Birkinshaw (1999, *Physics Reports* 310, 97) derives the thermal SZ (tSZ) spectrum, showing the intensity dip near 128 GHz, a null at ≃217 GHz, and the positive maximum around 370 GHz, while the kinematic SZ (kSZ) contribution remains frequency-independent in thermodynamic CMB units.\n",
    "- The review by Carlstrom, Holder & Reese (2002, *ARA&A* 40, 643) emphasises that multi-band SZ programs pair a low-frequency tSZ-dominated band (≲100 GHz) with bands on both sides of the 217 GHz null to isolate kSZ and foregrounds.\n",
    "- Operational surveys such as the SPT-SZ camera explicitly observe at 95 GHz, 150 GHz, and 220 GHz to separate tSZ and kSZ signals (South Pole Telescope instrument description, retrieved Jan 2026).\n",
    "\n",
    "Guided by those references, we probe 90, 150, 217, and 280 GHz below, bracketing the tSZ null and matching common survey bands."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Thermal SZ Comparisons\n",
    "\n",
    "Summary plots and statistics for the thermal component across WebSky and Agora realizations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Access data from /global/cfs/cdirs/cmb/www/pysm-data/websky/0.4/tsz_8192_hp.fits\n",
      "Access data from /global/cfs/cdirs/cmb/www/pysm-data/websky/0.4/tsz_8192_hp.fits\n",
      "Access data from /global/cfs/cdirs/cmb/www/pysm-data/agora/agora_utszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18.fits\n",
      "Access data from /global/cfs/cdirs/cmb/www/pysm-data/agora/agora_utszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18.fits\n",
      "Access data from /global/cfs/cdirs/cmb/www/pysm-data/agora/agora_ltszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18_lensed.fits\n",
      "Access data from /global/cfs/cdirs/cmb/www/pysm-data/agora/agora_ltszNG_bahamas80_bnd_unb_1.0e+12_1.0e+18_lensed.fits\n"
     ]
    }
   ],
   "source": [
    "models = [\"tsz1\", \"tsz2\", \"tsz3\"]\n",
    "model_titles = {\n",
    "    \"tsz1\": \"WebSky tSZ\",\n",
    "    \"tsz2\": \"Agora tSZ\",\n",
    "    \"tsz3\": \"Agora tSZ (lensed)\",\n",
    "}\n",
    "\n",
    "skys = {model: Sky(nside=512, preset_strings=[model]) for model in models}\n",
    "\n",
    "freqs_ghz = [90, 150, 217, 280]\n",
    "maps = {\n",
    "    freq: {model: skys[model].get_emission(freq * u.GHz) for model in models}\n",
    "    for freq in freqs_ghz\n",
    "}\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Monopoles removed (units in uK_RJ ):\n",
      "  90 GHz -> WebSky tSZ: -4.390e+00, Agora tSZ: -7.281e+00, Agora tSZ (lensed): -7.281e+00\n",
      "  150 GHz -> WebSky tSZ: -1.853e+00, Agora tSZ: -3.073e+00, Agora tSZ (lensed): -3.073e+00\n",
      "  217 GHz -> WebSky tSZ: -8.006e-03, Agora tSZ: -1.328e-02, Agora tSZ (lensed): -1.328e-02\n",
      "  280 GHz -> WebSky tSZ: 6.028e-01, Agora tSZ: 9.998e-01, Agora tSZ (lensed): 9.998e-01\n"
     ]
    },
    {
     "data": {