Skip to content

Commit

Permalink
Merge pull request #106 from galsci/synchrotron_logpoltens
Browse files Browse the repository at this point in the history
Synchrotron logpoltens
  • Loading branch information
zonca committed Mar 16, 2022
2 parents 9f51b6c + 94cb4d7 commit 29b5742
Show file tree
Hide file tree
Showing 13 changed files with 2,714 additions and 14 deletions.
6 changes: 6 additions & 0 deletions CHANGES.rst
Original file line number Diff line number Diff line change
@@ -1,3 +1,9 @@
3.4.0b2 (2022-03-15)
====================

- bugfix: fixed units error recently introduced in the CO Lines models, `reported by @mousset <https://github.com/galsci/pysm/issues/113>`_
- Synchrotron high resolution models `s4`, `s5`, `s6`, `based on Haslam and WMAP <https://github.com/galsci/pysm/pull/106>`_

3.4.0b1 (2022-03-03)
====================

Expand Down
10 changes: 8 additions & 2 deletions docs/models.rst
Original file line number Diff line number Diff line change
Expand Up @@ -50,11 +50,11 @@ Dust

- **d8**: simplified version of `d7` where the interstellar radiation field (ISRF) strength, instead of being a random realization, is fixed at 0.2. This corresponds reasonably well to a Modifield Black Body model with temperature of 20K and an index of 1.54.

- **d9**: simplified version of `d10` with a fixed spectral index of 1.48 and a fixed dust black body temperature of 19.6 K all over the sky, based on Planck 2018 results.
- **d9**: simplified version of **d10** with a fixed spectral index of 1.48 and a fixed dust black body temperature of 19.6 K all over the sky, based on Planck 2018 results.

- **d10**: single component modified black body model based on templates from the `GNILC needlet-based analysis of Planck data <https://wiki.cosmos.esa.int/planck-legacy-archive/index.php/Foreground_maps#GNILC_thermal_dust_maps>`_, with reduced contamination from CIB and point sources compared to the Commander maps used on **d1**. Small-scale fluctuations up to a $\ell_{max}$ of 16384 have been added to the templates in the `logpoltens` (Logarithm of the Polarization Fraction Tensor) formalism. Also the spectral index and dust temperature maps have random fluctuations at small scales. Available up to $N_{side}$ of 8192. Input templates are available `at NERSC <https://portal.nersc.gov/project/cmb/pysm-data/dust_gnilc/>`_. Input templates are available at 2048, 4096 and 8192; lower nside are simulated by running `hp.ud_grade` on the $N_{side}=2048$ maps. For more details about the pre-processing of the input data, see the `notebook about creating the GNILC-based template and adding small scales using the logpoltens formalism <preprocess-templates/gnilc_dust_logpoltens_templates.html>`_ and `the notebook about spectral index and dust temperature <preprocess-templates/gnilc_dust_spectralindex_Tdust.html>`_. The galactic region of about 30% of the sky within the GAL 70 Planck mask is set to the input GNILC map at 21.8' to avoid excess power in the spectra caused by the injected small scale signa, see `how the correction is applied <preprocess-templates/gnilc_dust_logpoltens_templates_galplane_prepare.html>`_ and `its impact on the power spectrum <preprocess-templates/gnilc_dust_logpoltens_templates_galplane_compare.html>`_.

- **d11**: like **d10** with stochastic small scales generated on-the-fly. It can reproduce **d10** if configured to run with a specific set of seeds and a specific $\ell_{max}$, see :py:class:`~pysm3.ModifiedBlackBodyRealization`.
- **d11**: like **d10** with stochastic small scales generated on-the-fly. It can reproduce **d10** if configured to run with a specific set of seeds and a specific $\ell_{max}$, see :py:class:`~pysm3.ModifiedBlackBodyRealization`. However reproducing **d10** is expensive because it needs to generate small scales with a $\ell_{max}=16384$ whatever output resolution is required, normally instead **d11** generates small scales just up to $\ell_{max}=3N_{side}-1$.

- **d12**: 3D model of polarized dust emission with 6 layers, based on the paper `"A 3-D model of polarised dust emission in the Milky Way" <https://arxiv.org/abs/1706.04162>`_, named MKD based on the names of the authors. Each layer has different templates, spectral index and dust temperature. All maps were generated at N_side 2048 with the Planck Sky Model (PSM) by Jacques Delabrouille.

Expand All @@ -67,6 +67,12 @@ Synchrotron

- **s3**: a power law with a curved index. The model uses the same index map as the nominal model, plus a curvature term. We use the best-fit curvature amplitude of -0.052 found in Kogut, A. 2012, ApJ, 753, 110, pivoted at 23 GHz.

- **s4**: simplified version of **s5** with a fixed spectral index of -3.1 all over the sky.

