Codes for making plots for O3 anisotropic stochastic paper

Importing packages and setting things up;

Before proceeding make sure that all the packages listed below are already installed

In [23]:
import numpy as np
import healpy as hp
import h5py
import scipy.io as sio
import matplotlib # requires matplotlib version >= 3.2.2
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
from scipy.io import loadmat
import matplotlib.patches as mpatches
from matplotlib.lines import Line2D
import warnings
warnings.filterwarnings("ignore", category=UserWarning) 

Parameters for plotting

In [24]:
fontsize = 24
matplotlib.rcParams.update({
    "font.size": fontsize,
    "axes.titlesize": fontsize,
    "axes.labelsize": fontsize,
    "xtick.labelsize": fontsize,
    "ytick.labelsize": fontsize,
    "xtick.major.size": fontsize * 0.8,
    "ytick.major.size": fontsize * 0.8,
    "legend.fontsize": fontsize * 0.9,
    "font.family": "Times new Roman",
    "figure.dpi": 100,
    "savefig.dpi": 300,
    "text.usetex": True,
    "path.simplify": True,
    "figure.figsize": (8, 6)
})

# Class to deal in an easy way with number of digits in colorbars (used for ULs and sigma maps)
# It allows to choose the number of digits for tickslabels and keeping 10^N at the right of the colorbar
# instead of having 1eN for each tickslabel.

class FormatScalarFormatter(matplotlib.ticker.ScalarFormatter):
            def __init__(self, fformat="%1.1f", offset=True, mathText=True):
                self.fformat = fformat
                matplotlib.ticker.ScalarFormatter.__init__(self,useOffset=offset,
                                                        useMathText=mathText)
            def _set_format(self):
                self.format = self.fformat
                if self._useMathText:
                    #self.format = '$%s$' % matplotlib.ticker._mathdefault(self.format)
                    self.format = '$%s$' % ('\\mathdefault{%s}' % self.format)

Variable definitions

In [25]:
# variable corresponding to spectral shapes of signal models 
alphas = ["0", "2over3", "3"]

Figure 2: Top row (BBR SNR maps)

In [26]:
# SNR sky-maps for BBR analysis

# Create figure
fig = plt.figure(figsize=(24, 8))

# Limits for the colorbars
vlims = [2.55, 3.72, 3.63]

# Loop over the *combined_maps.hdf5 files for the three spectral indices 0, 2/3 and 3
for alpha, vlim, index in zip (alphas, vlims, range(0,len(alphas))):
    
    # SNR *.hdf5 file
    snr_file = "./datafiles/bbr_a" + alpha + "_o1_o2_o3_combined_maps.hdf5"
    
    # 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)
    combined_snr = h5py.File(snr_file, "r")

    # Select the map of interest (SNR map)
    snr_map = np.array(combined_snr["map_snr_total"])
    print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))

    
    # Close the file
    combined_snr.close()
    
    # Plot the sky-map for a given alpha as subplot using healpy Mollweide projection
    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)
    
    # Add graticule
    hp.graticule(45,90)
    
    # Add useful annotations on the map for right ascension and declination
    plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) 
    
    # Trick to make colorbars all the same (healpy colorbar does not leave much freedom of customization and it is not used)
    fig = plt.gcf()
    ax = plt.gca()  
    image = ax.get_images()[0]
    cmap = fig.colorbar(image, ax=ax, shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40, ticks=[-2,0,2])
plt.show()
alpha:0, max snr:2.6
alpha:2over3, max snr:2.7
alpha:3, max snr:3.6

Figure 2: Bottom row (BBR Upper limit maps)

In [27]:
# Upper limit sky-maps for BBR analysis

fig = plt.figure(figsize=(24, 8))

# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for alpha, index in zip (alphas, range(0,len(alphas))):

    # upper limit *.mat file
    ul_dir = "./datafiles/bbr_a"+ alpha +"_zl_limits.mat"
    
    # 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
    combined_ul = sio.loadmat(ul_dir)
    ul_map = np.squeeze(np.array(combined_ul["flux_ul_combined"]))
    print("alpha:%s, min UL:%.1e, max UL:%.1e" % (alpha, np.min(ul_map), np.max(ul_map)))
    
    # Plot the sky-map for a given alpha as subplot using healpy Mollweide projection
    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])

    # Add graticule
    hp.graticule(45,90)
    
    # Add useful annotations on the map for right ascension and declination
    plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
    plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
    plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize) 
    
    
    fig = plt.gcf()
    ax = plt.gca()  
    image = ax.get_images()[0]
    
    # Further colorbar customization
    nticks = 4
    interval = 1./(nticks-1)
    diff= np.max(ul_map) -np.min(ul_map)
    ticks_no_round = [np.min(ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]
    
    cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=FormatScalarFormatter("%.1f"), shrink=0.7, orientation="horizontal", pad=0.15, aspect = 30)
    
    # Adding units as colorbar label
    cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20)
    cmap.ax.tick_params(labelsize=fontsize)
plt.show()
alpha:0, min UL:1.7e-08, max UL:7.6e-08
alpha:2over3, min UL:8.5e-09, max UL:4.1e-08
alpha:3, min UL:1.3e-10, max UL:1.1e-09

Figure 3: Top row (SHD SNR maps)

In [28]:
# Limits for the colorbars
vlims = [2.18, 3.72, 3.57]

fig = plt.figure(figsize=(24,8))
# Loop over the shd *.mat files for the three spectral indices 0, 2/3 and 3
for i, (alpha, vlim) in enumerate(zip(alphas, vlims)):
    ax = fig.add_subplot(130+(i+1), projection="mollweide")
    pproc_file = "./datafiles/shd_ul_a"+ alpha + "_o1_o2_o3_shd_map_files.mat"
    shd_data = loadmat(pproc_file)
    
    # Get the signal estimate (map0) and its sigma (sigmaPix0) for different pixels on the sky
    p0_map = shd_data["map0"]
    sigma_map = shd_data["sigmaPix0"]
    sig_map = sigma_map.reshape(360,181).T
    
    # Calculate SNR maps
    snr_map = (p0_map / sigma_map).reshape(360, 181).T
    print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))
    
    # Make the X, Y mesh and plot
    x = np.linspace(-np.pi,np.pi,360)
    y = np.linspace(-np.pi/2,np.pi/2,181)
    X,Y = np.meshgrid(x,y)
    image = ax.pcolormesh(X,Y,snr_map[:, ::-1], vmin=-vlim, vmax=vlim)
    
    # formate x and y tickes 
    raval = [-np.pi/2,0,np.pi/2]
    ralab = ['18h','12h','6h']
    ax.set_xticks(raval)
    xtick_labels = ax.set_xticklabels(ralab, color="w",fontweight='bold',fontsize=20)
    declval = [-np.pi/4,0,np.pi/4]
    decllab = [r'$-45^{\circ}$',r'$0^{\circ}$',r'$45^{\circ}$']
    ax.set_yticks(declval)
    ax.set_yticklabels(decllab,fontsize=20)
    plt.grid(True, color="k", linestyle=":")
    x = np.array([-np.radians(180), np.radians(180)])
    y = np.array([0, 0])
    plt.plot(x, y, ls='-', color="k", lw=0.8)
    cmp = plt.colorbar(image, shrink=0.7, ax=ax, orientation="horizontal", aspect=40)
    cmap.set_label('', fontsize = 10)
    for t in cmap.ax.get_xticklabels(): # reduce the font size of xlabel
        t.set_fontsize(10)
plt.show()
alpha:0, max snr:2.2
alpha:2over3, max snr:3.7
alpha:3, max snr:3.2

Figure 3: bottom row (SHD UL maps)

In [29]:
fig = plt.figure(figsize=(24,8))

# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for i, alpha in enumerate(alphas):
    ax = fig.add_subplot(130+(i+1), projection="mollweide")

    # upper limit *.mat file
    pproc_file = "./datafiles/shd_a" + alpha + "_limits_o1_o2_o3_HLV_ZL.mat"
    shd_data = loadmat(pproc_file)
    omega_ul_map = shd_data["omega_ul"].reshape(360, 181).T
    print("alpha:%s, min UL:%.1e, max UL:%.1e" % (alpha, np.min(omega_ul_map), np.max(omega_ul_map)))
    
    # Make the X, Y mesh and plot
    x = np.linspace(-np.pi,np.pi,360)
    y = np.linspace(-np.pi/2,np.pi/2,181)
    X,Y = np.meshgrid(x,y)
    image = ax.pcolormesh(X,Y,omega_ul_map[:, ::-1])

    # formate x and y tickes  
    raval = [-np.pi/2,0,np.pi/2]
    ralab = ['18h','12h','6h']
    ax.set_xticks(raval)
    xtick_labels = ax.set_xticklabels(ralab, color="w",fontweight='bold',fontsize=20)
    declval = [-np.pi/4,0,np.pi/4]
    decllab = ['$-45^{\circ}$','$0^{\circ}$','$45^{\circ}$']
    ax.set_yticks(declval)
    ax.set_yticklabels(decllab,fontsize=20)
    plt.grid(True, color="k", linestyle=":")
    x = np.array([-np.radians(180), np.radians(180)])
    y = np.array([0, 0])
    plt.plot(x, y, ls='-', color="k", lw=0.8)
    
    # Further colorbar customization
    nticks = 4
    interval = 1./(nticks-1)
    diff= np.max(omega_ul_map) -np.min(omega_ul_map)
    ticks_no_round = [np.min(omega_ul_map)+(diff)*interval*x for x in np.arange(0,nticks)]
    
    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)
    # Adding units as colorbar label
    cmap.set_label(r'$[\mathrm{sr}^{-1}]$', fontsize = fontsize,labelpad=20.)
    for t in cmap.ax.get_xticklabels():
        t.set_fontsize(20)
plt.show()
alpha:0, min UL:3.2e-09, max UL:9.3e-09
alpha:2over3, min UL:1.9e-09, max UL:9.7e-09
alpha:3, min UL:5.6e-10, max UL:3.4e-09

Figure 4: C_l upper limits

In [30]:
H0 = 2.2005e-18 # Hubble constant 
fref = 25.0 # reference frequency
const = 2.0 * np.pi**2 * fref**3 / (3 * H0**2) # normalization constant

colors = ['xkcd:black', 'xkcd:sky blue', 'xkcd:orange']
fig = plt.figure()
ax = fig.add_subplot(111)

# Loop over the upper limit *.mat files for the three spectral indices 0, 2/3 and 3
for alpha, color in zip(alphas, colors):
    pproc_file = "./datafiles/shd_ul_a" + alpha + "_o1_o2_o3_shd_map_files.mat"
    shd_data = loadmat(pproc_file)
    lmax = shd_data['L']
    Cell_ul = shd_data["omega_cl_ul"]
    Cell = shd_data['omega_cl']
    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", "/"))

ax.grid(which="major", linestyle="-")
ax.grid(which="minor", linestyle=":")
ax.legend(loc="upper right")
ax.set_xticks(np.arange(0, 17, 2))
ax.set_yticks([1e-8, 1e-9, 1e-10])
ax.set_ylim((0.8e-10, 0.2e-7))
ax.set_xlabel("$\ell$")
ax.set_ylabel("$C_\ell^{1/2}[\mathrm{sr}^{-1}]$")
plt.show()

Figure 5: Fisher Matrix comparison

In [31]:
baselines = ["HL", "HLV"]
colors = ['xkcd:black', 'xkcd:sky blue', 'xkcd:orange']
fkeep = 2/3 # keep 2/3 of dominant eigen values
yvals = []
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)

