diff --git a/docs/preprocess-templates/freefree/small_scale_freefree_f2_generate_template.ipynb b/docs/preprocess-templates/freefree/small_scale_freefree_f2_generate_template.ipynb new file mode 100644 index 00000000..4beb3caa --- /dev/null +++ b/docs/preprocess-templates/freefree/small_scale_freefree_f2_generate_template.ipynb @@ -0,0 +1,1401 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "c25cf661", + "metadata": {}, + "source": [ + "# Prepare PySM 'f2' free-free model template" + ] + }, + { + "cell_type": "markdown", + "id": "b8cf6188", + "metadata": {}, + "source": "This notebook implements the small-scale injection workflow for free-free intensity maps (intensity only, no polarization components).\n\n**Goal:** Generate a high-resolution free-free template by injecting statistically consistent small-scale structure while preserving the large-scale morphology.\n\n**Input data:** The emission-measure (EM) mean/variance map at nside=256 from the all-sky Galactic Faraday/electron-density reconstruction of [Hutschenreuter et al. (2024)](https://doi.org/10.1051/0004-6361/202346740). Throughout this notebook, we refer to this as the HE map, following Hutschenreuter & En\u00dflin (2020), which presented the original reconstruction framework. The FITS product `EM_mean_std.fits` is available on [Zenodo](https://zenodo.org/records/10523170).\n\n**Method:** The small-scale injection follows the logpoltens formalism used in modern PySM3 foreground templates, adapted here to intensity-only free-free: [PanEx Galactic Science Group et al. (2025)](https://arxiv.org/abs/2502.20452).\n\n**Workflow summary:**\n1. Download and convert Emission Measure map to brightness temperature at 30 GHz\n2. Transform to log-space: $i = \\ln(I + \\epsilon)$\n3. Fit power-law to angular power spectrum $\\mathcal{D}_\\ell$ in $\\ell = 30\\text{--}150$\n4. Construct extrapolated small-scale power spectrum with sigmoid high-pass filter\n5. Build spatially-varying modulation map from patch-based spectral analysis\n6. Filter large-scale map with $\\sqrt{1-w}$ low-pass; generate small-scale realization filtered with $\\sqrt{w}$ high-pass\n7. Combine: output = large-scale + modulated small-scale\n8. Transform back to linear space: $I = e^{i} - \\epsilon$" + }, + { + "cell_type": "markdown", + "id": "f98c9439", + "metadata": {}, + "source": [ + "## 1. Setup" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1e9083ac", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:41:56.260867Z", + "iopub.status.busy": "2026-03-24T23:41:56.260436Z", + "iopub.status.idle": "2026-03-24T23:41:58.842311Z", + "shell.execute_reply": "2026-03-24T23:41:58.841020Z" + } + }, + "outputs": [], + "source": [ + "import os\n", + "\n", + "os.environ.setdefault(\n", + " \"OMP_NUM_THREADS\",\n", + " os.environ.get(\"SLURM_CPUS_PER_TASK\") or str(os.cpu_count() or 1),\n", + ")\n", + "print(\"OMP_NUM_THREADS =\", os.environ[\"OMP_NUM_THREADS\"])\n", + "\n", + "from pathlib import Path\n", + "from datetime import date\n", + "import urllib.request\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import healpy as hp\n", + "from scipy.optimize import curve_fit\n", + "import pymaster as nmt\n", + "\n", + "import pysm3\n", + "import pysm3.units as u\n", + "from pysm3.utils import sigmoid, add_metadata\n", + "\n", + "%matplotlib inline\n", + "hp.disable_warnings()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ebf4a08", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:41:58.846494Z", + "iopub.status.busy": "2026-03-24T23:41:58.846112Z", + "iopub.status.idle": "2026-03-24T23:41:58.854650Z", + "shell.execute_reply": "2026-03-24T23:41:58.853426Z" + } + }, + "outputs": [], + "source": [ + "# ---- Configuration parameters ----\n", + "# All key parameters collected here for reproducibility.\n", + "\n", + "# Output resolution\n", + "nside_final = 2048\n", + "lmax_final = 3 * nside_final - 1\n", + "\n", + "# Power-law fit range (in multipole ell)\n", + "ell_fit_low = 30 # Lower bound for D_ell fit\n", + "ell_fit_high = 150 # Upper bound; also sets sigmoid transition center\n", + "\n", + "# Target small-scale spectral index (D_ell power law)\n", + "gamma_target = -1.5 # Steeper than fitted slope; controls small-scale falloff\n", + "\n", + "# Modulation map parameters\n", + "nsidepatches = 8 # Patch grid resolution (768 patches)\n", + "patch_radius_deg = 7.6 # Disc radius for each patch [deg]\n", + "smoothing_fwhm_deg = 5.5 # FWHM for smoothing modulation map [deg]\n", + "ell_pivot = 80 # Pivot multipole for modulation normalization\n", + "\n", + "# Random seed for reproducibility\n", + "RANDOM_SEED = 777\n", + "\n", + "# Version tag\n", + "version = date.today().strftime(\"%Y.%m.%d\")\n", + "\n", + "# Track output files\n", + "output_files = []\n", + "\n", + "# Paths\n", + "DATA_DIR = Path(\"ff_pysm_data\")\n", + "DATA_DIR.mkdir(parents=True, exist_ok=True)\n", + "\n", + "print(f\"Output: nside={nside_final}, lmax={lmax_final}\")\n", + "print(f\"Version: {version}\")" + ] + }, + { + "cell_type": "markdown", + "id": "33fd49ea", + "metadata": {}, + "source": "## 2. Download and read the input Emission Measure map\n\nThe input free-free tracer is the Emission Measure (EM) map at nside=256 from the all-sky Galactic Faraday/electron-density reconstruction of [Hutschenreuter et al. (2024)](https://doi.org/10.1051/0004-6361/202346740). The FITS file `EM_mean_std.fits` contains two columns: the posterior mean and standard deviation of EM in units of cm$^{-6}$ pc. We use the mean (first column)." + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c495589b", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:41:58.858148Z", + "iopub.status.busy": "2026-03-24T23:41:58.857830Z", + "iopub.status.idle": "2026-03-24T23:42:01.730888Z", + "shell.execute_reply": "2026-03-24T23:42:01.729667Z" + } + }, + "outputs": [], + "source": [ + "url = \"https://zenodo.org/records/10523170/files/EM_mean_std.fits?download=1\"\n", + "imapfile = DATA_DIR / \"HE_EM_mean_std.fits\"\n", + "\n", + "if not imapfile.exists():\n", + " print(f\"Downloading input map to {imapfile} ...\")\n", + " urllib.request.urlretrieve(url, imapfile)\n", + "else:\n", + " print(f\"Using existing file {imapfile}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f9e562ae", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:01.735086Z", + "iopub.status.busy": "2026-03-24T23:42:01.734568Z", + "iopub.status.idle": "2026-03-24T23:42:03.030277Z", + "shell.execute_reply": "2026-03-24T23:42:03.028816Z" + } + }, + "outputs": [], + "source": [ + "HE_EM = hp.read_map(imapfile) # EM mean, units: cm^-6 pc\n", + "\n", + "nside = hp.get_nside(HE_EM)\n", + "lmax = 3 * nside - 1\n", + "print(f\"Input map: nside={nside}, lmax={lmax}\")\n", + "\n", + "hp.mollview(HE_EM, title=\"HE EM Map\",\n", + " unit=r\"${\\rm cm}^{-6}\\,{\\rm pc}$\", cmap=\"turbo\", norm=\"hist\")\n", + "hp.graticule(dpar=30)" + ] + }, + { + "cell_type": "markdown", + "id": "764560c4", + "metadata": {}, + "source": [ + "The sharp horizontal boundary visible in the northern high-latitude region is inherited from the original H$\\alpha$ map used to construct this EM map, and results from the stitching of two independent H$\\alpha$ surveys." + ] + }, + { + "cell_type": "markdown", + "id": "b596eae8", + "metadata": {}, + "source": [ + "## 3. Convert Emission Measure to brightness temperature\n", + "\n", + "We convert the EM map to free-free brightness temperature in $\\mu\\mathrm{K_{RJ}}$ at a reference frequency of 30 GHz using the standard free-free emission formula.\n", + "\n", + "The optical depth is:\n", + "$$\n", + "\\tau_{\\rm ff} = 0.05468 \\; T_e^{-3/2} \\; \\mathrm{EM} \\; \\nu^{-2} \\; g_{\\rm ff}\n", + "$$\n", + "\n", + "where $g_{\\rm ff}$ is the Gaunt factor approximation from [Draine (2011)](https://ui.adsabs.harvard.edu/abs/2011piim.book.....D):\n", + "$$\n", + "g_{\\rm ff} = \\ln\\!\\left[\\exp\\!\\left(5.960 - \\frac{\\sqrt{3}}{\\pi}\\ln\\!\\left(\\nu_{\\rm GHz}\\,(T_e / 10^4\\,\\mathrm{K})^{-3/2}\\right)\\right) + e\\right]\n", + "$$\n", + "\n", + "The brightness temperature is:\n", + "$$\n", + "T_{\\rm RJ} = T_e \\,(1 - e^{-\\tau_{\\rm ff}})\n", + "$$\n", + "\n", + "We adopt $T_e = 7000\\,\\mathrm{K}$ as a representative mean electron temperature for the diffuse ISM, following [Dickinson et al. (2003)](https://doi.org/10.1046/j.1365-8711.2003.06605.x) and [Planck 2015 X](https://doi.org/10.1051/0004-6361/201525967)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d74ac2bc", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:03.042117Z", + "iopub.status.busy": "2026-03-24T23:42:03.041669Z", + "iopub.status.idle": "2026-03-24T23:42:03.048952Z", + "shell.execute_reply": "2026-03-24T23:42:03.047646Z" + } + }, + "outputs": [], + "source": [ + "def convert_EM_to_TRJ(EM_map, freq_ghz, Te=7000):\n", + " \"\"\"Convert Emission Measure to Rayleigh-Jeans brightness temperature.\n", + "\n", + " Parameters\n", + " ----------\n", + " EM_map : array\n", + " Emission measure in cm^-6 pc.\n", + " freq_ghz : float\n", + " Frequency in GHz.\n", + " Te : float\n", + " Electron temperature in K (default 7000 K).\n", + "\n", + " Returns\n", + " -------\n", + " T_RJ : array\n", + " Brightness temperature in uK_RJ.\n", + " \"\"\"\n", + " gff = np.log(\n", + " np.exp(5.960 - (np.sqrt(3) / np.pi) * np.log(freq_ghz * (Te * 1e-4) ** (-1.5)))\n", + " + np.exp(1)\n", + " )\n", + " tau = 0.05468 * Te ** (-1.5) * EM_map * freq_ghz ** (-2) * gff\n", + " return 1e6 * Te * (1 - np.exp(-tau)) # uK_RJ" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3004ad66", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:03.052866Z", + "iopub.status.busy": "2026-03-24T23:42:03.052437Z", + "iopub.status.idle": "2026-03-24T23:42:03.810144Z", + "shell.execute_reply": "2026-03-24T23:42:03.808777Z" + } + }, + "outputs": [], + "source": [ + "HE_ff_T = convert_EM_to_TRJ(HE_EM, freq_ghz=30.0)\n", + "print(f\"T_RJ range: [{np.min(HE_ff_T):.4f}, {np.max(HE_ff_T):.4f}] uK_RJ\")\n", + "\n", + "hp.mollview(HE_ff_T, cmap=\"turbo\",\n", + " title=\"HE Free-free Intensity at 30 GHz\",\n", + " unit=r\"$T\\;[\\mu\\mathrm{K_{RJ}}]$\", norm=\"hist\")\n", + "hp.graticule(dpar=30)" + ] + }, + { + "cell_type": "markdown", + "id": "1dba7d6a", + "metadata": {}, + "source": [ + "## 4. Log-space transformation\n", + "\n", + "We work in log-space to make the field more Gaussian and enable additive small-scale injection. The forward transform is:\n", + "\n", + "$$\n", + "i = \\ln(I + \\epsilon)\n", + "$$\n", + "\n", + "where $\\epsilon$ is a small positive floor chosen from the data to handle zeros without distorting positive values. The inverse is:\n", + "\n", + "$$\n", + "I = e^{i} - \\epsilon\n", + "$$\n", + "\n", + "**Note:** This differs from the full log-polarization-tensor formalism used for polarized components (synchrotron/dust IQU), which involves $\\ln(\\det T)$ and hyperbolic functions. Here we only have intensity, so a simple $\\ln(I + \\epsilon)$ suffices." + ] + }, + { + "cell_type": "markdown", + "id": "63c1cbd9", + "metadata": {}, + "source": [ + "### Choosing epsilon\n", + "\n", + "The `choose_eps` function selects a log-floor epsilon from the data itself as the maximum of three candidates:\n", + "- `abs_floor` (default $10^{-30}$): absolute minimum to avoid numerical issues\n", + "- `frac_min` $\\times$ min(positive values) (default 0.1): fraction of the smallest positive pixel\n", + "- `frac_median` $\\times$ median(positive values) (default $10^{-12}$): fraction of the median\n", + "\n", + "This ensures $\\epsilon$ is large enough to regularize zeros but small enough not to distort the dynamic range." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3bca7758", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:03.821024Z", + "iopub.status.busy": "2026-03-24T23:42:03.820609Z", + "iopub.status.idle": "2026-03-24T23:42:03.828702Z", + "shell.execute_reply": "2026-03-24T23:42:03.827575Z" + } + }, + "outputs": [], + "source": [ + "def choose_eps(m, abs_floor=1e-30, frac_min=0.1, frac_median=1e-12):\n", + " \"\"\"Choose a log-floor eps from the data.\n", + "\n", + " Returns max(abs_floor, frac_min * min(positive), frac_median * median(positive)).\n", + " \"\"\"\n", + " m = np.asarray(m, dtype=float)\n", + " pos = m[np.isfinite(m) & (m > 0)]\n", + " if pos.size == 0:\n", + " return abs_floor\n", + " return max(abs_floor, frac_min * float(pos.min()), frac_median * float(np.median(pos)))\n", + "\n", + "\n", + "def map_to_log(m, eps):\n", + " \"\"\"Forward log transform: i = ln(I + eps).\"\"\"\n", + " return np.log(np.asarray(m, dtype=float) + eps)\n", + "\n", + "\n", + "def log_to_map(log_m, eps):\n", + " \"\"\"Inverse log transform: I = exp(i) - eps.\"\"\"\n", + " return np.exp(np.asarray(log_m, dtype=float)) - eps" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6b12862", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:03.832294Z", + "iopub.status.busy": "2026-03-24T23:42:03.831996Z", + "iopub.status.idle": "2026-03-24T23:42:03.867367Z", + "shell.execute_reply": "2026-03-24T23:42:03.866076Z" + } + }, + "outputs": [], + "source": [ + "eps = choose_eps(HE_ff_T)\n", + "\n", + "n_pix = HE_ff_T.size\n", + "n_zero = int(np.sum(HE_ff_T == 0))\n", + "n_neg = int(np.sum(HE_ff_T < 0))\n", + "n_nonfinite = int(np.sum(~np.isfinite(HE_ff_T)))\n", + "print(f\"eps = {eps:.4e}\")\n", + "print(f\"Total pixels: {n_pix}, zeros: {n_zero}, negative: {n_neg}, non-finite: {n_nonfinite}\")\n", + "\n", + "log_he_ff = map_to_log(HE_ff_T, eps)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0fd2ddcd", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:03.870298Z", + "iopub.status.busy": "2026-03-24T23:42:03.870028Z", + "iopub.status.idle": "2026-03-24T23:42:04.882507Z", + "shell.execute_reply": "2026-03-24T23:42:04.881220Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 5))\n", + "hp.mollview(HE_ff_T, cmap=\"turbo\", norm=\"hist\",\n", + " title=\"Original Free-Free Map ($I$)\", sub=(1, 2, 1))\n", + "hp.mollview(log_he_ff, cmap=\"turbo\", norm=\"hist\",\n", + " title=\"Log-space Free-Free Map ($i$)\", sub=(1, 2, 2))" + ] + }, + { + "cell_type": "markdown", + "id": "329731fc", + "metadata": {}, + "source": [ + "## 5. Power spectrum and power-law fit\n", + "\n", + "We compute the angular power spectrum $C_\\ell$ of the log-space map and fit a power-law model\n", + "\n", + "$$\n", + "\\mathcal{D}_\\ell \\equiv \\frac{\\ell(\\ell+1)}{2\\pi} C_\\ell = A\\,\\ell^{\\gamma}\n", + "$$\n", + "\n", + "in the multipole range $\\ell \\in [30, 150]$. This range avoids:\n", + "- The lowest multipoles ($\\ell < 30$) dominated by large-scale Galactic plane morphology and sampling variance\n", + "- Multipoles above $\\ell = 150$ where the input map (nside=256) is affected by pixel smoothing\n", + "\n", + "The fitted amplitude and slope are then used to construct an extrapolated power law for small scales. We adopt a **target spectral index** $\\gamma_{\\rm target} = -1.5$, steeper than the fitted slope. The extrapolated amplitude is anchored at $\\ell_{\\rm fit,high}$ to ensure continuity:\n", + "\n", + "$$\n", + "A_{\\rm extrap} = |A_{\\rm fit}| \\cdot \\ell_{\\rm fit,high}^{\\;\\gamma_{\\rm fit} - \\gamma_{\\rm target}}\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ce5e0d3", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:04.893282Z", + "iopub.status.busy": "2026-03-24T23:42:04.892886Z", + "iopub.status.idle": "2026-03-24T23:42:04.900126Z", + "shell.execute_reply": "2026-03-24T23:42:04.898872Z" + } + }, + "outputs": [], + "source": [ + "def run_anafast(m, lmax):\n", + " \"\"\"Compute power spectrum and return ell, cl_norm, cl, dl.\"\"\"\n", + " cl = hp.anafast(m, lmax=lmax, use_pixel_weights=True)\n", + " ell = np.arange(lmax + 1, dtype=float)\n", + " cl_norm = ell * (ell + 1) / (2 * np.pi)\n", + " cl_norm[0] = 1\n", + " dl = cl_norm * cl\n", + " dl[0] = 1\n", + " return ell, cl_norm, cl, dl\n", + "\n", + "\n", + "def powerlaw_model(ell, A, gamma):\n", + " \"\"\"Power-law model: D_ell = A * ell^gamma.\"\"\"\n", + " return A * ell ** gamma" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c20ba38d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:04.903185Z", + "iopub.status.busy": "2026-03-24T23:42:04.902888Z", + "iopub.status.idle": "2026-03-24T23:42:05.369185Z", + "shell.execute_reply": "2026-03-24T23:42:05.368076Z" + } + }, + "outputs": [], + "source": [ + "ell_HE, cl_norm_HE, cl_HE, dl_HE = run_anafast(log_he_ff, lmax)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5e1555e9", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:05.372306Z", + "iopub.status.busy": "2026-03-24T23:42:05.372042Z", + "iopub.status.idle": "2026-03-24T23:42:07.150804Z", + "shell.execute_reply": "2026-03-24T23:42:07.149686Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 5))\n", + "\n", + "plt.subplot(1, 2, 1)\n", + "plt.loglog(ell_HE, cl_HE, color=\"darkorange\")\n", + "plt.grid()\n", + "plt.title(r\"$C_\\ell$ for log-space free-free\")\n", + "plt.xlabel(r\"$\\ell$\")\n", + "plt.ylabel(r\"$C_\\ell$\")\n", + "\n", + "plt.subplot(1, 2, 2)\n", + "plt.loglog(ell_HE, dl_HE, color=\"darkblue\")\n", + "plt.xlim([2, lmax])\n", + "plt.grid()\n", + "plt.title(r\"$\\mathcal{D}_\\ell$ for log-space free-free\")\n", + "plt.xlabel(r\"$\\ell$\")\n", + "plt.ylabel(r\"$\\ell(\\ell+1)C_\\ell/2\\pi$\")\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a096a325", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:07.154839Z", + "iopub.status.busy": "2026-03-24T23:42:07.154502Z", + "iopub.status.idle": "2026-03-24T23:42:07.165669Z", + "shell.execute_reply": "2026-03-24T23:42:07.164632Z" + } + }, + "outputs": [], + "source": [ + "# Fit D_ell power law in [ell_fit_low, ell_fit_high]\n", + "xdata = np.arange(ell_fit_low, ell_fit_high)\n", + "ydata = xdata * (xdata + 1) / (2 * np.pi) * cl_HE[xdata]\n", + "(A_fit, gamma_fit), cov = curve_fit(powerlaw_model, xdata, ydata)\n", + "print(f\"Fitted: A={A_fit:.4e}, gamma={gamma_fit:.4f}\")\n", + "\n", + "# Extrapolated amplitude anchored at ell_fit_high with target slope\n", + "A_extrap = np.abs(A_fit) * ell_fit_high ** (gamma_fit - gamma_target)\n", + "print(f\"Extrapolated: A_extrap={A_extrap:.4e}, gamma_target={gamma_target}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "aa4c9667", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:07.169127Z", + "iopub.status.busy": "2026-03-24T23:42:07.168719Z", + "iopub.status.idle": "2026-03-24T23:42:08.249644Z", + "shell.execute_reply": "2026-03-24T23:42:08.248543Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(10, 5))\n", + "plt.loglog(ell_HE, dl_HE, color=\"darkblue\", label=\"Data\")\n", + "\n", + "scaling = np.zeros_like(ell_HE)\n", + "scaling[2:] = powerlaw_model(ell_HE[2:], A_extrap, gamma_target)\n", + "plt.plot(ell_HE[2:], scaling[2:], color=\"tomato\", ls=\"--\",\n", + " label=rf\"Extrapolated $\\gamma={gamma_target}$\")\n", + "\n", + "plt.axvline(ell_fit_low, ls=\":\", color=\"gray\", label=rf\"$\\ell={ell_fit_low}$\")\n", + "plt.axvline(ell_fit_high, ls=\":\", color=\"gray\", label=rf\"$\\ell={ell_fit_high}$\")\n", + "plt.xlim([2, 3 * lmax])\n", + "plt.legend()\n", + "plt.title(\"Power-law fit to log-space free-free TT spectrum\")\n", + "plt.xlabel(r\"$\\ell$\")\n", + "plt.ylabel(r\"$\\ell(\\ell+1)C_\\ell/2\\pi$\")" + ] + }, + { + "cell_type": "markdown", + "id": "b1360279", + "metadata": {}, + "source": [ + "## 6. Sigmoid filter and small-scale power spectrum\n", + "\n", + "We use a sigmoid function (imported from `pysm3.utils`) to smoothly separate large and small scales in harmonic space. The sigmoid transitions from 0 to 1 around $\\ell_0 = \\ell_{\\rm fit,high} = 150$ with width $\\Delta = \\ell_0/10 = 15$.\n", + "\n", + "The harmonic-space filters are designed to be **power-conserving**:\n", + "- Large-scale filter: $f_{\\rm LS}(\\ell) = \\sqrt{1 - w(\\ell)}$\n", + "- Small-scale filter: $f_{\\rm SS}(\\ell) = \\sqrt{w(\\ell)}$\n", + "\n", + "Since $f_{\\rm LS}^2 + f_{\\rm SS}^2 = 1$, the total power is preserved in the cross-fade between the original large-scale map and the synthetic small-scale realization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdf642f7", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:08.253035Z", + "iopub.status.busy": "2026-03-24T23:42:08.252741Z", + "iopub.status.idle": "2026-03-24T23:42:08.555632Z", + "shell.execute_reply": "2026-03-24T23:42:08.554443Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(8, 5))\n", + "plt.plot(ell_HE, sigmoid(ell_HE, ell_fit_high, ell_fit_high / 10),\n", + " color=\"darkviolet\", linewidth=2)\n", + "plt.xlabel(r\"Multipole $\\ell$\", fontsize=14)\n", + "plt.ylabel(r\"Sigmoid weight $w(\\ell)$\", fontsize=14)\n", + "plt.title(f\"Sigmoid transition at $\\\\ell_0={ell_fit_high}$, width={ell_fit_high // 10}\")\n", + "plt.grid(True, linestyle=\"--\", alpha=0.6)\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2480b2a0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:08.558779Z", + "iopub.status.busy": "2026-03-24T23:42:08.558451Z", + "iopub.status.idle": "2026-03-24T23:42:08.566299Z", + "shell.execute_reply": "2026-03-24T23:42:08.565039Z" + } + }, + "outputs": [], + "source": [ + "ell_final = np.arange(lmax_final + 1)\n", + "cl_norm_final = ell_final * (ell_final + 1) / (2 * np.pi)\n", + "cl_norm_final[0] = 1\n", + "\n", + "# Extrapolated power-law D_ell, converted to C_ell\n", + "powerlawscaling = np.zeros(lmax_final + 1)\n", + "powerlawscaling[2:] = powerlaw_model(ell_final[2:], A_extrap, gamma_target)\n", + "cl_pl = powerlawscaling / cl_norm_final\n", + "\n", + "# Sigmoid filter at output resolution\n", + "w_final = sigmoid(ell_final, ell_fit_high, ell_fit_high / 10)\n", + "f_ss_final = np.sqrt(w_final)\n", + "\n", + "# Filtered small-scale C_ell (filter^2 = w_final)\n", + "cl_ss = cl_pl * w_final" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dfaddca7", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:08.569895Z", + "iopub.status.busy": "2026-03-24T23:42:08.569578Z", + "iopub.status.idle": "2026-03-24T23:42:09.202664Z", + "shell.execute_reply": "2026-03-24T23:42:09.201356Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(8, 5))\n", + "plt.loglog(ell_final, cl_ss * cl_norm_final, color=\"orangered\", linewidth=2,\n", + " label=\"Sigmoid-filtered (small scales only)\")\n", + "plt.loglog(ell_final, powerlawscaling, color=\"darkgreen\", ls=\"--\",\n", + " label=\"Full extrapolated power law\")\n", + "plt.xlabel(r\"Multipole $\\ell$\", fontsize=14)\n", + "plt.ylabel(r\"$\\mathcal{D}_\\ell$\", fontsize=14)\n", + "plt.title(\"Small-scale power spectrum\")\n", + "plt.legend()\n", + "plt.grid(True, linestyle=\"--\", alpha=0.6)\n", + "plt.xscale(\"log\")\n", + "plt.yscale(\"log\")\n", + "plt.xlim([30, 3 * lmax_final])\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "071e8f81", + "metadata": {}, + "source": [ + "## 7. Spatially-varying modulation map\n", + "\n", + "The small-scale realization from `synfast` is statistically isotropic. To capture the spatial variation of small-scale power across the sky (stronger near the Galactic plane, weaker at high latitudes), we construct a modulation map.\n", + "\n", + "**Procedure:**\n", + "1. Downgrade the log-space map to nside=128 for faster spectral estimation.\n", + "2. Define 768 sky patches on an nside=8 grid, each a disc of radius 7.6\u00b0.\n", + "3. For each patch, compute the masked power spectrum using NaMaster C2 apodization.\n", + "4. Fit a power law $A_d\\,(\\ell/100)^\\alpha$ to the binned spectrum in $\\ell \\in [50, 100]$.\n", + "5. Evaluate the fitted amplitude at a pivot multipole $\\ell_{\\rm pivot}=80$ and normalize by the full-sky spectrum at that multipole.\n", + "6. Accumulate $\\sqrt{C_\\ell^{\\rm patch}(\\ell_{\\rm pivot}) / C_\\ell^{\\rm full}(\\ell_{\\rm pivot})}$ weighted by the apodized patch mask.\n", + "7. Divide by effective coverage count and smooth with FWHM=5.5\u00b0.\n", + "\n", + "The patch radius of 7.6\u00b0 is large enough to estimate the spectrum at $\\ell \\sim 50\\text{--}100$ while being small enough to capture spatial variation. The smoothing FWHM of 5.5\u00b0 suppresses patch-to-patch noise." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5731e561", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:09.206515Z", + "iopub.status.busy": "2026-03-24T23:42:09.206248Z", + "iopub.status.idle": "2026-03-24T23:42:09.390125Z", + "shell.execute_reply": "2026-03-24T23:42:09.389150Z" + } + }, + "outputs": [], + "source": [ + "# Mask Galactic plane regions (log(T) > 4)\n", + "maskmod = np.zeros_like(log_he_ff)\n", + "maskmod[log_he_ff > 4] = 1.0\n", + "\n", + "# Downgrade map and mask to nside=128\n", + "log_he_ff2 = hp.ud_grade(log_he_ff, nside_out=128)\n", + "maskmod = hp.ud_grade(maskmod, nside_out=hp.get_nside(log_he_ff2))\n", + "\n", + "lmaxmod = 100\n", + "apo_mask = nmt.mask_apodization(maskmod, 5, apotype=\"C2\")\n", + "ell0, norm, cl0, dl0 = run_anafast(log_he_ff2 * maskmod, lmax=lmaxmod)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "20515882", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:09.393170Z", + "iopub.status.busy": "2026-03-24T23:42:09.392883Z", + "iopub.status.idle": "2026-03-24T23:42:09.398205Z", + "shell.execute_reply": "2026-03-24T23:42:09.397203Z" + } + }, + "outputs": [], + "source": [ + "centers = np.vstack(\n", + " hp.pix2vec(ipix=np.arange(hp.nside2npix(nsidepatches)), nside=nsidepatches)\n", + ").T\n", + "print(f\"Number of patches: {len(centers)}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "259d5d41", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:09.401038Z", + "iopub.status.busy": "2026-03-24T23:42:09.400760Z", + "iopub.status.idle": "2026-03-24T23:42:27.015176Z", + "shell.execute_reply": "2026-03-24T23:42:27.014184Z" + } + }, + "outputs": [], + "source": [ + "fit_model = lambda x, Ad, alpha: Ad * (x / 100) ** alpha\n", + "\n", + "\n", + "def bin_cell(cell, dig):\n", + " \"\"\"Bin power spectrum values by digitized ell bins.\"\"\"\n", + " cb, errb = [], []\n", + " for i in np.unique(dig):\n", + " msk = dig == i\n", + " cb.append(cell[msk].mean())\n", + " errb.append(cell[msk].std())\n", + " return np.array(cb), np.array(errb)\n", + "\n", + "\n", + "def bin_ell(ells, dig):\n", + " \"\"\"Compute mean ell and half-width for each bin.\"\"\"\n", + " lb, dl = [], []\n", + " for i in np.unique(dig):\n", + " msk = dig == i\n", + " lb.append(ells[msk].mean())\n", + " dl.append((ells[msk].max() - ells[msk].min()) / 2)\n", + " return np.array(lb), np.array(dl)\n", + "\n", + "\n", + "tmod = np.zeros_like(log_he_ff2)\n", + "n_eff = np.zeros_like(log_he_ff2)\n", + "\n", + "for ipix, c in enumerate(centers):\n", + " patch = np.zeros_like(log_he_ff2)\n", + " maskpixs = hp.query_disc(\n", + " nside=hp.get_nside(log_he_ff2), vec=c,\n", + " radius=np.radians(patch_radius_deg),\n", + " )\n", + " patch[maskpixs] = 1\n", + " apo_patch = nmt.mask_apodization(patch, 5, apotype=\"C2\")\n", + " fsky = apo_patch.sum() / apo_patch.size\n", + "\n", + " ellp, norm, clp, dlp = run_anafast(log_he_ff2 * apo_patch, lmax=lmaxmod)\n", + " digi = np.digitize(ellp, np.linspace(0, lmaxmod, 10))\n", + " dtt, errtt = bin_cell(clp, digi) / fsky\n", + " lb, delta_l = bin_ell(ellp, digi)\n", + "\n", + " lmaskt = np.logical_and(lb < 100, lb > 50)\n", + "\n", + " param_tt, _ = curve_fit(\n", + " fit_model, xdata=lb[lmaskt], ydata=dtt[lmaskt], sigma=errtt[lmaskt],\n", + " )\n", + "\n", + " tmod += np.sqrt(fit_model(ell_pivot, *param_tt) / cl0[ell_pivot]) * apo_patch\n", + " n_eff += apo_patch\n", + "\n", + " if ipix % 100 == 0:\n", + " print(f\" Patch {ipix}/{len(centers)}\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f247602f", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:27.018374Z", + "iopub.status.busy": "2026-03-24T23:42:27.018059Z", + "iopub.status.idle": "2026-03-24T23:42:27.856579Z", + "shell.execute_reply": "2026-03-24T23:42:27.855200Z" + } + }, + "outputs": [], + "source": [ + "tsm_raw = np.divide(tmod, n_eff, out=np.zeros_like(tmod), where=n_eff > 0)\n", + "tsm = hp.smoothing(tsm_raw, fwhm=np.radians(smoothing_fwhm_deg))\n", + "\n", + "hp.mollview(tsm,\n", + " title=rf\"Modulation map (patch radius={patch_radius_deg}$^\\circ$)\",\n", + " cmap=\"turbo\", norm=\"hist\")\n", + "hp.graticule(dpar=30)" + ] + }, + { + "cell_type": "markdown", + "id": "148fc907", + "metadata": {}, + "source": [ + "## 8. Filter the large-scale map\n", + "\n", + "We apply a low-pass filter $f_{\\rm LS}(\\ell) = \\sqrt{1 - w(\\ell)}$ to the original log-space map in harmonic space. This removes modes above the sigmoid transition, leaving only the large-scale structure that will be combined with the synthetic small-scale realization." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "931eb25b", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:27.865286Z", + "iopub.status.busy": "2026-03-24T23:42:27.864854Z", + "iopub.status.idle": "2026-03-24T23:42:27.919268Z", + "shell.execute_reply": "2026-03-24T23:42:27.918209Z" + } + }, + "outputs": [], + "source": [ + "np.random.seed(RANDOM_SEED)\n", + "\n", + "alm_ff_fullsky = hp.map2alm(log_he_ff, lmax=lmax, use_pixel_weights=True)\n", + "\n", + "# Low-pass filter: sqrt(1 - w)\n", + "w_he = sigmoid(ell_HE, x0=ell_fit_high, width=ell_fit_high / 10)\n", + "f_ls_he = np.sqrt(1.0 - w_he)\n", + "ii_LS_alm = hp.almxfl(alm_ff_fullsky, f_ls_he)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6d14a43", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:27.922389Z", + "iopub.status.busy": "2026-03-24T23:42:27.922117Z", + "iopub.status.idle": "2026-03-24T23:42:28.669077Z", + "shell.execute_reply": "2026-03-24T23:42:28.667912Z" + } + }, + "outputs": [], + "source": [ + "plt.loglog(hp.alm2cl(alm_ff_fullsky), color=\"gold\", label=\"Full-sky\")\n", + "plt.loglog(hp.alm2cl(ii_LS_alm), color=\"teal\", label=\"Low-pass filtered\")\n", + "plt.axvline(ell_fit_high, color=\"grey\", ls=\"--\", label=rf\"$\\ell={ell_fit_high}$\")\n", + "plt.xlabel(r\"Multipole $\\ell$\")\n", + "plt.ylabel(r\"$C_\\ell$\")\n", + "plt.legend()\n", + "plt.grid()" + ] + }, + { + "cell_type": "markdown", + "id": "77260c9c", + "metadata": {}, + "source": [ + "## 9. Generate small-scale realization and combine\n", + "\n", + "We generate a Gaussian realization from the extrapolated power-law $C_\\ell$ at the output resolution (nside=$N_{\\rm side,final}$), apply the high-pass sigmoid filter $f_{\\rm SS} = \\sqrt{w}$ to keep only the small-scale modes, and modulate the result with the spatially-varying modulation map.\n", + "\n", + "The final combined log-space map is:\n", + "\n", + "$$\n", + "i_{\\rm out}(\\hat{n}) = \\underbrace{f_{\\rm LS}(\\ell) \\ast i_{\\rm input}(\\hat{n})}_{\\text{large scales}} \\;+\\; \\underbrace{M(\\hat{n}) \\cdot \\bigl[f_{\\rm SS}(\\ell) \\ast s(\\hat{n})\\bigr]}_{\\text{modulated small scales}}\n", + "$$\n", + "\n", + "where $s(\\hat{n})$ is the synfast realization and $M(\\hat{n})$ is the modulation map normalized to unit RMS." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ccd1df0", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:42:28.672916Z", + "iopub.status.busy": "2026-03-24T23:42:28.672628Z", + "iopub.status.idle": "2026-03-24T23:44:05.619195Z", + "shell.execute_reply": "2026-03-24T23:44:05.617986Z" + } + }, + "outputs": [], + "source": [ + "# Generate small-scale realization from extrapolated power law\n", + "map_ss_full = hp.synfast(cl_pl, lmax=lmax_final, new=True, nside=nside_final)\n", + "\n", + "# High-pass filter: keep only small scales\n", + "alm_ss_full = hp.map2alm(map_ss_full, lmax=lmax_final)\n", + "alm_ss = hp.almxfl(alm_ss_full, f_ss_final)\n", + "map_ss = hp.alm2map(alm_ss, nside=nside_final)\n", + "\n", + "# Modulate with spatially-varying modulation map (normalized to unit RMS)\n", + "map_ss_modulated = map_ss * hp.ud_grade(\n", + " tsm / np.sqrt(np.mean(tsm ** 2)), nside_out=nside_final,\n", + ")\n", + "\n", + "# Large-scale map at output resolution\n", + "map_ls = hp.alm2map(ii_LS_alm, nside=nside_final)\n", + "\n", + "# Combine\n", + "ii_map_out = map_ss_modulated + map_ls" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2359ed1d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:44:05.622832Z", + "iopub.status.busy": "2026-03-24T23:44:05.622530Z", + "iopub.status.idle": "2026-03-24T23:44:06.632648Z", + "shell.execute_reply": "2026-03-24T23:44:06.631305Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 5))\n", + "hp.mollview(log_he_ff, cmap=\"turbo\", norm=\"hist\",\n", + " title=f\"Original (log-space, nside={nside})\", sub=(1, 2, 1))\n", + "hp.mollview(ii_map_out, cmap=\"turbo\", norm=\"hist\",\n", + " title=f\"With small scales (log-space, nside={nside_final})\", sub=(1, 2, 2))" + ] + }, + { + "cell_type": "markdown", + "id": "b7ea33b4", + "metadata": {}, + "source": [ + "## 10. Diagnostic visualizations\n", + "\n", + "### Gnomonic projections at different processing stages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06e1ea79", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:44:06.641875Z", + "iopub.status.busy": "2026-03-24T23:44:06.641515Z", + "iopub.status.idle": "2026-03-24T23:44:10.433443Z", + "shell.execute_reply": "2026-03-24T23:44:10.432480Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(15, 5))\n", + "hp.gnomview(ii_map_out, reso=3.75, xsize=320, rot=[30, -58],\n", + " sub=155, cmap=\"turbo\", norm=\"hist\", title=\"Combined $i$\")\n", + "hp.gnomview(map_ls, reso=3.75, xsize=320, rot=[30, -58],\n", + " sub=154, cmap=\"turbo\", norm=\"hist\", title=\"Large-scale $i$\")\n", + "hp.gnomview(map_ss, reso=3.75, xsize=320, rot=[30, -58],\n", + " sub=151, cmap=\"turbo\", norm=\"hist\", title=\"Small-scale $i$\")\n", + "hp.gnomview(map_ss_modulated, reso=3.75, xsize=320, rot=[30, -58],\n", + " sub=153, cmap=\"turbo\", norm=\"hist\", title=\"Modulated SS $i$\")\n", + "hp.gnomview(tsm, reso=3.75, xsize=320, rot=[30, -58],\n", + " sub=152, cmap=\"turbo\", norm=\"hist\", title=\"Modulation\")\n", + "\n", + "lon, lat = 130, 48\n", + "plt.figure(figsize=(15, 5))\n", + "hp.gnomview(ii_map_out, reso=3.75, xsize=320, rot=[lon, lat],\n", + " sub=155, cmap=\"turbo\", norm=\"hist\", title=\"Combined $i$\")\n", + "hp.gnomview(map_ls, reso=3.75, xsize=320, rot=[lon, lat],\n", + " sub=154, cmap=\"turbo\", norm=\"hist\", title=\"Large-scale $i$\")\n", + "hp.gnomview(map_ss, reso=3.75, xsize=320, rot=[lon, lat],\n", + " sub=151, cmap=\"turbo\", norm=\"hist\", title=\"Small-scale $i$\")\n", + "hp.gnomview(map_ss_modulated, reso=3.75, xsize=320, rot=[lon, lat],\n", + " sub=153, cmap=\"turbo\", norm=\"hist\", title=\"Modulated SS $i$\")\n", + "hp.gnomview(tsm, reso=3.75, xsize=320, rot=[lon, lat],\n", + " sub=152, cmap=\"turbo\", norm=\"hist\", title=\"Modulation\")" + ] + }, + { + "cell_type": "markdown", + "id": "487943a8", + "metadata": {}, + "source": [ + "### Power spectra at different processing stages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3cbce4bd", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:44:10.446148Z", + "iopub.status.busy": "2026-03-24T23:44:10.445717Z", + "iopub.status.idle": "2026-03-24T23:44:47.860556Z", + "shell.execute_reply": "2026-03-24T23:44:47.859368Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(10, 7))\n", + "\n", + "# Input small-scale Cl (fitted)\n", + "plt.loglog(cl_ss * cl_norm_final, color=\"green\", label=r\"Small-scale $C_\\ell$ (fitted)\")\n", + "\n", + "# Synfast-generated small-scale Cl\n", + "ell_map_ss, cl_norm_map_ss, cl_map_ss, dl_map_ss = run_anafast(map_ss, lmax_final)\n", + "plt.loglog(ell_map_ss, cl_map_ss * cl_norm_map_ss, color=\"orangered\",\n", + " label=r\"Small-scale $C_\\ell$ (synfast)\")\n", + "\n", + "# Modulated small-scale Cl\n", + "ell_map_ss_mod, cl_norm_map_ss_mod, cl_map_ss_mod, _ = run_anafast(map_ss_modulated, lmax_final)\n", + "plt.loglog(ell_map_ss_mod, cl_map_ss_mod * cl_norm_map_ss_mod, color=\"blue\",\n", + " label=r\"Small-scale $C_\\ell$ (modulated)\")\n", + "\n", + "# Large-scale Cl\n", + "plt.loglog(ell_HE * (ell_HE + 1) / (2 * np.pi) * hp.alm2cl(ii_LS_alm),\n", + " color=\"magenta\", label=r\"Large-scale $C_\\ell$ (filtered)\")\n", + "\n", + "# Total output Cl\n", + "ell_map_tot, cl_norm_map_tot, cl_map_tot, _ = run_anafast(ii_map_out, lmax_final)\n", + "plt.loglog(ell_map_tot, cl_map_tot * cl_norm_map_tot, color=\"gold\",\n", + " label=r\"Output $C_\\ell$\")\n", + "\n", + "# Input HE Cl\n", + "plt.loglog(ell_HE, ell_HE * (ell_HE + 1) / (2 * np.pi) * cl_HE,\n", + " color=\"black\", label=\"Input (HE)\")\n", + "\n", + "plt.legend()\n", + "plt.xlabel(r\"$\\ell$\")\n", + "plt.ylabel(r\"$\\mathcal{D}_\\ell$\")\n", + "plt.title(\"Angular power spectra comparison\")\n", + "plt.grid(True, linestyle=\"--\", alpha=0.6)\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "4a7efad0", + "metadata": {}, + "source": [ + "### Log-space power spectrum: input vs output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e77c3e53", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:44:47.864539Z", + "iopub.status.busy": "2026-03-24T23:44:47.864261Z", + "iopub.status.idle": "2026-03-24T23:45:00.067467Z", + "shell.execute_reply": "2026-03-24T23:45:00.066253Z" + } + }, + "outputs": [], + "source": [ + "ell_out, cl_norm_out, cl_out, dl_out = run_anafast(ii_map_out, lmax_final)\n", + "\n", + "plt.figure(figsize=(10, 5))\n", + "plt.loglog(ell_out, dl_out, label=\"Output (small-scales added)\")\n", + "plt.loglog(ell_HE, dl_HE, color=\"limegreen\", label=\"Input (HE)\")\n", + "plt.axvline(ell_fit_high, ls=\"--\", color=\"gray\", label=rf\"$\\ell={ell_fit_high}$\")\n", + "plt.axvline(ell_fit_low, ls=\"--\", color=\"gray\", label=rf\"$\\ell={ell_fit_low}$\")\n", + "plt.legend()\n", + "plt.grid()\n", + "plt.title(rf\"Log-space TT spectrum, $\\gamma_{{\\rm target}}={gamma_target:.2f}$\")\n", + "plt.ylabel(r\"$\\ell(\\ell+1)C_\\ell/2\\pi$\")\n", + "plt.xlabel(r\"$\\ell$\")" + ] + }, + { + "cell_type": "markdown", + "id": "5a3b4810", + "metadata": {}, + "source": [ + "## 11. Transform back to linear space\n", + "\n", + "Apply the inverse log transform $I = e^{i} - \\epsilon$ to recover the physical brightness temperature map." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1cf4d18", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:45:00.070487Z", + "iopub.status.busy": "2026-03-24T23:45:00.070218Z", + "iopub.status.idle": "2026-03-24T23:45:01.947766Z", + "shell.execute_reply": "2026-03-24T23:45:01.946519Z" + } + }, + "outputs": [], + "source": [ + "output_map = log_to_map(ii_map_out, eps)\n", + "\n", + "hp.mollview(output_map,\n", + " title=\"Final Free-Free Intensity at 30 GHz (with small scales)\",\n", + " cmap=\"turbo\", unit=r\"$T\\;[\\mu\\mathrm{K_{RJ}}]$\", norm=\"hist\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6d5bd6d", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:45:01.956440Z", + "iopub.status.busy": "2026-03-24T23:45:01.956032Z", + "iopub.status.idle": "2026-03-24T23:45:02.865458Z", + "shell.execute_reply": "2026-03-24T23:45:02.864306Z" + } + }, + "outputs": [], + "source": [ + "plt.figure(figsize=(12, 5))\n", + "hp.mollview(HE_ff_T, cmap=\"turbo\", norm=\"hist\",\n", + " title=f\"Original (nside={nside})\", sub=(1, 2, 1),\n", + " unit=r\"$T\\;[\\mu\\mathrm{K_{RJ}}]$\")\n", + "hp.mollview(output_map, cmap=\"turbo\", norm=\"hist\",\n", + " title=f\"With small scales (nside={nside_final})\", sub=(1, 2, 2),\n", + " unit=r\"$T\\;[\\mu\\mathrm{K_{RJ}}]$\")" + ] + }, + { + "cell_type": "markdown", + "id": "e2af8511", + "metadata": {}, + "source": [ + "### Gnomonic comparison at different stages" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d2300ed", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:45:02.876583Z", + "iopub.status.busy": "2026-03-24T23:45:02.876191Z", + "iopub.status.idle": "2026-03-24T23:45:05.504288Z", + "shell.execute_reply": "2026-03-24T23:45:05.502740Z" + } + }, + "outputs": [], + "source": [ + "lat = 15\n", + "plt.figure(figsize=(15, 10))\n", + "\n", + "hp.gnomview(HE_ff_T, title=\"$I$ (no small scales)\",\n", + " rot=[0, lat], reso=3.75, xsize=320, norm=\"hist\",\n", + " sub=231, cmap=\"turbo\")\n", + "hp.gnomview(log_he_ff, title=\"$i$ (no small scales)\",\n", + " rot=[0, lat], reso=3.75, xsize=320, norm=\"hist\",\n", + " sub=232, cmap=\"turbo\")\n", + "hp.gnomview(map_ss, title=\"Small scales\",\n", + " rot=[0, lat], reso=3.75, xsize=320, norm=\"hist\",\n", + " sub=233, cmap=\"turbo\")\n", + "hp.gnomview(map_ss_modulated, title=\"Small scales (modulated)\",\n", + " rot=[0, lat], reso=3.75, xsize=320, norm=\"hist\",\n", + " sub=234, cmap=\"turbo\")\n", + "hp.gnomview(ii_map_out, title=\"$i$ with small scales\",\n", + " rot=[0, lat], reso=3.75, xsize=320, norm=\"hist\",\n", + " sub=235, cmap=\"turbo\")\n", + "hp.gnomview(output_map, title=\"$I$ with small scales\",\n", + " rot=[0, lat], reso=3.75, xsize=320, norm=\"hist\",\n", + " sub=236, cmap=\"turbo\")" + ] + }, + { + "cell_type": "markdown", + "id": "65306c4a", + "metadata": {}, + "source": [ + "### Linear-space power spectrum: input vs output" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4dbbdaa5", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:45:05.525595Z", + "iopub.status.busy": "2026-03-24T23:45:05.525308Z", + "iopub.status.idle": "2026-03-24T23:45:18.252343Z", + "shell.execute_reply": "2026-03-24T23:45:18.251473Z" + } + }, + "outputs": [], + "source": [ + "ell_lin, _, cl_lin_out, _ = run_anafast(output_map, lmax_final)\n", + "ell_lin_in, _, cl_lin_in, _ = run_anafast(HE_ff_T, lmax)\n", + "\n", + "plt.figure(figsize=(10, 5))\n", + "plt.loglog(ell_lin_in, ell_lin_in * (ell_lin_in + 1) / (2 * np.pi) * cl_lin_in,\n", + " label=\"Input (HE)\")\n", + "plt.loglog(ell_lin, ell_lin * (ell_lin + 1) / (2 * np.pi) * cl_lin_out,\n", + " label=\"Output (with small scales)\")\n", + "plt.axvline(ell_fit_high, ls=\"--\", color=\"gray\", label=rf\"$\\ell={ell_fit_high}$\")\n", + "plt.axvline(ell_fit_low, ls=\"--\", color=\"gray\", label=rf\"$\\ell={ell_fit_low}$\")\n", + "plt.legend()\n", + "plt.grid()\n", + "plt.ylabel(r\"$\\ell(\\ell+1)C_\\ell/2\\pi\\;[\\mu\\mathrm{K_{RJ}}^2]$\")\n", + "plt.xlabel(r\"$\\ell$\")\n", + "plt.xlim(2, lmax_final + 1000)" + ] + }, + { + "cell_type": "markdown", + "id": "mvudnlgjb3", + "source": "## 12. Comparison with the f1 model\n\nThe existing PySM `f1` free-free model uses a power-law frequency scaling with constant spectral index $\\alpha = -2.14$ and a template at 30 GHz derived from PySM 2 (nside=512, no small-scale injection). Here we compare the new `f2` model (this notebook) with `f1` in both power spectrum and map space.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "vdypbslipk", + "source": "# Load f1 model at the output nside\nsky_f1 = pysm3.Sky(nside=nside_final, preset_strings=[\"f1\"])\nf1_map = sky_f1.get_emission(30 * u.GHz)[0].value # I component in uK_RJ\n\n# Downgrade f2 output to same nside for fair comparison\nf2_map = output_map if hp.get_nside(output_map) == nside_final else hp.ud_grade(output_map, nside_out=nside_final)\n\nprint(f\"f1 map: nside={hp.get_nside(f1_map)}, range=[{f1_map.min():.2f}, {f1_map.max():.2f}] uK_RJ\")\nprint(f\"f2 map: nside={hp.get_nside(f2_map)}, range=[{f2_map.min():.2f}, {f2_map.max():.2f}] uK_RJ\")", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "enefvnu5z8a", + "source": "### Full-sky power spectrum: f1 vs f2", + "metadata": {} + }, + { + "cell_type": "code", + "id": "1lr5pwuk03t", + "source": "lmax_comp = min(3 * nside_final - 1, 6000)\n\nell_f1, _, cl_f1, _ = run_anafast(f1_map, lmax_comp)\nell_f2, _, cl_f2, _ = run_anafast(f2_map, lmax_comp)\n\nplt.figure(figsize=(10, 5))\nplt.loglog(ell_f1, ell_f1 * (ell_f1 + 1) / (2 * np.pi) * cl_f1,\n label=\"f1 (PySM 2 template)\", color=\"C0\")\nplt.loglog(ell_f2, ell_f2 * (ell_f2 + 1) / (2 * np.pi) * cl_f2,\n label=\"f2 (this notebook)\", color=\"C1\")\nplt.axvline(ell_fit_high, ls=\"--\", color=\"gray\", label=rf\"$\\ell={ell_fit_high}$\")\nplt.legend()\nplt.grid()\nplt.title(\"Full-sky $\\\\mathcal{D}_\\\\ell$ comparison: f1 vs f2 at 30 GHz\")\nplt.ylabel(r\"$\\ell(\\ell+1)C_\\ell/2\\pi\\;[\\mu\\mathrm{K_{RJ}}^2]$\")\nplt.xlabel(r\"$\\ell$\")\nplt.xlim(2, lmax_comp)", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "y1itm009pe", + "source": "### Masked-sky power spectrum: f1 vs f2\n\nWe apply the Planck GAL080 Galactic plane mask (retaining ~80% of the sky, excluding the brightest Galactic emission) to compare the spectra away from the Galactic plane.", + "metadata": {} + }, + { + "cell_type": "code", + "id": "elolch2kpdr", + "source": "# Download Planck Galactic plane mask\nplanck_mask_filename = DATA_DIR / \"HFI_Mask_GalPlane-apo2_2048_R2.00.fits\"\nplanck_mask_url = \"http://pla.esac.esa.int/pla/aio/product-action?MAP.MAP_ID=HFI_Mask_GalPlane-apo2_2048_R2.00.fits\"\n\nif not planck_mask_filename.exists():\n print(f\"Downloading Planck mask to {planck_mask_filename} ...\")\n urllib.request.urlretrieve(planck_mask_url, planck_mask_filename)\nelse:\n print(f\"Using existing file {planck_mask_filename}\")\n\n# Read GAL080 mask and upgrade to output nside\nplanck_mask_2048 = hp.read_map(planck_mask_filename, [\"GAL080\"])\nplanck_mask = hp.ud_grade(planck_mask_2048, nside_out=nside_final)\nplanck_mask = np.int_(np.ma.masked_not_equal(planck_mask, 0.0).mask)\nfsky_mask = planck_mask.sum() / planck_mask.size\nprint(f\"GAL080 mask f_sky = {fsky_mask:.3f}\")\n\nell_f1m, _, cl_f1m, _ = run_anafast(f1_map * planck_mask, lmax_comp)\nell_f2m, _, cl_f2m, _ = run_anafast(f2_map * planck_mask, lmax_comp)\n\nplt.figure(figsize=(10, 5))\nplt.loglog(ell_f1m, ell_f1m * (ell_f1m + 1) / (2 * np.pi) * cl_f1m / fsky_mask,\n label=\"f1 (PySM 2 template)\", color=\"C0\")\nplt.loglog(ell_f2m, ell_f2m * (ell_f2m + 1) / (2 * np.pi) * cl_f2m / fsky_mask,\n label=\"f2 (this notebook)\", color=\"C1\")\nplt.axvline(ell_fit_high, ls=\"--\", color=\"gray\", label=rf\"$\\ell={ell_fit_high}$\")\nplt.legend()\nplt.grid()\nplt.title(r\"Masked-sky (Planck GAL080) $\\mathcal{D}_\\ell$ comparison: f1 vs f2 at 30 GHz\")\nplt.ylabel(r\"$\\ell(\\ell+1)C_\\ell/2\\pi\\;[\\mu\\mathrm{K_{RJ}}^2]$\")\nplt.xlabel(r\"$\\ell$\")\nplt.xlim(2, lmax_comp)", + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "markdown", + "id": "5puhtl3wyh", + "source": "### Map comparison: f1 vs f2 (gnomonic projections)", + "metadata": {} + }, + { + "cell_type": "code", + "id": "bap4ih6sy8n", + "source": [ + "# Remove monopole globally for fair comparison\n", + "f1_map_nomon = f1_map - np.mean(f1_map)\n", + "f2_map_nomon = f2_map - np.mean(f2_map)\n", + "f1_at_f2_nside = hp.ud_grade(f1_map_nomon, nside_out=hp.get_nside(f2_map_nomon))\n", + "\n", + "# Determine common color range\n", + "vmin_gc = min(np.percentile(f1_at_f2_nside, 1), np.percentile(f2_map_nomon, 1))\n", + "vmax_gc = max(np.percentile(f1_at_f2_nside, 99), np.percentile(f2_map_nomon, 99))\n", + "\n", + "# Galactic plane region\n", + "plt.figure(figsize=(12, 5))\n", + "hp.gnomview(f1_at_f2_nside, reso=3.75, xsize=400, rot=[0, 0],\n", + " sub=121, cmap=\"turbo\", title=\"f1 (Galactic center)\",\n", + " min=vmin_gc, max=vmax_gc)\n", + "hp.gnomview(f2_map_nomon, reso=3.75, xsize=400, rot=[0, 0],\n", + " sub=122, cmap=\"turbo\", title=\"f2 (Galactic center)\",\n", + " min=vmin_gc, max=vmax_gc)\n", + "\n", + "# High-latitude region\n", + "vmin_hl = min(np.percentile(f1_at_f2_nside, 5), np.percentile(f2_map_nomon, 5))\n", + "vmax_hl = max(np.percentile(f1_at_f2_nside, 95), np.percentile(f2_map_nomon, 95))\n", + "\n", + "plt.figure(figsize=(12, 5))\n", + "hp.gnomview(f1_at_f2_nside, reso=3.75, xsize=400, rot=[0, 60],\n", + " sub=121, cmap=\"turbo\", title=\"f1 (high latitude)\",\n", + " min=vmin_hl, max=vmax_hl)\n", + "hp.gnomview(f2_map_nomon, reso=3.75, xsize=400, rot=[0, 60],\n", + " sub=122, cmap=\"turbo\", title=\"f2 (high latitude)\",\n", + " min=vmin_hl, max=vmax_hl)\n", + "\n", + "# Intermediate latitude\n", + "plt.figure(figsize=(12, 5))\n", + "hp.gnomview(f1_at_f2_nside, reso=3.75, xsize=400, rot=[130, 30],\n", + " sub=121, cmap=\"turbo\", title=\"f1 (mid latitude)\",\n", + " min=vmin_hl, max=vmax_hl)\n", + "hp.gnomview(f2_map_nomon, reso=3.75, xsize=400, rot=[130, 30],\n", + " sub=122, cmap=\"turbo\", title=\"f2 (mid latitude)\",\n", + " min=vmin_hl, max=vmax_hl)\n" + ], + "metadata": {}, + "execution_count": null, + "outputs": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Full-sky mollview comparison (monopole removed, same color scale)\n", + "vmin_mol = min(np.percentile(f1_at_f2_nside, 1), np.percentile(f2_map_nomon, 1))\n", + "vmax_mol = max(np.percentile(f1_at_f2_nside, 99), np.percentile(f2_map_nomon, 99))\n", + "\n", + "plt.figure(figsize=(12, 5))\n", + "hp.mollview(f1_at_f2_nside, cmap=\"turbo\",\n", + " title=\"f1 (monopole removed)\", sub=(1, 2, 1),\n", + " unit=r\"$T\\;[\\mu\\mathrm{K_{RJ}}]$\",\n", + " min=vmin_mol, max=vmax_mol)\n", + "hp.mollview(f2_map_nomon, cmap=\"turbo\",\n", + " title=\"f2 (monopole removed)\", sub=(1, 2, 2),\n", + " unit=r\"$T\\;[\\mu\\mathrm{K_{RJ}}]$\",\n", + " min=vmin_mol, max=vmax_mol)\n" + ] + }, + { + "cell_type": "markdown", + "id": "bd6cd343", + "metadata": {}, + "source": "## 13. Save output FITS files\n\nWe save three files:\n1. The high-resolution free-free intensity map at 30 GHz (with small scales injected)\n2. The modulation map (dimensionless, at nside=128)\n3. The original HE free-free map at 30 GHz (nside=256) for reference" + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ecd84cd", + "metadata": { + "execution": { + "iopub.execute_input": "2026-03-24T23:45:18.255601Z", + "iopub.status.busy": "2026-03-24T23:45:18.255302Z", + "iopub.status.idle": "2026-03-24T23:45:20.735364Z", + "shell.execute_reply": "2026-03-24T23:45:20.734007Z" + } + }, + "outputs": [], + "source": [ + "# Output free-free map with small-scales\n", + "fname = DATA_DIR / f\"ff_highres_map_nside{nside_final}_f2_30GHz.fits\"\n", + "hp.write_map(fname, output_map, overwrite=True, dtype=np.float64,\n", + " column_units=\"uK_RJ\")\n", + "add_metadata([str(fname)], coord=\"G\", unit=\"uK_RJ\")\n", + "output_files.append(fname)\n", + "\n", + "# Modulation map (dimensionless)\n", + "nside_mod = hp.get_nside(tsm)\n", + "fname = DATA_DIR / f\"ff_modulation_map_nside{nside_mod}_f2_30GHz.fits\"\n", + "hp.write_map(fname, tsm, overwrite=True, dtype=np.float64,\n", + " column_units=\"dimensionless\")\n", + "add_metadata([str(fname)], coord=\"G\", unit=\"dimensionless\")\n", + "output_files.append(fname)\n", + "\n", + "# Input HE free-free map (without small-scales)\n", + "fname = DATA_DIR / f\"ff_HE_map_nside{nside}_30GHz.fits\"\n", + "hp.write_map(fname, HE_ff_T, overwrite=True, dtype=np.float64,\n", + " column_units=\"uK_RJ\")\n", + "add_metadata([str(fname)], coord=\"G\", unit=\"uK_RJ\")\n", + "output_files.append(fname)\n", + "\n", + "print(\"Output files:\")\n", + "for f in output_files:\n", + " print(f\" {f}\")" + ] + } + ], + "metadata": { + "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.9" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} \ No newline at end of file