{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Haslam index" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from healpy.newvisufunc import projview\n", "from healpy.rotator import Rotator\n", "import healpy as hp\n", "import traceback\n", "\n", "from pylfmap import LFmap\n", "from SSM import SSM\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import matplotlib.pyplot as plt\n", "import itertools\n", "\n", "import os\n", "import sys\n", "from radiocalibrationtoolkit import *" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# some global plot settings\n", "plt.rcParams[\"axes.labelweight\"] = \"bold\"\n", "plt.rcParams[\"font.weight\"] = \"bold\"\n", "plt.rcParams['font.size'] = 16\n", "plt.rcParams['legend.fontsize']= 14\n", "\n", "plt.rcParams['xtick.major.width']= 2\n", "plt.rcParams['ytick.major.width']= 2\n", "\n", "plt.rcParams['xtick.major.size']= 5\n", "plt.rcParams['ytick.major.size']= 5\n", "\n", "plt.rcParams['xtick.labelsize']= 14\n", "plt.rcParams['ytick.labelsize']= 14" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Prepare objects\n", "lfmap = LFmap()\n", "lfss = LowFrequencySkyModel(freq_unit=\"MHz\")\n", "gsm2008 = GlobalSkyModel(freq_unit=\"MHz\")\n", "gsm2016 = GlobalSkyModel2016(freq_unit=\"MHz\")\n", "ssm = SSM()\n", "gmoss = GMOSS()\n", "ulsa_fdi = ULSA(index_type='freq_dependent_index')\n", "ulsa_ci = ULSA(index_type='constant_index')\n", "ulsa_dpi = ULSA(index_type='direction_dependent_index')" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "spectral_incides = np.linspace(-2.8,-2.2,31)\n", "frequency_range = np.arange(30, 85, 5)\n", "\n", "# frequency_MHz = 45\n", "\n", "diff_dict = {}\n", "for frequency_MHz in frequency_range:\n", "\n", " lfmap_map = hp.rotator.Rotator.rotate_map_pixel(\n", " Rotator(coord=[\"C\", \"G\"]), lfmap.generate(frequency_MHz)\n", " )\n", " lfmap_map = lfmap.generate(frequency_MHz)\n", " lfss_map = lfss.generate(frequency_MHz)\n", " gsm2008_map = gsm2008.generate(frequency_MHz)\n", " gsm2016_map = gsm2016.generate(frequency_MHz)\n", " ssm_map = ssm.generate(frequency_MHz)\n", " gmoss_map = gmoss.generate(frequency_MHz)\n", " ulsa_fdi_map = ulsa_fdi.generate(frequency_MHz)\n", " ulsa_ci_map = ulsa_ci.generate(frequency_MHz)\n", " ulsa_dpi_map = ulsa_dpi.generate(frequency_MHz)\n", "\n", " # convert to same NSIDE\n", " new_nside = 64\n", " lfmap_map_N = hp.ma(hp.pixelfunc.ud_grade(lfmap_map, new_nside))\n", " lfss_map_N = hp.ma(hp.pixelfunc.ud_grade(lfss_map, new_nside))\n", " gsm2008_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2008_map, new_nside))\n", " gsm2016_map_N = hp.ma(hp.pixelfunc.ud_grade(gsm2016_map, new_nside))\n", " ssm_map_N = hp.ma(hp.pixelfunc.ud_grade(ssm_map, new_nside))\n", " gmoss_map_N = hp.ma(hp.pixelfunc.ud_grade(gmoss_map, new_nside))\n", " ulsa_fdi_map_N = hp.ma(hp.pixelfunc.ud_grade(ulsa_fdi_map, new_nside))\n", "\n", " map_dict = {\n", " \"LFSS\": lfss_map_N,\n", " \"GSM08\": gsm2008_map_N,\n", " \"GSM16\": gsm2016_map_N,\n", " \"LFmap\": lfmap_map_N,\n", " \"SSM\": ssm_map_N,\n", " \"GMOSS\": gmoss_map_N,\n", " \"ULSA\": ulsa_fdi_map_N,\n", " }\n", "\n", " mean_diffs_list = []\n", " for spectral_index in spectral_incides:\n", " haslam = HaslamSkyModel(freq_unit=\"MHz\", spectral_index=spectral_index)\n", " haslam_map = haslam.generate(frequency_MHz)\n", " haslam_map_N = hp.ma(hp.pixelfunc.ud_grade(haslam_map, new_nside))\n", "\n", " diff = []\n", " for key in map_dict:\n", " diff.append(np.sum(haslam_map_N)/np.sum(map_dict[key]))\n", "\n", " mean_diffs_list.append(np.mean(diff))\n", "\n", " diff_dict[frequency_MHz] = mean_diffs_list" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, ax = plt.subplots()\n", "best = []\n", "for key in diff_dict:\n", " ax.plot(spectral_incides, np.asarray(diff_dict[key]), label='{} MHz'.format(key))\n", " best.append(spectral_incides[np.argmin( np.abs(np.asarray(diff_dict[key])-1 ))])\n", "\n", "ax.axes.axhline(1, linestyle='--', color='g')\n", "# ax.axes.axvspan(np.min(best), np.max(best), alpha=0.5, color='g')\n", "ax.axes.axvline(np.mean(best), alpha=0.5, color='r', linestyle='--')\n", "\n", "ax.set_xlabel('spectral index')\n", "ax.set_ylabel('average ratio')\n", "ax.legend()\n", "\n", "plt.draw() # this is required, or the ticklabels may not exist (yet) at the next step\n", "labels = [w.get_text() for w in ax.get_xticklabels()]\n", "locs=list(ax.get_xticks())\n", "labels+=['\\nbest:{:.2f}'.format(np.mean(best))]\n", "locs+=[np.mean(best)]\n", "ax.set_xticklabels(labels)\n", "ax.set_xticks(locs)\n", "plt.gca().get_xticklabels()[-1].set_color('red')\n", "\n", "ax.xaxis.get_major_ticks()[-1].tick1line.set_markeredgewidth(3)\n", "ax.xaxis.get_major_ticks()[-1].tick1line.set_markersize(15)\n", "ax.xaxis.get_major_ticks()[-1].tick1line.set_markeredgecolor(\"red\")\n", "\n", "plt.draw()\n", "\n", "print('Best: {}'.format(np.mean(best)))" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "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.6" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": true, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false }, "vscode": { "interpreter": { "hash": "e7370f93d1d0cde622a1f8e1c04877d8463912d04d973331ad4851f04de6915a" } } }, "nbformat": 4, "nbformat_minor": 2 }