# Loop over the *.mat files for the three spectral indices 0, 2/3 and 3; for HL and HLV network
for i, (alpha, color) in enumerate(zip(alphas, colors)):    
    for ii, (baseline, style) in enumerate(zip(baselines, ["dotted", "solid"])):
        matfile = "./datafiles/shd_ul_a" + alpha + "_shd_map_files_" + baseline + ".mat"
        Fisher = loadmat(matfile)["fisher0"]
        eigs, _ = np.linalg.eigh(Fisher) # calculate eigen values
        eigs = eigs[::-1]
        yvals.append(eigs/eigs[0])
        if alpha == '2over3':
            lab = r'$\alpha=2/3$' + ', ' + baseline
        else:   
            lab = r'$\alpha=$%s'%alpha + ', ' + baseline
        ax.loglog(np.arange(len(eigs)), eigs/eigs[0], color=color, linestyle=style, label=lab)
    ax.axvline(len(eigs) * fkeep, linestyle="--", color=color)
_, _ = ax.get_legend_handles_labels()
legend_elements1 = [mpatches.Patch(color=color, label=r"$\alpha=%s$" % alpha.replace("over", "/")) for color, alpha in zip(colors, alphas)]
legend_elements2 = [Line2D([0], [0], color="k", lw=2, label=label, linestyle=ls) for ls, label in zip(["-", ":"], ["HLV", "HL"])]
leg1 = ax.legend(handles=legend_elements1, loc='lower left', frameon=True, fontsize=18)
leg2 = ax.legend(handles=legend_elements2, loc=(0.03, 0.3), frameon=True, fontsize=18)
ax.add_artist(leg1)
for line in leg2.get_lines():
    line.set_linewidth(2)
ax.set_xlabel("Mode index",fontsize=20)
ax.set_ylabel("Eigenvalues",fontsize=20)
plt.show()
findfont: Font family ['Times new Roman'] not found. Falling back to DejaVu Sans.

Figure 6: NBR plots

In [32]:
# function to NBR plot upper limits for the three directions of interest (Sco X-1, SN 1987a, Galactic Center )
def plotUL(point):
    # point = {scox1,sn1987a,gc}
    nbr_file = "./datafiles/NBR_" + point + "_ul.hdf5"
    struct = h5py.File(nbr_file, "r") # load the corresponding file
    freqs = np.array(struct["f_all"])
    ul = np.array(struct["ul_" + point])
    onesigma = np.array(struct["onesigma_" + point])
    
    plt.xlim([np.amin(freqs),np.amax(freqs)]);
    plt.ylim([5e-26,1e-22]);
    plt.grid(which = "both", lw = 0.2, ls = "--")
    plt.loglog(freqs, ul, c = "dimgrey", lw =  0.3, label = "HLV: O1+O2+O3: 95 % UL" )
    plt.loglog(freqs, onesigma, c = "k", lw = 0.5, label = "$1\sigma$ sensitivity")

    plt.ylabel("Strain amplitude ($h_0$)")
    plt.xlabel("Frequency [Hz]")
    plt.legend()
    if point == "scox1":
        plt.title("Sco X-1")
        plt.savefig('ul_and_onesig-scox1-HLV.png',bbox_inches='tight')
    if point == "sn1987a":
        plt.title("SN 1987a")
        plt.savefig('ul_and_onesig-sn-HLV.png',bbox_inches='tight')
    if point == "gc":
        plt.title("Galactic Center")
        plt.savefig('ul_and_onesig-gc-HLV.png',bbox_inches='tight')
    return "Plotting " + point
In [33]:
plotUL("scox1")
Out[33]:
'Plotting scox1'
In [34]:
plotUL("sn1987a")
Out[34]:
'Plotting sn1987a'
In [35]:
plotUL("gc")
Out[35]:
'Plotting gc'

Figure 7: Sigma maps for single baselines (BBR)

In [36]:
# HL, HV, LV 1-sigma maps