- **s5**: power law model based on the same templates of **s1**, Haslam in temperature and WMAP 9 year 23 GHz in polarization. Small-scale fluctuations up to a $\ell_{max}$ of 16384 have been added to the templates in the `logpoltens` (Logarithm of the Polarization Fraction Tensor) formalism. The spectral index map from **s1** has been rescaled `based on S-PASS <https://arxiv.org/abs/1802.01145>`_ and had small scales added to upgrade it up to $N_{side}$ 8192. Input templates are available `at NERSC <https://portal.nersc.gov/project/cmb/pysm-data/synch/>`_ at 2048, 4096 and 8192; lower nside are simulated by running `hp.ud_grade` on the $N_{side}=2048$ maps. For more details about the pre-processing of the input data, see the `notebook about creating the templates and adding small scales using the logpoltens formalism <preprocess-templates/synchrotron_template_logpoltens.html>`_ and `the notebook about spectral index <preprocess-templates/synchrotron_beta.html>`_.

- **s6**: like **s5** with stochastic small scales generated on-the-fly. It can reproduce **s5** if configured to run with a specific set of seeds and a specific $\ell_{max}$, see :py:class:`~pysm3.PowerLawRealization`. However reproducing **s5** is expensive because it needs to generate small scales with a $\ell_{max}=16384$ whatever output resolution is required, normally instead **s6** generates small scales just up to $\ell_{max}=3N_{side}-1$.


AME
===
Expand Down
385 changes: 385 additions & 0 deletions docs/preprocess-templates/synchrotron_beta.ipynb

Large diffs are not rendered by default.

1,552 changes: 1,552 additions & 0 deletions docs/preprocess-templates/synchrotron_template_logpoltens.ipynb

Large diffs are not rendered by default.

258 changes: 258 additions & 0 deletions docs/preprocess-templates/utils_synch_generate_map.ipynb
Original file line number Diff line number Diff line change
@@ -0,0 +1,258 @@
{
"cells": [
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import healpy as hp\n",
"from pathlib import Path"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"datadir=Path(\"data/\")\n",
"output_dir = Path(\"production-data/synch\")\n",
"output_dir_raw = output_dir / \"raw\""
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"parameters"
]
},
"outputs": [],
"source": [
"output_nside = 8192"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"output_lmax = min(3*output_nside - 1, 8192*2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"comp = \"synch\"\n",
"sub = \"template\""
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Large scales"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"largescale_filename = list(output_dir_raw.glob(f\"{comp}*largescale*{sub}*.fits\"))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"assert len(largescale_filename) == 1"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"largescale_filename = largescale_filename[0]"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"largescale_filename"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"alm_log_pol_tens_large_scale = hp.read_alm(largescale_filename, hdu=(1,2,3))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"map_log_pol_tens_large_scale = hp.alm2map(alm_log_pol_tens_large_scale.astype(np.complex128), nside=output_nside)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Small scales modulation"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"modulate_alm = { k:hp.read_alm(output_dir_raw/f\"synch_{k}_modulation_alms_lmax1535.fits\").astype(np.complex128) for k in [\"temperature\",\"polarization\"] }"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Small scales"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"cl_small_scale = hp.read_cl(output_dir_raw / \"synch_small_scales_logpoltens_cl_lmax16384.fits\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"synalm_lmax = 8192*2 # for reproducibility\n",
"# synalm_lmax = output_lmax\n",
"np.random.seed(555)\n",
"\n",
"alm_log_pol_tens_small_scale = hp.synalm(\n",
" list(cl_small_scale) + [np.zeros_like(cl_small_scale[0])] * 3,\n",
" lmax=synalm_lmax,\n",
" new=True,\n",
")\n",
"\n",
"alm_log_pol_tens_small_scale = [hp.almxfl(each, np.ones(3*output_nside-1)) for each in alm_log_pol_tens_small_scale]\n",
"map_log_pol_tens_small_scale = hp.alm2map(alm_log_pol_tens_small_scale, nside=output_nside)\n",
"map_log_pol_tens_small_scale[0] *= hp.alm2map(modulate_alm[\"temperature\"], output_nside)\n",
"map_log_pol_tens_small_scale[1:] *= hp.alm2map(modulate_alm[\"polarization\"], output_nside)\n",
"assert np.isnan(map_log_pol_tens_small_scale).sum() == 0"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Combine scales\n",
"\n",
"* Combine small and large scale maps\n",
"* Transform from logpoltens to IQU\n",
"* Write output map"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"map_log_pol_tens = map_log_pol_tens_large_scale + map_log_pol_tens_small_scale"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from pysm3.utils import log_pol_tens_to_map, add_metadata"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"output_map = log_pol_tens_to_map(map_log_pol_tens)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"hp.write_map(output_dir / f\"synch_template_nside{output_nside}.fits\", output_map, dtype=np.float32, overwrite=True)"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"add_metadata([output_dir / f\"synch_template_nside{output_nside}.fits\"], coord=\"G\", unit=\"uK_RJ\", ref_freq=\"23 GHz\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "condanamaster2",
"language": "python",
"name": "condanamaster2"
},
"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.8.0"
}
},
"nbformat": 4,
"nbformat_minor": 4
}

0 comments on commit 29b5742

Please sign in to comment.