{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Codes for making plots for O3 anisotropic stochastic paper" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Importing packages and setting things up; \n", "## Before proceeding make sure that all the packages listed below are already installed " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import healpy as hp\n", "import h5py\n", "import scipy.io as sio\n", "import matplotlib # requires matplotlib version >= 3.2.2\n", "import matplotlib.pyplot as plt\n", "import matplotlib.ticker as ticker\n", "from scipy.io import loadmat\n", "import matplotlib.patches as mpatches\n", "from matplotlib.lines import Line2D\n", "import warnings\n", "warnings.filterwarnings(\"ignore\", category=UserWarning) " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Parameters for plotting" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fontsize = 24\n", "matplotlib.rcParams.update({\n", " \"font.size\": fontsize,\n", " \"axes.titlesize\": fontsize,\n", " \"axes.labelsize\": fontsize,\n", " \"xtick.labelsize\": fontsize,\n", " \"ytick.labelsize\": fontsize,\n", " \"xtick.major.size\": fontsize * 0.8,\n", " \"ytick.major.size\": fontsize * 0.8,\n", " \"legend.fontsize\": fontsize * 0.9,\n", " \"font.family\": \"Times new Roman\",\n", " \"figure.dpi\": 100,\n", " \"savefig.dpi\": 300,\n", " \"text.usetex\": True,\n", " \"path.simplify\": True,\n", " \"figure.figsize\": (8, 6)\n", "})\n", "\n", "# Class to deal in an easy way with number of digits in colorbars (used for ULs and sigma maps)\n", "# It allows to choose the number of digits for tickslabels and keeping 10^N at the right of the colorbar\n", "# instead of having 1eN for each tickslabel.\n", "\n", "class FormatScalarFormatter(matplotlib.ticker.ScalarFormatter):\n", " def __init__(self, fformat=\"%1.1f\", offset=True, mathText=True):\n", " self.fformat = fformat\n", " matplotlib.ticker.ScalarFormatter.__init__(self,useOffset=offset,\n", " useMathText=mathText)\n", " def _set_format(self):\n", " self.format = self.fformat\n", " if self._useMathText:\n", " #self.format = '$%s$' % matplotlib.ticker._mathdefault(self.format)\n", " self.format = '$%s$' % ('\\\\mathdefault{%s}' % self.format)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Variable definitions" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# variable corresponding to spectral shapes of signal models \n", "alphas = [\"0\", \"2over3\", \"3\"]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 2: Top row (BBR SNR maps)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# SNR sky-maps for BBR analysis\n", "\n", "# Create figure\n", "fig = plt.figure(figsize=(24, 8))\n", "\n", "# Limits for the colorbars\n", "vlims = [2.55, 3.72, 3.63]\n", "\n", "# Loop over the *combined_maps.hdf5 files for the three spectral indices 0, 2/3 and 3\n", "for alpha, vlim, index in zip (alphas, vlims, range(0,len(alphas))):\n", " \n", " # SNR *.hdf5 file\n", " snr_file = \"./datafiles/bbr_a\" + alpha + \"_o1_o2_o3_combined_maps.hdf5\"\n", " \n", " # Read the file and create a dictionary-like variable; the file has three variables SNR (map_snr_total), signal estimate (ptEst_total) and its sigma (sig_total)\n", " combined_snr = h5py.File(snr_file, \"r\")\n", "\n", " # Select the map of interest (SNR map)\n", " snr_map = np.array(combined_snr[\"map_snr_total\"])\n", " print(\"alpha:%s, max snr:%.1f\" % (alpha, np.max(snr_map)))\n", "\n", " \n", " # Close the file\n", " combined_snr.close()\n", " \n", " # Plot the sky-map for a given alpha as subplot using healpy Mollweide projection\n", " hp.mollview(np.real(np.squeeze(snr_map)),rot=(180,0,0),flip='astro', title=\"\",nest=False, format = \"%.2f\", cbar=False, sub=(1, 3, index+1), margins=[0.0125,0.0125, 0.0125, 0.0125], min=-vlim, max=vlim)\n", " \n", " # Add graticule\n", " hp.graticule(45,90)\n", " \n", " # Add useful annotations on the map for right ascension and declination\n", " plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate(r'$45^{\\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)\n", " plt.annotate(r'$0^{\\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)\n", " plt.annotate(r'$-45^{\\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) \n", " \n", " # Trick to make colorbars all the same (healpy colorbar does not leave much freedom of customization and it is not used)\n", " fig = plt.gcf()\n", " ax = plt.gca() \n", " image = ax.get_images()[0]\n", " cmap = fig.colorbar(image, ax=ax, shrink=0.7, orientation=\"horizontal\", pad=0.15, aspect = 40, ticks=[-2,0,2])\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 2: Bottom row (BBR Upper limit maps)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Upper limit sky-maps for BBR analysis\n", "\n", "fig = plt.figure(figsize=(24, 8))\n", "\n", "# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3\n", "for alpha, index in zip (alphas, range(0,len(alphas))):\n", "\n", " # upper limit *.mat file\n", " ul_dir = \"./datafiles/bbr_a\"+ alpha +\"_zl_limits.mat\"\n", " \n", " # Read the file and create a dictionary-like variable; the file has variables correponding to upplimits from O1+O2, just O3 (for three baselines) and O1+O2+O3\n", " combined_ul = sio.loadmat(ul_dir)\n", " ul_map = np.squeeze(np.array(combined_ul[\"flux_ul_combined\"]))\n", " print(\"alpha:%s, min UL:%.1e, max UL:%.1e\" % (alpha, np.min(ul_map), np.max(ul_map)))\n", " \n", " # Plot the sky-map for a given alpha as subplot using healpy Mollweide projection\n", " hp.mollview(np.real(np.squeeze(ul_map)),rot=(180,0,0),flip='astro', title=\"\",nest=False, format = \"%.2g\", cbar=False, sub=(1, 3, index+1), margins=[0.0125,0.0125, 0.0125, 0.0125])\n", "\n", " # Add graticule\n", " hp.graticule(45,90)\n", " \n", " # Add useful annotations on the map for right ascension and declination\n", " plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate(r'$45^{\\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)\n", " plt.annotate(r'$0^{\\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)\n", " plt.annotate(r'$-45^{\\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) \n", " \n", " \n", " fig = plt.gcf()\n", " ax = plt.gca() \n", " image = ax.get_images()[0]\n", " \n", " # Further colorbar customization\n", " nticks = 4\n", " interval = 1./(nticks-1)\n", " diff= np.max(ul_map) -np.min(ul_map)\n", " ticks_no_round = [np.min(ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]\n", " \n", " cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter(\"%.1f\"), shrink=0.7, orientation=\"horizontal\", pad=0.15, aspect = 30)\n", " \n", " # Adding units as colorbar label\n", " cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20)\n", " cmap.ax.tick_params(labelsize=fontsize)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 3: Top row (SHD SNR maps)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Limits for the colorbars\n", "vlims = [2.18, 3.72, 3.57]\n", "\n", "fig = plt.figure(figsize=(24,8))\n", "# Loop over the shd *.mat files for the three spectral indices 0, 2/3 and 3\n", "for i, (alpha, vlim) in enumerate(zip(alphas, vlims)):\n", " ax = fig.add_subplot(130+(i+1), projection=\"mollweide\")\n", " pproc_file = \"./datafiles/shd_ul_a\"+ alpha + \"_o1_o2_o3_shd_map_files.mat\"\n", " shd_data = loadmat(pproc_file)\n", " \n", " # Get the signal estimate (map0) and its sigma (sigmaPix0) for different pixels on the sky\n", " p0_map = shd_data[\"map0\"]\n", " sigma_map = shd_data[\"sigmaPix0\"]\n", " sig_map = sigma_map.reshape(360,181).T\n", " \n", " # Calculate SNR maps\n", " snr_map = (p0_map / sigma_map).reshape(360, 181).T\n", " print(\"alpha:%s, max snr:%.1f\" % (alpha, np.max(snr_map)))\n", " \n", " # Make the X, Y mesh and plot\n", " x = np.linspace(-np.pi,np.pi,360)\n", " y = np.linspace(-np.pi/2,np.pi/2,181)\n", " X,Y = np.meshgrid(x,y)\n", " image = ax.pcolormesh(X,Y,snr_map[:, ::-1], vmin=-vlim, vmax=vlim)\n", " \n", " # formate x and y tickes \n", " raval = [-np.pi/2,0,np.pi/2]\n", " ralab = ['18h','12h','6h']\n", " ax.set_xticks(raval)\n", " xtick_labels = ax.set_xticklabels(ralab, color=\"w\",fontweight='bold',fontsize=20)\n", " declval = [-np.pi/4,0,np.pi/4]\n", " decllab = [r'$-45^{\\circ}$',r'$0^{\\circ}$',r'$45^{\\circ}$']\n", " ax.set_yticks(declval)\n", " ax.set_yticklabels(decllab,fontsize=20)\n", " plt.grid(True, color=\"k\", linestyle=\":\")\n", " x = np.array([-np.radians(180), np.radians(180)])\n", " y = np.array([0, 0])\n", " plt.plot(x, y, ls='-', color=\"k\", lw=0.8)\n", " cmp = plt.colorbar(image, shrink=0.7, ax=ax, orientation=\"horizontal\", aspect=40)\n", " cmap.set_label('', fontsize = 10)\n", " for t in cmap.ax.get_xticklabels(): # reduce the font size of xlabel\n", " t.set_fontsize(10)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 3: bottom row (SHD UL maps)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig = plt.figure(figsize=(24,8))\n", "\n", "# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3\n", "for i, alpha in enumerate(alphas):\n", " ax = fig.add_subplot(130+(i+1), projection=\"mollweide\")\n", "\n", " # upper limit *.mat file\n", " pproc_file = \"./datafiles/shd_a\" + alpha + \"_limits_o1_o2_o3_HLV_ZL.mat\"\n", " shd_data = loadmat(pproc_file)\n", " omega_ul_map = shd_data[\"omega_ul\"].reshape(360, 181).T\n", " print(\"alpha:%s, min UL:%.1e, max UL:%.1e\" % (alpha, np.min(omega_ul_map), np.max(omega_ul_map)))\n", " \n", " # Make the X, Y mesh and plot\n", " x = np.linspace(-np.pi,np.pi,360)\n", " y = np.linspace(-np.pi/2,np.pi/2,181)\n", " X,Y = np.meshgrid(x,y)\n", " image = ax.pcolormesh(X,Y,omega_ul_map[:, ::-1])\n", "\n", " # formate x and y tickes \n", " raval = [-np.pi/2,0,np.pi/2]\n", " ralab = ['18h','12h','6h']\n", " ax.set_xticks(raval)\n", " xtick_labels = ax.set_xticklabels(ralab, color=\"w\",fontweight='bold',fontsize=20)\n", " declval = [-np.pi/4,0,np.pi/4]\n", " decllab = ['$-45^{\\circ}$','$0^{\\circ}$','$45^{\\circ}$']\n", " ax.set_yticks(declval)\n", " ax.set_yticklabels(decllab,fontsize=20)\n", " plt.grid(True, color=\"k\", linestyle=\":\")\n", " x = np.array([-np.radians(180), np.radians(180)])\n", " y = np.array([0, 0])\n", " plt.plot(x, y, ls='-', color=\"k\", lw=0.8)\n", " \n", " # Further colorbar customization\n", " nticks = 4\n", " interval = 1./(nticks-1)\n", " diff= np.max(omega_ul_map) -np.min(omega_ul_map)\n", " ticks_no_round = [np.min(omega_ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]\n", " \n", " cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter(\"%.1f\", mathText=False), shrink=0.7, orientation=\"horizontal\", pad=0.15, aspect = 40)\n", " # Adding units as colorbar label\n", " cmap.set_label(r'$[\\mathrm{sr}^{-1}]$', fontsize = fontsize,labelpad=20.)\n", " for t in cmap.ax.get_xticklabels():\n", " t.set_fontsize(20)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 4: C_l upper limits" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "H0 = 2.2005e-18 # Hubble constant \n", "fref = 25.0 # reference frequency\n", "const = 2.0 * np.pi**2 * fref**3 / (3 * H0**2) # normalization constant\n", "\n", "colors = ['xkcd:black', 'xkcd:sky blue', 'xkcd:orange']\n", "fig = plt.figure()\n", "ax = fig.add_subplot(111)\n", "\n", "# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3\n", "for alpha, color in zip(alphas, colors):\n", " pproc_file = \"./datafiles/shd_ul_a\" + alpha + \"_o1_o2_o3_shd_map_files.mat\"\n", " shd_data = loadmat(pproc_file)\n", " lmax = shd_data['L']\n", " Cell_ul = shd_data[\"omega_cl_ul\"]\n", " Cell = shd_data['omega_cl']\n", " ax.semilogy(np.arange(lmax + 1), np.sqrt(Cell_ul)[0], color=color, linestyle=\"None\", marker=\"v\", markersize = 12, label=r\"$\\alpha=%s$\" % alpha.replace(\"over\", \"/\"))\n", "\n", "ax.grid(which=\"major\", linestyle=\"-\")\n", "ax.grid(which=\"minor\", linestyle=\":\")\n", "ax.legend(loc=\"upper right\")\n", "ax.set_xticks(np.arange(0, 17, 2))\n", "ax.set_yticks([1e-8, 1e-9, 1e-10])\n", "ax.set_ylim((0.8e-10, 0.2e-7))\n", "ax.set_xlabel(\"$\\ell$\")\n", "ax.set_ylabel(\"$C_\\ell^{1/2}[\\mathrm{sr}^{-1}]$\")\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 5: Fisher Matrix comparison" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "baselines = [\"HL\", \"HLV\"]\n", "colors = ['xkcd:black', 'xkcd:sky blue', 'xkcd:orange']\n", "fkeep = 2/3 # keep 2/3 of dominant eigen values\n", "yvals = []\n", "fig = plt.figure(figsize=(6,6))\n", "ax = fig.add_subplot(111)\n", "\n", "# Loop over the *.mat files for the three spectral indices 0, 2/3 and 3; for HL and HLV network\n", "for i, (alpha, color) in enumerate(zip(alphas, colors)): \n", " for ii, (baseline, style) in enumerate(zip(baselines, [\"dotted\", \"solid\"])):\n", " matfile = \"./datafiles/shd_ul_a\" + alpha + \"_shd_map_files_\" + baseline + \".mat\"\n", " Fisher = loadmat(matfile)[\"fisher0\"]\n", " eigs, _ = np.linalg.eigh(Fisher) # calculate eigen values\n", " eigs = eigs[::-1]\n", " yvals.append(eigs/eigs[0])\n", " if alpha == '2over3':\n", " lab = r'$\\alpha=2/3$' + ', ' + baseline\n", " else: \n", " lab = r'$\\alpha=$%s'%alpha + ', ' + baseline\n", " ax.loglog(np.arange(len(eigs)), eigs/eigs[0], color=color, linestyle=style, label=lab)\n", " ax.axvline(len(eigs) * fkeep, linestyle=\"--\", color=color)\n", "_, _ = ax.get_legend_handles_labels()\n", "legend_elements1 = [mpatches.Patch(color=color, label=r\"$\\alpha=%s$\" % alpha.replace(\"over\", \"/\")) for color, alpha in zip(colors, alphas)]\n", "legend_elements2 = [Line2D([0], [0], color=\"k\", lw=2, label=label, linestyle=ls) for ls, label in zip([\"-\", \":\"], [\"HLV\", \"HL\"])]\n", "leg1 = ax.legend(handles=legend_elements1, loc='lower left', frameon=True, fontsize=18)\n", "leg2 = ax.legend(handles=legend_elements2, loc=(0.03, 0.3), frameon=True, fontsize=18)\n", "ax.add_artist(leg1)\n", "for line in leg2.get_lines():\n", " line.set_linewidth(2)\n", "ax.set_xlabel(\"Mode index\",fontsize=20)\n", "ax.set_ylabel(\"Eigenvalues\",fontsize=20)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 6: NBR plots" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# function to NBR plot upper limits for the three directions of interest (Sco X-1, SN 1987a, Galactic Center )\n", "def plotUL(point):\n", " # point = {scox1,sn1987a,gc}\n", " nbr_file = \"./datafiles/NBR_\" + point + \"_ul.hdf5\"\n", " struct = h5py.File(nbr_file, \"r\") # load the corresponding file\n", " freqs = np.array(struct[\"f_all\"])\n", " ul = np.array(struct[\"ul_\" + point])\n", " onesigma = np.array(struct[\"onesigma_\" + point])\n", " \n", " plt.xlim([np.amin(freqs),np.amax(freqs)]);\n", " plt.ylim([5e-26,1e-22]);\n", " plt.grid(which = \"both\", lw = 0.2, ls = \"--\")\n", " plt.loglog(freqs, ul, c = \"dimgrey\", lw = 0.3, label = \"HLV: O1+O2+O3: 95 % UL\" )\n", " plt.loglog(freqs, onesigma, c = \"k\", lw = 0.5, label = \"$1\\sigma$ sensitivity\")\n", "\n", " plt.ylabel(\"Strain amplitude ($h_0$)\")\n", " plt.xlabel(\"Frequency [Hz]\")\n", " plt.legend()\n", " if point == \"scox1\":\n", " plt.title(\"Sco X-1\")\n", " plt.savefig('ul_and_onesig-scox1-HLV.png',bbox_inches='tight')\n", " if point == \"sn1987a\":\n", " plt.title(\"SN 1987a\")\n", " plt.savefig('ul_and_onesig-sn-HLV.png',bbox_inches='tight')\n", " if point == \"gc\":\n", " plt.title(\"Galactic Center\")\n", " plt.savefig('ul_and_onesig-gc-HLV.png',bbox_inches='tight')\n", " return \"Plotting \" + point" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotUL(\"scox1\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotUL(\"sn1987a\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotUL(\"gc\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 7: Sigma maps for single baselines (BBR)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# HL, HV, LV 1-sigma maps\n", "\n", "# Normalization constant such that the 1-sigma quantity corresponds to GW flux\n", "G = 6.6726e-8\n", "c = 2.9979e10\n", "fref = 25.0\n", "flux_factor = (c**3) *np.pi*(fref**2)/(4*G)\n", "\n", "# Same as for ULs maps, but with 9 subplots of 1-sigma sensitivity for three baselines and three spectral indices\n", "fig = plt.figure(figsize=(24,18))\n", "baselines = [\"HL\", \"HV\", \"LV\"]\n", "for ii, baseline in enumerate(baselines):\n", " print(\"Baseline:%s \" % (baseline))\n", " for jj, alpha in enumerate (alphas): \n", " map_dir = \"./datafiles/bbr_a{0}_{1}_Map_dirty32.hdf5\".format(alpha, baseline)\n", " maps = h5py.File(map_dir, \"r\") # load file\n", " print(\"alpha:%s, max snr:%.1f\" % (alpha, np.max(np.array(maps['map_snr_total'])))) \n", " map_sigma = flux_factor*np.array(maps['sig_total']) # express 1-sigma in GW flux\n", " hp.mollview(np.real(np.squeeze(map_sigma)),rot=(180,0,0),flip='astro', title=\"\",nest=False, format = \"%.2f\", cbar=False, sub=(3, 3, 3*ii+ jj+1), margins=[0.0125,0.0125, 0.0125, 0.0125])\n", " hp.graticule(45,90)\n", " plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)\n", " plt.annotate(r'$45^{\\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)\n", " plt.annotate(r'$0^{\\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)\n", " plt.annotate(r'$-45^{\\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) \n", " if alpha=='0':\n", " if baseline=='HL':\n", " plt.annotate('$\\\\bf{HL}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)\n", " if baseline=='HV':\n", " plt.annotate('$\\\\bf{HV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)\n", " if baseline=='LV':\n", " plt.annotate('$\\\\bf{LV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)\n", " \n", " fig = plt.gcf()\n", " ax = plt.gca() \n", " image = ax.get_images()[0]\n", " nticks = 4\n", " interval = 1./(nticks-1)\n", " diff= np.max(map_sigma) -np.min(map_sigma)\n", " ticks_no_round = [np.min(map_sigma)+(diff)*interval*x for x in np.arange(0,nticks)]\n", " fmt = FormatScalarFormatter(\"%.1f\")\n", " cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=fmt, shrink=0.7, orientation=\"horizontal\", pad=0.15, aspect = 40)\n", " cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20) \n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 8: Sigma maps for single baselines (SHD)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "baselines = {'HL','HV','LV'}\n", "fig = plt.figure(figsize=(24,18))\n", "\n", "# Same as for ULs maps, but with 9 subplots of 1-sigma sensitivity for three baselines and three spectral indices\n", "for ii, baseline in enumerate(baselines):\n", " print(\"Baseline:%s \" % (baseline))\n", " for i, alpha in enumerate(alphas):\n", " ax = fig.add_subplot(330+(i+1)+3*ii, projection=\"mollweide\")\n", " pproc_file = \"./datafiles/shd_ul_a\" + alpha + \"_shd_map_files_\" + baseline + \".mat\"\n", " shd_data = loadmat(pproc_file) # load data file\n", " p0_map = shd_data[\"map0\"]\n", " sigma_map = shd_data[\"sigmaPix0\"]\n", "\n", " sig_map = sigma_map.reshape(360,181).T * const # express 1-sigma as a normalized quantity\n", " snr_map = (p0_map / sigma_map).reshape(360, 181).T\n", " print(\"alpha:%s, max snr:%.1f\" % (alpha, np.max(snr_map)))\n", " x = np.linspace(-np.pi,np.pi,360)\n", " y = np.linspace(-np.pi/2,np.pi/2,181)\n", " X,Y = np.meshgrid(x,y)\n", "\n", " # make plot\n", " image = ax.pcolormesh(X,Y,sig_map[:, ::-1], vmin=np.min(sig_map), vmax=np.max(sig_map))\n", " raval = [-np.pi/2,0,np.pi/2]\n", " ralab = ['18h','12h','6h']\n", " ax.set_xticks(raval)\n", " xtick_labels = ax.set_xticklabels(ralab, color=\"w\",fontsize=20)\n", " declval = [-np.pi/4,0,np.pi/4]\n", " decllab = ['$-45^{\\circ}$','$0^{\\circ}$','$45^{\\circ}$']\n", " ax.set_yticks(declval)\n", " ax.set_yticklabels(decllab,fontsize=20)\n", " plt.grid(True, color=\"k\", linestyle=\":\")\n", " x = np.array([-np.radians(180), np.radians(180)])\n", " y = np.array([0, 0])\n", " plt.plot(x, y, ls='-', color=\"k\", lw=0.8)\n", " \n", " nticks = 4\n", " interval = 1./(nticks-1)\n", " diff= np.max(sig_map) -np.min(sig_map)\n", " ticks_no_round = [np.min(sig_map)+(diff)*interval*x for x in np.arange(0,nticks)]\n", "\n", " cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter(\"%.1f\", mathText=False), shrink=0.7, orientation=\"horizontal\", pad=0.15, aspect = 40)\n", " cmap.set_label(r'$[\\mathrm{sr}^{-1}]$', fontsize = fontsize, labelpad=20.)\n", " for t in cmap.ax.get_xticklabels():\n", " t.set_fontsize(20)\n", " if alpha=='0':\n", " if baseline=='HL':\n", " plt.annotate('$\\\\bf{HL}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)\n", " if baseline=='HV':\n", " plt.annotate('$\\\\bf{HV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)\n", " if baseline=='LV':\n", " plt.annotate('$\\\\bf{LV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Figure 9: Sigma maps for single baselines (NBR)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# function to NBR 1-sigma plots for the three directions of interest (Sco X-1, SN 1987a, Galactic Center) and three baselines\n", "def plotSigma(point):\n", " # point = {scox1,sn1987a,gc}\n", " maps_mat_output = h5py.File(\"./datafiles/NBR_\" + point + '_sigma_HLV.hdf5', \"r\")\n", " tempstruct = h5py.File(\"./datafiles/NBR_scox1_ul.hdf5\", \"r\")\n", " f_all = np.array(tempstruct[\"f_all\"])\n", " plt.loglog(f_all,np.array(maps_mat_output[\"sigma_HL\"]), label='HL', lw = 2,color='black',ls = '-')\n", " plt.loglog(f_all, np.array(maps_mat_output[\"sigma_HV\"]), label='HV',lw = 2,color='skyblue',ls = '-.')\n", " plt.loglog(f_all, np.array(maps_mat_output[\"sigma_LV\"]), label='LV',lw = 2,color='orange',ls=':')\n", " plt.legend(bbox_to_anchor=(1,1), loc=\"upper right\")\n", " plt.grid(which = \"both\", lw = 0.2, ls = \"--\")\n", " plt.xlabel('Frequency [Hz]')\n", " plt.ylabel('Sigma(f)')\n", " if point == \"scox1\":\n", " plt.title(\"Sco X-1\")\n", " if point == \"sn1987a\":\n", " plt.title(\"SN 1987a\")\n", " if point == \"gc\":\n", " plt.title(\"Galactic Center\")\n", " plt.tight_layout()\n", " plt.savefig('NBR_sigma_' + point + '.png',bbox_inches='tight')\n", " plt.show()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotSigma(\"scox1\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotSigma(\"sn1987a\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plotSigma(\"gc\")" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "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.9.2" } }, "nbformat": 4, "nbformat_minor": 2 }