# Normalization constant such that the 1-sigma quantity corresponds to GW flux
G = 6.6726e-8
c = 2.9979e10
fref = 25.0
flux_factor = (c**3) *np.pi*(fref**2)/(4*G)

# Same as for ULs maps, but with 9 subplots of 1-sigma sensitivity for three baselines and three spectral indices
fig = plt.figure(figsize=(24,18))
baselines = ["HL", "HV", "LV"]
for ii, baseline in enumerate(baselines):
    print("Baseline:%s " % (baseline))
    for jj, alpha in enumerate (alphas): 
        map_dir = "./datafiles/bbr_a{0}_{1}_Map_dirty32.hdf5".format(alpha, baseline)
        maps = h5py.File(map_dir, "r") # load file
        print("alpha:%s, max snr:%.1f" % (alpha, np.max(np.array(maps['map_snr_total']))))        
        map_sigma = flux_factor*np.array(maps['sig_total']) # express 1-sigma in GW flux
        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])
        hp.graticule(45,90)
        plt.annotate('12h',(0.46,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
        plt.annotate('6h',(0.725,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
        plt.annotate('18h',(0.21,0.52),xycoords='axes fraction', color='w', fontsize=fontsize)
        plt.annotate(r'$45^{\circ}$',(0.03,0.8),xycoords='axes fraction', fontsize=fontsize)
        plt.annotate(r'$0^{\circ}$',(-0.045,0.47),xycoords='axes fraction', fontsize=fontsize)
        plt.annotate(r'$-45^{\circ}$',(-0.002,0.14),xycoords='axes fraction', fontsize=fontsize)  
        if alpha=='0':
            if baseline=='HL':
                plt.annotate('$\\bf{HL}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)
            if baseline=='HV':
                plt.annotate('$\\bf{HV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)
            if baseline=='LV':
                plt.annotate('$\\bf{LV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=30)
        
        fig = plt.gcf()
        ax = plt.gca()  
        image = ax.get_images()[0]
        nticks = 4
        interval = 1./(nticks-1)
        diff= np.max(map_sigma) -np.min(map_sigma)
        ticks_no_round = [np.min(map_sigma)+(diff)*interval*x for x in np.arange(0,nticks)]
        fmt = FormatScalarFormatter("%.1f")
        cmap = fig.colorbar(image, ax=ax, ticks=ticks_no_round, format=fmt, shrink=0.7, orientation="horizontal", pad=0.15, aspect = 40)
        cmap.set_label(r'[erg cm$^{-2}$ Hz$^{-1}$ s$^{-1}$]', fontsize = 1.2*fontsize, labelpad = 20)      
plt.show()
Baseline:HL 
alpha:0, max snr:2.3
alpha:2over3, max snr:2.5
alpha:3, max snr:3.7
Baseline:HV 
alpha:0, max snr:3.4
alpha:2over3, max snr:3.7
alpha:3, max snr:3.6
Baseline:LV 
alpha:0, max snr:3.1
alpha:2over3, max snr:3.1
alpha:3, max snr:4.1

Figure 8: Sigma maps for single baselines (SHD)

In [37]:
baselines = {'HL','HV','LV'}
fig = plt.figure(figsize=(24,18))

# Same as for ULs maps, but with 9 subplots of 1-sigma sensitivity for three baselines and three spectral indices
for ii, baseline in enumerate(baselines):
    print("Baseline:%s " % (baseline))
    for i, alpha in enumerate(alphas):
        ax = fig.add_subplot(330+(i+1)+3*ii, projection="mollweide")
        pproc_file = "./datafiles/shd_ul_a" + alpha + "_shd_map_files_" + baseline + ".mat"
        shd_data = loadmat(pproc_file) # load data file
        p0_map = shd_data["map0"]
        sigma_map = shd_data["sigmaPix0"]

        sig_map = sigma_map.reshape(360,181).T  * const # express 1-sigma as a normalized quantity
        snr_map = (p0_map / sigma_map).reshape(360, 181).T
        print("alpha:%s, max snr:%.1f" % (alpha, np.max(snr_map)))
        x = np.linspace(-np.pi,np.pi,360)
        y = np.linspace(-np.pi/2,np.pi/2,181)
        X,Y = np.meshgrid(x,y)

        # make plot
        image = ax.pcolormesh(X,Y,sig_map[:, ::-1], vmin=np.min(sig_map), vmax=np.max(sig_map))
        raval = [-np.pi/2,0,np.pi/2]
        ralab = ['18h','12h','6h']
        ax.set_xticks(raval)
        xtick_labels = ax.set_xticklabels(ralab, color="w",fontsize=20)
        declval = [-np.pi/4,0,np.pi/4]
        decllab = ['$-45^{\circ}$','$0^{\circ}$','$45^{\circ}$']
        ax.set_yticks(declval)
        ax.set_yticklabels(decllab,fontsize=20)
        plt.grid(True, color="k", linestyle=":")
        x = np.array([-np.radians(180), np.radians(180)])
        y = np.array([0, 0])
        plt.plot(x, y, ls='-', color="k", lw=0.8)
        
        nticks = 4
        interval = 1./(nticks-1)
        diff= np.max(sig_map) -np.min(sig_map)
        ticks_no_round = [np.min(sig_map)+(diff)*interval*x for x in np.arange(0,nticks)]

        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)
        cmap.set_label(r'$[\mathrm{sr}^{-1}]$', fontsize = fontsize, labelpad=20.)
        for t in cmap.ax.get_xticklabels():
            t.set_fontsize(20)
        if alpha=='0':
            if baseline=='HL':
                plt.annotate('$\\bf{HL}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)
            if baseline=='HV':
                plt.annotate('$\\bf{HV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)
            if baseline=='LV':
                plt.annotate('$\\bf{LV}$',(-0.18,0.45),xycoords='axes fraction',fontsize=24)
plt.show()
Baseline:HL 
alpha:0, max snr:1.6
alpha:2over3, max snr:3.0
alpha:3, max snr:3.9
Baseline:LV 
alpha:0, max snr:1.5
alpha:2over3, max snr:1.9
alpha:3, max snr:3.9
Baseline:HV 
alpha:0, max snr:2.1
alpha:2over3, max snr:3.9
alpha:3, max snr:4.0

Figure 9: Sigma maps for single baselines (NBR)

In [38]:
# function to NBR 1-sigma plots for the three directions of interest (Sco X-1, SN 1987a, Galactic Center) and three baselines
def plotSigma(point):
    # point = {scox1,sn1987a,gc}
    maps_mat_output = h5py.File("./datafiles/NBR_" + point + '_sigma_HLV.hdf5', "r")
    tempstruct = h5py.File("./datafiles/NBR_scox1_ul.hdf5", "r")
    f_all = np.array(tempstruct["f_all"])
    plt.loglog(f_all,np.array(maps_mat_output["sigma_HL"]), label='HL', lw =  2,color='black',ls = '-')
    plt.loglog(f_all, np.array(maps_mat_output["sigma_HV"]), label='HV',lw =  2,color='skyblue',ls = '-.')
    plt.loglog(f_all, np.array(maps_mat_output["sigma_LV"]), label='LV',lw =  2,color='orange',ls=':')
    plt.legend(bbox_to_anchor=(1,1), loc="upper right")
    plt.grid(which = "both", lw = 0.2, ls = "--")
    plt.xlabel('Frequency [Hz]')
    plt.ylabel('Sigma(f)')
    if point == "scox1":
        plt.title("Sco X-1")
    if point == "sn1987a":
        plt.title("SN 1987a")
    if point == "gc":
        plt.title("Galactic Center")
    plt.tight_layout()
    plt.savefig('NBR_sigma_' + point + '.png',bbox_inches='tight')
    plt.show()
In [39]:
plotSigma("scox1")
In [40]:
plotSigma("sn1987a")
In [41]:
plotSigma("gc")