diff --git a/docs/low_frequency_synchrotron_model.ipynb b/docs/low_frequency_synchrotron_model.ipynb
new file mode 100644
index 00000000..d868003d
--- /dev/null
+++ b/docs/low_frequency_synchrotron_model.ipynb
@@ -0,0 +1,731 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "9e63c846",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "import healpy as hp\n",
+    "from scipy.odr import Model, Data, ODR\n",
+    "from pygdsm import gsm16\n",
+    "import pysm3 as pysm\n",
+    "import pysm3.units as u\n",
+    "import cmasher as cmr\n",
+    "import wget\n",
+    "from pathlib import Path\n",
+    "import zipfile \n",
+    "import gdown"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "906175d8",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def sform(param, xxx):\n",
+    "    '''linear reggression'''\n",
+    "    return xxx * param[0] + param[1]\n",
+    "\n",
+    "def remove_off(map, hasmap, maperr, haserr):\n",
+    "    nside_superpix=8\n",
+    "\n",
+    "    super_pixs=  np.arange(hp.nside2npix(nside_superpix )) #added\n",
+    "    offsets = np.zeros_like(super_pixs)*1. \n",
+    "    offmap = np.full(hp.nside2npix(nside_superpix), np.nan)  \n",
+    "   \n",
+    "    nansize =0 \n",
+    "    for jj, ipix in enumerate(super_pixs) : \n",
+    "\n",
+    "        super_map = np.zeros(hp.nside2npix(nside_superpix))\n",
+    "        super_map [ipix] =1\n",
+    "        patch =  (hp.ud_grade(super_map, nside_out=hp.get_nside(map))) .astype('bool')\n",
+    "     \n",
+    "        try:\n",
+    "\n",
+    "            errx = haserr[patch] * hasmap[patch]\n",
+    "            erry =  maperr * map[patch]\n",
+    "            mydata = Data(hasmap[patch], map[patch], wd=errx**-2, we=erry**-2)\n",
+    "            myodr = ODR(mydata, Model(sform), beta0=[0.0, 0.0])\n",
+    "            myoutput = myodr.run()\n",
+    "            cfit = myoutput.beta[1]\n",
+    "            if np.corrcoef(hasmap[patch], map[patch])[0,1] > 0.75:\n",
+    "                offsets[jj] = cfit\n",
+    "            else:\n",
+    "                offsets[jj] = np.nan\n",
+    "        except  ValueError: \n",
+    "            nansize+=1 \n",
+    "            offsets[jj]= np.nan \n",
+    "            \n",
+    "    offmap[super_pixs] = offsets \n",
+    "\n",
+    "    avec = np.nanmedian( offsets )\n",
+    "    \n",
+    "    return avec\n",
+    "\n",
+    "input_dir = Path(\"../input_data\")\n",
+    "beta_path = input_dir / \"synchrotron_spectral_parameters/\"                                 \n",
+    "  \n",
+    "if not (beta_path.exists()):                             \n",
+    "    input_dir.mkdir(parents=True, exist_ok=True)\n",
+    "    zip_path = input_dir / \"synchrotron_spectral_parameters.zip\"              \n",
+    "                                                                                \n",
+    "    gdown.download(\n",
+    "        id=\"11N8u-Vrd3stFNN_v6UirGjLbcYrJm5GC\",                               \n",
+    "        output=str(zip_path),                                                 \n",
+    "        quiet=False,\n",
+    "    )                                                                         \n",
+    "                  \n",
+    "    with zipfile.ZipFile(zip_path) as zf:                                     \n",
+    "        zf.extractall(input_dir)\n",
+    "\n",
+    "beta_path = input_dir / \"synchrotron_spectral_parameters\" / \"final_beta_1deg.npy\"\n",
+    "cs_path   = input_dir / \"synchrotron_spectral_parameters\" / \"final_cs_1deg.npy\"  \n",
+    "                                                                                \n",
+    "beta_map = np.load(beta_path)\n",
+    "cs_map   = np.load(cs_path)\n",
+    "\n",
+    "nside = 256\n",
+    "\n",
+    "# Show average spectral index value\n",
+    "move_nu = np.arange(45, 2300, 10)\n",
+    "move_beta = np.array([beta_map + cs_map * np.log(move_nu[xx] / 45) for xx in range(len(move_nu))])\n",
+    "plt.figure(figsize=(10, 5))\n",
+    "plt.plot(move_nu, np.mean(move_beta, axis=1))\n",
+    "plt.fill_between(move_nu, np.mean(move_beta, axis=1) - np.std(move_beta, axis=1), np.mean(move_beta, axis=1) + np.std(move_beta, axis=1), alpha=0.3)\n",
+    "plt.xlabel('Frequency (MHz)')\n",
+    "plt.ylabel('Average Spectral Index')\n",
+    "plt.title('Average Spectral Index vs Frequency')\n",
+    "plt.show()\n",
+    "\n",
+    "move_nu = np.arange(0.045, 70, 0.10)\n",
+    "move_beta = np.array([beta_map + cs_map * np.log(move_nu[xx] / 0.045) for xx in range(len(move_nu))])\n",
+    "plt.figure(figsize=(10, 5))\n",
+    "plt.plot(move_nu, np.mean(move_beta, axis=1))\n",
+    "plt.fill_between(move_nu, np.mean(move_beta, axis=1) - np.std(move_beta, axis=1), np.mean(move_beta, axis=1) + np.std(move_beta, axis=1), alpha=0.3)\n",
+    "plt.xlabel('Frequency (GHz)')\n",
+    "plt.ylabel('Average Spectral Index')\n",
+    "plt.title('Average Spectral Index vs Frequency')\n",
+    "plt.show()\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "fa5bae32",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def rotate_map(mapin, coord_in , coord_out ) : \n",
+    "    alm = hp.map2alm (mapin )\n",
+    "    R=  hp.Rotator(coord=[coord_in, coord_out ])\n",
+    "    alm =  R.rotate_alm (alm )\n",
+    "    map_out = hp.alm2map (alm , nside= hp.get_nside(mapin ))\n",
+    "    return map_out\n",
+    "\n",
+    "#empirical data set to nside 256\n",
+    "fstring=\"haslam408_ds_Remazeilles2014.fits\"\n",
+    "\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    hasmap =  hp.read_map(f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/foregrounds/haslam_2014/{fstring}\", out =\"../input_data\")\n",
+    "    hasmap =  hp.read_map(f\"../input_data/{fstring}\"  )\n",
+    "\n",
+    "hasmap = hp.ud_grade(hasmap, nside)\n",
+    "hasmap = hasmap - 8.9\n",
+    "\n",
+    "#read in the ff maps\n",
+    "fstring=\"COM_CompMap_freefree-commander_0256_R2.00.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    tmp =  hp.read_map(f\"../input_data/{fstring}\" , field= ['EM_ML', 'TEMP_ML']) \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/foregrounds/haslam_2014/{fstring}\", out =\"../input_data\")\n",
+    "    tmp =  hp.read_map(f\"../input_data/{fstring}\", field= ['EM_ML', 'TEMP_ML'])\n",
+    "planck_te = tmp [1] #in K\n",
+    "\n",
+    "fstring = \"EM_mean_std.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    hust=  hp.read_map(filename=f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://zenodo.org/records/10523170/files/{fstring}\" , out =\"../input_data\")\n",
+    "    hust =  hp.read_map(f\"../input_data/{fstring}\")\n",
+    "\n",
+    "Tff = lambda Te , nu,EM  : Te * (1.0 - np.exp(-tauff(Te,nu, EM ))) \n",
+    "tauff = lambda Te,nu , EM : 0.05468 *   (Te)**-1.5 * (nu)**-2 * EM *gff (Te,nu) \n",
+    "Zi =1 \n",
+    "gff =lambda Te ,nu :  np.log (np.exp(1.0) +np.exp (5.960 -np.sqrt(3)/np.pi *np.log (Zi *nu *(Te/(1e4)  )**-1.5  )) )  \n",
+    "\n",
+    "pysm_sync_mod = 's7'"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "adc60742",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# eda map is at 159 Mhz and 3.1 degree\n",
+    "goal_f = 159.0\n",
+    "fstring = \"EDA2_159MHz_I_noPrior_HPXbin.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    edamap =  hp.read_map(f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/foregrounds/EDA2/{fstring}\" , out =\"../input_data\")\n",
+    "    edamap =  hp.read_map(f\"../input_data/{fstring}\")\n",
+    "eda2 = hp.ud_grade(edamap, nside)\n",
+    "eda2 = rotate_map(eda2, 'C', 'G')\n",
+    "thevec = hp.pixelfunc.ang2vec(np.pi/2.7, np.pi/1.5, lonlat=False)\n",
+    "centre = hp.query_disc(nside, thevec, 0.7)\n",
+    "eda2[centre] = np.nan\n",
+    "\n",
+    "eda2[np.where(np.isnan(eda2))[0]] = hp.UNSEEN\n",
+    "eda2_5 = hp.smoothing(eda2, fwhm=np.radians(np.sqrt((5.0)**2 - (3.1)**2)))\n",
+    "has_5 = hp.smoothing(hasmap, fwhm=np.radians(np.sqrt((5.0)**2 - (56./60.)**2)))\n",
+    "haserr = np.sqrt((0.1 * hasmap)**2 + 1.3**2) /  hasmap\n",
+    "eda2_5[np.where(eda2_5 == hp.UNSEEN)] = np.nan\n",
+    "eda2[np.where(eda2 == hp.UNSEEN)] = np.nan\n",
+    "eda2_off = remove_off(eda2_5, has_5, 0.05, haserr)\n",
+    "\n",
+    "hp.mollview(eda2 - eda2_off, min=130, max=1000)\n",
+    "plt.title('EDA2')\n",
+    "plt.show()\n",
+    "\n",
+    "beta_159 = (beta_map + cs_map * np.log(goal_f/45.))\n",
+    "ff159 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "param_159 = hasmap * (goal_f/408.)**beta_159 + hp.smoothing(ff159, fwhm=np.radians(np.sqrt((3.1)**2 - 1.0**2)))\n",
+    "param_159[np.where(np.isnan(param_159))[0]] = hp.UNSEEN\n",
+    "param_159 = hp.smoothing(param_159, fwhm=np.radians(np.sqrt((3.1)**2 - (56./60.)**2)))\n",
+    "param_159[np.where(param_159 == hp.UNSEEN)[0]] = np.nan\n",
+    "param_159[np.where(np.isnan(eda2))[0]] = np.nan\n",
+    "\n",
+    "hp.mollview(param_159,  min=130, max=1000)\n",
+    "plt.title('Parametric 159 model')\n",
+    "plt.show()\n",
+    "\n",
+    "sky = pysm.Sky(nside=nside, preset_strings=[pysm_sync_mod])\n",
+    "pysm159 = sky.get_emission((goal_f/1000.) * u.GHz)[0] * 1.e-6 #gives Stokes I, Q and U in microK\n",
+    "pysm159.value[np.where(np.isnan(eda2))[0]] = hp.UNSEEN\n",
+    "map_resol = hp.nside2resol(nside, arcmin=True) / 60.\n",
+    "pysm159 = hp.smoothing(pysm159.value, fwhm=np.radians(np.sqrt((3.1)**2 - map_resol**2))) + hp.smoothing(ff159, fwhm=np.radians(np.sqrt((3.1)**2 - 1.0**2)))\n",
+    "pysm159[np.where(pysm159 == hp.UNSEEN)] = np.nan\n",
+    "\n",
+    "hp.mollview(pysm159,  min=130, max=1000)\n",
+    "plt.title('PSM 159')\n",
+    "plt.show()\n",
+    "\n",
+    "# galactic plane mask\n",
+    "fstring = f\"HFI_Mask_GalPlane-apo0_2048_R2.00.fits\"\n",
+    "try : \n",
+    "    print(f\"reading {fstring} \")\n",
+    "    gmask = hp.read_map(f\"../input_data/{fstring}\" , field=4)\n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"http://pla.esac.esa.int/pla/aio/product-action?SOURCE_LIST.NAME={fstring}\", out =\"../input_data\") \n",
+    "    gmask = hp.read_map(f\"../input_data/{fstring}\" , field=4)  \n",
+    "\n",
+    "gmask = hp.ud_grade(gmask, nside)\n",
+    "gmask = np.where(gmask < 0.99, 0.0, 1.1)\n",
+    "gmask[np.where(gmask == 0.0)[0]] = np.nan   \n",
+    "\n",
+    "diff = eda2 - eda2_off - param_159\n",
+    "diff_pysm = eda2 - eda2_off - pysm159\n",
+    "\n",
+    "cmapEC = cmr.eclipse  \n",
+    "\n",
+    "fig, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1, figsize=(6,8))\n",
+    "plt.axes(ax1)\n",
+    "hp.mollview(diff / (eda2 - eda2_off),  min=-1.0, max=1.0, hold=True, title='Param 159 MHz', cmap=cmapEC)\n",
+    "plt.axes(ax2)\n",
+    "hp.mollview(diff_pysm / (eda2 - eda2_off),  min=-1.0, max=1.0, hold=True, title='PySM 159 MHz', cmap=cmapEC)\n",
+    "plt.show()\n",
+    "\n",
+    "scat_x = (eda2 - eda2_off) * gmask\n",
+    "notnans = np.where(~np.isnan(scat_x))[0]\n",
+    "\n",
+    "fig, ((ax1, ax2, ax3)) = plt.subplots(nrows=1, ncols=3, figsize=(24,6))\n",
+    "plt.axes(ax1)\n",
+    "plt.scatter(scat_x, pysm159* gmask, s=20, edgecolors='teal', facecolors='none', label='PYSM3 Model', alpha=0.4)\n",
+    "plt.scatter(scat_x, param_159* gmask, s=20, edgecolors='darkslateblue', facecolors='none', label='Parametric Model', alpha=0.4)\n",
+    "ax1.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax1.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlim([0.0, np.nanmax(pysm159*gmask)+200.])\n",
+    "plt.ylim([0.0, np.nanmax(pysm159*gmask)+200.])\n",
+    "plt.plot(scat_x, scat_x, 'k-')\n",
+    "plt.text(870, 2700, f'Parametric corr: {np.corrcoef(scat_x[notnans], (param_159* gmask)[notnans])[0,1]:.3f} K', color='darkslateblue', fontsize=14)\n",
+    "plt.text(870, 2400, f'PySM corr: {np.corrcoef(scat_x[notnans], (pysm159* gmask)[notnans])[0,1]:.3f} K', color='teal', fontsize=14)\n",
+    "plt.xlabel('Data 159 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.ylabel('Model 159 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.axes(ax2)\n",
+    "plt.hist(abs(diff/(eda2 - eda2_off))*100., bins=40, range=[0, 100], label='Parametric (Full)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(eda2 - eda2_off))*100., bins=40, range=[0, 100], label='PYSM3 (Full)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(eda2 - eda2_off))*100.), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(eda2 - eda2_off))*100.), color='teal', linestyle='--')\n",
+    "ax2.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax2.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12)\n",
+    "plt.axes(ax3)\n",
+    "plt.hist(abs(diff/(eda2 - eda2_off))*100.*gmask, bins=40, range=[0, 100], label='Parametric (NoGal)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(eda2 - eda2_off))*100.*gmask, bins=40, range=[0, 100], label='PYSM3 (NoGal)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(eda2 - eda2_off))*100.*gmask), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(eda2 - eda2_off))*100.*gmask), color='teal', linestyle='--')\n",
+    "ax3.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax3.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12)\n",
+    "plt.show()\n",
+    "\n",
+    "print(np.nanmedian(abs(diff/(eda2 - eda2_off))*100.*gmask), np.nanpercentile(abs(diff/(eda2 - eda2_off))*100.*gmask, 25), np.nanpercentile(abs(diff/(eda2 - eda2_off))*100.*gmask, 75))\n",
+    "print(np.nanmedian(abs(diff_pysm/(eda2 - eda2_off))*100.*gmask), np.nanpercentile(abs(diff_pysm/(eda2 - eda2_off))*100.*gmask, 25), np.nanpercentile(abs(diff_pysm/(eda2 - eda2_off))*100.*gmask, 75))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "63e663ba",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# dwingoloo is at 820 Mhz and 1.2 degree\n",
+    "goal_f = 820.0\n",
+    "fstring = \"Dwingeloo_Kelvins_1_256.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    dwing =  hp.read_map(filename=f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/foregrounds/dwingeloo_820/{fstring}\" , out =\"../input_data\")\n",
+    "    dwing =  hp.read_map(f\"../input_data/{fstring}\")\n",
+    "dwing[np.where(dwing < -1000)[0]] = np.nan\n",
+    "dwing[np.where(np.isnan(dwing))[0]] = hp.UNSEEN\n",
+    "dwing_5 = hp.smoothing(dwing, fwhm=np.radians(np.sqrt((5.0)**2 - (1.2)**2)))\n",
+    "dwing_5[np.where(dwing_5 == hp.UNSEEN)] = np.nan\n",
+    "dwing[np.where(dwing == hp.UNSEEN)] = np.nan\n",
+    "dwing_off = remove_off(dwing_5, has_5, 0.06, haserr)\n",
+    "\n",
+    "hp.mollview(dwing, norm='hist')\n",
+    "plt.title('Dwingeloo')\n",
+    "plt.show()\n",
+    "\n",
+    "beta_820 = (beta_map + cs_map * np.log(goal_f/45.))\n",
+    "ff820 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "param_820 = hasmap * (goal_f/408.)**beta_820 + hp.smoothing(ff820, fwhm=np.radians(np.sqrt((1.2)**2 - 1.0**2)))\n",
+    "param_820[np.where(np.isnan(param_820))[0]] = hp.UNSEEN\n",
+    "param_820 = hp.smoothing(param_820, fwhm=np.radians(np.sqrt((1.2)**2 - (56./60.)**2)))\n",
+    "param_820[np.where(param_820 == hp.UNSEEN)[0]] = np.nan\n",
+    "param_820[np.where(np.isnan(dwing))] = np.nan\n",
+    "\n",
+    "hp.mollview(param_820, norm='hist')\n",
+    "plt.title('Parametric 820 model')\n",
+    "plt.show()\n",
+    "\n",
+    "ff820 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "sky = pysm.Sky(nside=nside, preset_strings=[pysm_sync_mod])\n",
+    "pysm820 = sky.get_emission((goal_f/1000.) * u.GHz)[0] * 1.e-6 #gives Stokes I, Q and U in microK\n",
+    "pysm820.value[np.where(np.isnan(dwing))[0]] = hp.UNSEEN\n",
+    "map_resol = hp.nside2resol(nside, arcmin=True) / 60.\n",
+    "pysm820 = hp.smoothing(pysm820.value, fwhm=np.radians(np.sqrt((1.2)**2 - map_resol**2))) + hp.smoothing(ff820, fwhm=np.radians(np.sqrt((1.2)**2 - 1.0**2)))\n",
+    "pysm820[np.where(pysm820 == hp.UNSEEN)] = np.nan\n",
+    "\n",
+    "hp.mollview(pysm820, norm='hist')\n",
+    "plt.title('PYSM 820')\n",
+    "plt.show()\n",
+    "\n",
+    "diff = dwing - dwing_off - param_820\n",
+    "diff_pysm = dwing - dwing_off - pysm820\n",
+    "\n",
+    "fig, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1, figsize=(6,8))\n",
+    "plt.axes(ax1)\n",
+    "hp.mollview(diff / (dwing - dwing_off),  min=-1.0, max=1.0, hold=True, title='Param 820 MHz', cmap=cmapEC)\n",
+    "plt.axes(ax2)\n",
+    "hp.mollview(diff_pysm / (dwing - dwing_off ),  min=-1.0, max=1.0, hold=True, title='PySM 820 MHz', cmap=cmapEC)\n",
+    "plt.show()\n",
+    "\n",
+    "scat_x = (dwing - dwing_off) * gmask\n",
+    "notnans = np.where(~np.isnan(scat_x))[0]\n",
+    "\n",
+    "fig, ((ax1, ax2, ax3)) = plt.subplots(nrows=1, ncols=3, figsize=(24,6))\n",
+    "plt.axes(ax1)\n",
+    "plt.scatter(scat_x, pysm820*gmask, s=20, edgecolors='teal', facecolors='none', label='PYSM3 Model', alpha=0.4)\n",
+    "plt.scatter(scat_x, param_820*gmask, s=20, edgecolors='darkslateblue', facecolors='none', label='Parametric Model', alpha=0.4)\n",
+    "plt.plot(scat_x, scat_x, 'k-')\n",
+    "plt.xlim([0.0, 1+np.nanmax(pysm820*gmask)])\n",
+    "plt.ylim([0.0, 1+np.nanmax(pysm820*gmask)])\n",
+    "ax1.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax1.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.text(1, 10, f'Parametric corr: {np.corrcoef(scat_x[notnans], (param_820*gmask)[notnans])[0,1]:.3f} K', color='darkslateblue', fontsize=14)\n",
+    "plt.text(1, 8.5, f'PYSM corr: {np.corrcoef(scat_x[notnans], (pysm820*gmask)[notnans])[0,1]:.3f} K', color='teal', fontsize=14)\n",
+    "plt.xlabel('Data 820 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.ylabel('Model 820 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.axes(ax2)\n",
+    "plt.hist(abs(diff/(dwing - dwing_off))*100., bins=40, range=[0, 300], label='Parametric (Full)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(dwing - dwing_off))*100., bins=40, range=[0, 200], label='PYSM3 (Full)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(dwing - dwing_off))*100.), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(dwing - dwing_off))*100.), color='teal', linestyle='--')\n",
+    "ax2.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax2.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12)\n",
+    "plt.axes(ax3)\n",
+    "plt.hist(abs(diff/(dwing - dwing_off))*100.*gmask, bins=40, range=[0, 300], label='Parametric (No Gal)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(dwing - dwing_off))*100.*gmask, bins=40, range=[0, 200], label='PYSM3 (No Gal)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(dwing - dwing_off))*100.*gmask), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(dwing - dwing_off))*100.*gmask), color='teal', linestyle='--')\n",
+    "ax3.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax3.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12)\n",
+    "plt.show()\n",
+    "\n",
+    "print(np.nanmedian(abs(diff/(dwing - dwing_off))*100.*gmask), np.nanpercentile(abs(diff/(dwing - dwing_off))*100.*gmask, 25), np.nanpercentile(abs(diff/(dwing - dwing_off))*100.*gmask, 75))\n",
+    "print(np.nanmedian(abs(diff_pysm/(dwing - dwing_off))*100.*gmask), np.nanpercentile(abs(diff_pysm/(dwing - dwing_off))*100.*gmask, 25), np.nanpercentile(abs(diff_pysm/(dwing - dwing_off))*100.*gmask, 75))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "bcd678f0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# jonas is at 2326 Mhz and 20 arcmin\n",
+    "goal_f = 2326.0\n",
+    "fstring = \"lambda_23de_hea.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    jonas =  hp.read_map(filename=f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/foregrounds/rhodes_2326/{fstring}\" , out =\"../input_data\")\n",
+    "    jonas =  hp.read_map(f\"../input_data/{fstring}\")\n",
+    "jonas = rotate_map(jonas, 'C', 'G')\n",
+    "jonas[np.where(jonas < 0.03)[0]] = np.nan\n",
+    "jonas[np.where(np.isnan(jonas))[0]] = hp.UNSEEN\n",
+    "jonas_5 = hp.smoothing(jonas, fwhm=np.radians(np.sqrt((5.0)**2 - (20./60.)**2)))\n",
+    "jonas_5[np.where(jonas_5 == hp.UNSEEN)] = np.nan\n",
+    "jonas = hp.smoothing(jonas, fwhm=np.radians(np.sqrt((56./60.)**2 - (20./60.)**2)))\n",
+    "jonas[np.where(jonas == hp.UNSEEN)] = np.nan\n",
+    "jonas_off = remove_off(jonas_5, has_5, 0.06, haserr)\n",
+    "\n",
+    "hp.mollview(jonas, norm='hist')\n",
+    "plt.title('Jonas')\n",
+    "plt.show()\n",
+    "\n",
+    "beta_2326 = (beta_map + cs_map * np.log(goal_f/45.))\n",
+    "ff2326 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "param_2326 = hasmap * (goal_f/408.)**beta_2326 + ff2326\n",
+    "param_2326[np.where(np.isnan(param_2326))[0]] = hp.UNSEEN\n",
+    "param_2326[np.where(param_2326 == hp.UNSEEN)[0]] = np.nan\n",
+    "param_2326[np.where(np.isnan(jonas))] = np.nan\n",
+    "\n",
+    "hp.mollview(param_2326, norm='hist')\n",
+    "plt.title('Parametric 2326 model')\n",
+    "plt.show()\n",
+    "\n",
+    "ff2326 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "sky = pysm.Sky(nside=nside, preset_strings=[pysm_sync_mod])\n",
+    "pysm2326 = sky.get_emission((goal_f/1000.) * u.GHz)[0] * 1.e-6 #gives Stokes I, Q and U in microK\n",
+    "pysm2326.value[np.where(np.isnan(jonas))[0]] = hp.UNSEEN\n",
+    "map_resol = hp.nside2resol(nside, arcmin=True) / 60.\n",
+    "pysm2326 = hp.smoothing(pysm2326.value, fwhm=np.radians(np.sqrt((56./60.)**2 - map_resol**2))) + ff2326\n",
+    "pysm2326[np.where(pysm2326 == hp.UNSEEN)] = np.nan\n",
+    "\n",
+    "hp.mollview(pysm2326, norm='hist')\n",
+    "plt.title('PYSM 2326')\n",
+    "plt.show()\n",
+    "\n",
+    "diff = jonas - jonas_off - param_2326\n",
+    "diff_pysm = jonas - jonas_off - pysm2326\n",
+    "\n",
+    "fig, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1, figsize=(6,8))\n",
+    "plt.axes(ax1)\n",
+    "hp.mollview(diff / (jonas - jonas_off),  min=-1.0, max=1.0, hold=True, title='Param 2326 MHz', cmap=cmapEC)\n",
+    "plt.axes(ax2)\n",
+    "hp.mollview(diff_pysm / (jonas - jonas_off),  min=-1.0, max=1.0, hold=True, title='PySM 2326 MHz', cmap=cmapEC)\n",
+    "plt.show()\n",
+    "\n",
+    "scat_x = (jonas - jonas_off) * gmask\n",
+    "notnans = np.where(~np.isnan(scat_x))[0]\n",
+    "\n",
+    "fig, ((ax1, ax2, ax3)) = plt.subplots(nrows=1, ncols=3, figsize=(24,6))\n",
+    "plt.axes(ax1)\n",
+    "plt.scatter(scat_x, pysm2326*gmask, s=20, edgecolors='teal', facecolors='none', label='PYSM3 Model', alpha=0.4)\n",
+    "plt.scatter(scat_x, param_2326*gmask, s=20, edgecolors='darkslateblue', facecolors='none', label='Parametric Model', alpha=0.4)\n",
+    "plt.plot(scat_x, scat_x, 'k-')\n",
+    "plt.xlim([0.0, 1+np.nanmax(param_2326*gmask)])\n",
+    "plt.ylim([0.0, 1+np.nanmax(param_2326*gmask)])\n",
+    "ax1.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax1.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.text(4.5, 8.1, f'PYSM corr: {np.corrcoef(scat_x [notnans], (pysm2326*gmask)[notnans])[0,1]:.3f} K', color='teal', fontsize=14)\n",
+    "plt.text(4.5, 7.4, f'Parametric corr: {np.corrcoef(scat_x [notnans], (param_2326*gmask)[notnans])[0,1]:.3f} K', color='darkslateblue', fontsize=14)\n",
+    "plt.xlabel('Data 2326 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.ylabel('Model 2326 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.axes(ax2)\n",
+    "plt.hist(abs(diff/(jonas - jonas_off))*100., bins=40, range=[0, 120], label='Parametric (Full)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(jonas - jonas_off))*100., bins=40, range=[0, 120], label='PYSM3 (Full)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(jonas - jonas_off))*100.), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(jonas - jonas_off))*100.), color='teal', linestyle='--')\n",
+    "ax2.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax2.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12)\n",
+    "plt.axes(ax3)\n",
+    "plt.hist(abs(diff/(jonas - jonas_off))*100., bins=40, range=[0, 120], label='Parametric (NoGal)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(jonas - jonas_off))*100.*gmask, bins=40, range=[0, 120], label='PYSM3  (NoGal)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(jonas - jonas_off))*100.*gmask), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(jonas - jonas_off))*100*gmask), color='teal', linestyle='--')\n",
+    "ax3.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax3.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12)\n",
+    "plt.show()\n",
+    "\n",
+    "print(np.nanmedian(abs(diff/(jonas - jonas_off))*100.*gmask), np.nanpercentile(abs(diff/(jonas - jonas_off))*100.*gmask, 25), np.nanpercentile(abs(diff/(jonas - jonas_off))*100.*gmask, 75))\n",
+    "print(np.nanmedian(abs(diff_pysm/(jonas - jonas_off))*100.*gmask), np.nanpercentile(abs(diff_pysm/(jonas - jonas_off))*100.*gmask, 25), np.nanpercentile(abs(diff_pysm/(jonas - jonas_off))*100.*gmask, 75))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "c5019d1d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# spass is at 2303 Mhz and 8.9 arcmin\n",
+    "goal_f = 2303.0\n",
+    "fstring = \"spass_dr1_1902_healpix_Tb.i.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    spass =  hp.read_map(filename=f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/foregrounds/spass/{fstring}\" , out =\"../input_data\")\n",
+    "    spass =  hp.read_map(f\"../input_data/{fstring}\")\n",
+    "spass = hp.ud_grade(spass, 256)\n",
+    "spass[np.where(spass < -1000)] = np.nan\n",
+    "spass[np.where(np.isnan(spass))[0]] = hp.UNSEEN\n",
+    "spass_5 = hp.smoothing(spass, fwhm=np.radians(np.sqrt((5.0)**2 - (8.9/60.)**2)))\n",
+    "spass_5[np.where(spass_5 == hp.UNSEEN)] = np.nan\n",
+    "spass= hp.smoothing(spass, fwhm=np.radians(np.sqrt((56./60.)**2 - (8.9/60.)**2)))\n",
+    "spass[np.where(spass == hp.UNSEEN)] = np.nan\n",
+    "spass_off = remove_off(spass_5, has_5, 0.05, haserr)\n",
+    "\n",
+    "hp.mollview(spass, norm='hist')\n",
+    "plt.title('SPASS')\n",
+    "plt.show()\n",
+    "\n",
+    "beta_2303 = (beta_map + cs_map * np.log(goal_f/45.))\n",
+    "ff2303 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "param_2303 = hasmap * (goal_f/408.)**beta_2303 + ff2303\n",
+    "param_2303[np.where(np.isnan(param_2303))[0]] = hp.UNSEEN\n",
+    "param_2303[np.where(param_2303 == hp.UNSEEN)[0]] = np.nan\n",
+    "param_2303[np.where(np.isnan(spass))] = np.nan\n",
+    "\n",
+    "hp.mollview(param_2303, norm='hist')\n",
+    "plt.title('Parametric 2303 model')\n",
+    "plt.show()\n",
+    "\n",
+    "sky = pysm.Sky(nside=nside, preset_strings=[pysm_sync_mod])\n",
+    "pysm2303= sky.get_emission((goal_f/1000.) * u.GHz)[0] * 1.e-6 #gives Stokes I, Q and U in microK\n",
+    "pysm2303.value[np.where(np.isnan(spass))[0]] = hp.UNSEEN\n",
+    "map_resol = hp.nside2resol(nside, arcmin=True) / 60.\n",
+    "pysm2303 = hp.smoothing(pysm2303.value, fwhm=np.radians(np.sqrt((56./60.)**2 - map_resol**2))) + ff2303\n",
+    "pysm2303[np.where(pysm2303 == hp.UNSEEN)] = np.nan\n",
+    "\n",
+    "hp.mollview(pysm2303, norm='hist')\n",
+    "plt.title('PYSM 2303')\n",
+    "plt.show()\n",
+    "\n",
+    "diff = spass - spass_off - param_2303\n",
+    "diff_pysm = spass - spass_off - pysm2303\n",
+    "\n",
+    "fig, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1, figsize=(6,8))\n",
+    "plt.axes(ax1)\n",
+    "hp.mollview(diff / (spass - spass_off),  min=-1.0, max=1.0, hold=True, title='Param 2303 MHz', cmap=cmapEC)\n",
+    "plt.axes(ax2)\n",
+    "hp.mollview(diff_pysm / (spass - spass_off),  min=-1.0, max=1.0, hold=True, title='PySM 2303 MHz', cmap=cmapEC)\n",
+    "plt.show()\n",
+    "\n",
+    "\n",
+    "scat_x = (spass - spass_off) * gmask\n",
+    "notnans = np.where(~np.isnan(scat_x))[0]\n",
+    "\n",
+    "fig, ((ax1, ax2, ax3)) = plt.subplots(nrows=1, ncols=3, figsize=(24,6))\n",
+    "plt.axes(ax1)\n",
+    "plt.scatter(scat_x, pysm2303*gmask, s=20, edgecolors='teal', facecolors='none', label='PYSM3 Model', alpha=0.4)\n",
+    "plt.scatter(scat_x , param_2303*gmask, s=20, edgecolors='darkslateblue', facecolors='none', label='Parametric Model', alpha=0.4)\n",
+    "plt.plot(scat_x, scat_x, 'k-')\n",
+    "plt.xlim([0.0, 1+np.nanmax(pysm2303*gmask)])\n",
+    "plt.ylim([0.0, 1+np.nanmax(pysm2303*gmask)])\n",
+    "ax1.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax1.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.text(0.2, 1.76, f'Parametric corr: {np.corrcoef(scat_x[notnans], (param_2303*gmask)[notnans])[0,1]:.3f} K', color='darkslateblue', fontsize=14)\n",
+    "plt.text(0.2, 1.46, f'PYSM corr: {np.corrcoef(scat_x[notnans], (pysm2303*gmask)[notnans])[0,1]:.3f} K', color='teal', fontsize=14)\n",
+    "plt.xlabel('Data 2303 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.ylabel('Model 2303 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.axes(ax2)\n",
+    "plt.hist(abs(diff/(spass - spass_off))*100., bins=40, range=[0, 120], label='Parametric (Full)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(spass - spass_off))*100., bins=40, range=[0, 120], label='PYSM3 (Full)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(spass - spass_off))*100.), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(spass - spass_off))*100.), color='teal', linestyle='--')\n",
+    "ax2.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax2.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12, loc='upper right')\n",
+    "plt.axes(ax3)\n",
+    "plt.hist(abs(diff/(spass - spass_off))*100.*gmask, bins=40, range=[0, 120], label='Parametric (NoGal)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(spass - spass_off))*100.*gmask, bins=40, range=[0, 120], label='PYSM3 (NoGal)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(spass - spass_off))*100.*gmask), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(spass - spass_off))*100.*gmask), color='teal', linestyle='--')\n",
+    "ax3.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax3.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12, loc='upper right')\n",
+    "plt.show()\n",
+    "\n",
+    "print(np.nanmedian(abs(diff/(spass - spass_off))*100.*gmask), np.nanpercentile(abs(diff/(spass - spass_off))*100.*gmask, 25), np.nanpercentile(abs(diff/(spass - spass_off))*100.*gmask, 75))\n",
+    "print(np.nanmedian(abs(diff_pysm/(spass - spass_off))*100.*gmask), np.nanpercentile(abs(diff_pysm/(spass - spass_off))*100.*gmask, 25), np.nanpercentile(abs(diff_pysm/(spass - spass_off))*100.*gmask, 75))"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "51d19e91",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# q11 is at 11000 Mhz and 1 deg\n",
+    "goal_f = 11000.0\n",
+    "fstring = \"quijote_mfi_smth_skymap_11ghz_512_dr1.fits\"\n",
+    "try :\n",
+    "    print(f\"reading {fstring} \")\n",
+    "    q11 =  hp.read_map(filename=f\"../input_data/{fstring}\")  \n",
+    "except FileNotFoundError: \n",
+    "    filename = wget.download( f\"https://lambda.gsfc.nasa.gov/data/suborbital/QUIJOTE/QUIJOTE-MFI-Release1/{fstring}\" , out =\"../input_data\")\n",
+    "    q11 =  hp.read_map(f\"../input_data/{fstring}\")\n",
+    "q11 = hp.ud_grade(q11, 256) #convert from mk to K\n",
+    "q11_5 = hp.smoothing(q11, fwhm=np.radians(np.sqrt((5.0)**2 - (60./60.)**2)))\n",
+    "q11_5[np.where(q11_5 == hp.UNSEEN)] = np.nan\n",
+    "q11[np.where(q11 == hp.UNSEEN)] = np.nan\n",
+    "q11_off = remove_off(q11_5, has_5, 0.05, haserr)\n",
+    "q11 = q11 / 1000.\n",
+    "q11_off= q11_off / 1000.\n",
+    "\n",
+    "hp.mollview(q11, norm='hist')\n",
+    "plt.title('q11')\n",
+    "plt.show()\n",
+    "\n",
+    "beta_11 = (beta_map + cs_map * np.log(goal_f/45.))\n",
+    "ff11 = Tff(Te=planck_te, nu=goal_f/1000., EM=hust)\n",
+    "param_11 = hasmap * (goal_f/408.)**beta_11 + ff11\n",
+    "param_11[np.where(np.isnan(param_11))[0]] = hp.UNSEEN\n",
+    "param_11[np.where(param_11 == hp.UNSEEN)[0]] = np.nan\n",
+    "param_11[np.where(np.isnan(q11))] = np.nan\n",
+    "\n",
+    "hp.mollview(param_11, norm='hist')\n",
+    "plt.title('Parametric 11000 model')\n",
+    "plt.show()\n",
+    "\n",
+    "sky = pysm.Sky(nside=nside, preset_strings=[pysm_sync_mod])\n",
+    "pysm11= sky.get_emission((goal_f/1000.) * u.GHz)[0] * 1.e-6 #gives Stokes I, Q and U in microK\n",
+    "pysm11.value[np.where(np.isnan(q11))[0]] = hp.UNSEEN\n",
+    "map_resol = hp.nside2resol(nside, arcmin=True) / 60.\n",
+    "pysm11 = hp.smoothing(pysm11.value, fwhm=np.radians(np.sqrt((56./60.)**2 - map_resol**2))) + ff11\n",
+    "pysm11[np.where(pysm11 == hp.UNSEEN)] = np.nan\n",
+    "\n",
+    "hp.mollview(pysm11, norm='hist')\n",
+    "plt.title('PYSM 11000')\n",
+    "plt.show()\n",
+    "\n",
+    "\n",
+    "diff = q11 - q11_off - param_11\n",
+    "diff_pysm = q11 - q11_off - pysm11\n",
+    "\n",
+    "fig, ((ax1, ax2)) = plt.subplots(nrows=2, ncols=1, figsize=(6,8))\n",
+    "plt.axes(ax1)\n",
+    "hp.mollview(diff / (q11 - q11_off),  min=-1.0, max=1.0, hold=True, title='Param 11000 MHz', cmap=cmapEC)\n",
+    "plt.axes(ax2)\n",
+    "hp.mollview(diff_pysm / (q11 - q11_off),  min=-1.0, max=1.0, hold=True, title='PySM 11000 MHz', cmap=cmapEC)\n",
+    "plt.show()\n",
+    "\n",
+    "\n",
+    "scat_x = (q11 - q11_off) * gmask\n",
+    "notnans = np.where(~np.isnan(scat_x))[0]\n",
+    "\n",
+    "fig, ((ax1, ax2, ax3)) = plt.subplots(nrows=1, ncols=3, figsize=(24,6))\n",
+    "plt.axes(ax1)\n",
+    "plt.scatter(scat_x , param_11*gmask, s=20, edgecolors='darkslateblue', facecolors='none', label='Parametric Model')\n",
+    "plt.scatter(scat_x, pysm11*gmask, s=20, edgecolors='teal', facecolors='none', label='PYSM3 Model')\n",
+    "plt.plot(scat_x, scat_x, 'k-')\n",
+    "plt.xlim([0.0, 0.03+np.nanmax(pysm11*gmask)])\n",
+    "plt.ylim([0.0, 0.03+np.nanmax(pysm11*gmask)])\n",
+    "ax1.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax1.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.text(0.02, 0.04, f'Parametric corr: {np.corrcoef(scat_x[notnans], (param_11*gmask)[notnans])[0,1]:.3f} K', color='darkslateblue', fontsize=14)\n",
+    "plt.text(0.02, 0.035, f'PYSM corr: {np.corrcoef(scat_x[notnans], (pysm11*gmask)[notnans])[0,1]:.3f} K', color='teal', fontsize=14)\n",
+    "plt.xlabel('Data 11000 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.ylabel('Model 11000 MHz Brightness Temperature (K)', fontsize=14)\n",
+    "plt.axes(ax2)\n",
+    "plt.hist(abs(diff/(q11 - q11_off))*100., bins=40, range=[0, 150], label='Parametric (Full)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(q11 - q11_off))*100., bins=40, range=[0, 150], label='PYSM3 (Full)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(q11 - q11_off))*100.), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(q11 - q11_off))*100.), color='teal', linestyle='--')\n",
+    "ax2.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax2.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12, loc='upper right')\n",
+    "plt.axes(ax3)\n",
+    "plt.hist(abs(diff/(q11 - q11_off))*100.*gmask, bins=40, range=[0, 150], label='Parametric (NoGal)', histtype='step', color='darkslateblue')\n",
+    "plt.hist(abs(diff_pysm/(q11 - q11_off))*100.*gmask, bins=40, range=[0, 150], label='PYSM3 (NoGal)', histtype='step', color='teal')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff/(q11 - q11_off))*100.*gmask), color='darkslateblue', linestyle='--')\n",
+    "plt.axvline(x=np.nanmedian(abs(diff_pysm/(q11 - q11_off))*100.*gmask), color='teal', linestyle='--')\n",
+    "ax3.tick_params(axis='both', which='major', labelsize=14)\n",
+    "ax3.tick_params(axis='both', which='minor', labelsize=14)\n",
+    "plt.xlabel('Absolute Percentage Difference (%)', fontsize=14)\n",
+    "plt.ylabel('Number of Pixels', fontsize=14)\n",
+    "plt.legend(fontsize=12, loc='upper right')\n",
+    "plt.show()\n",
+    "\n",
+    "print(np.nanmedian(abs(diff/(q11 - q11_off))*100.*gmask), np.nanpercentile(abs(diff/(q11 - q11_off))*100.*gmask, 25), np.nanpercentile(abs(diff/(q11 - q11_off))*100.*gmask, 75))\n",
+    "print(np.nanmedian(abs(diff_pysm/(q11 - q11_off))*100.*gmask), np.nanpercentile(abs(diff_pysm/(q11 - q11_off))*100.*gmask, 25), np.nanpercentile(abs(diff_pysm/(q11 - q11_off))*100.*gmask, 75))"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "env_healpy",
+   "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.10.19"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}