# Examine or display the 3PC sensitivity map. # Reproduce 3PC Fig 25. # Obtain sensitivity limit for a given point on the sky. # Toby Burnett, July 2023. import matplotlib.pyplot as plt import healpy as hp filename='3PC_SensitivityMap_20230629.fits' # Generate plot similar to 3PC Figure 25. sensitivity_map = hp.read_map(filename) hp.mollview(sensitivity_map, norm='log', title=filename, unit=r'$\mathrm{G_{100}\ limit\ (erg\ cm^{-2}\ s^{-1})}$', xsize=1200, cmap='plasma', badcolor='white',min=4e-13, max=2e-12); plt.show() def get_lb(name): """ funcion which uses astropy to look up source names returns Galactic longitude, l, and latitude, b, of known source, in degrees. """ from astropy.coordinates import SkyCoord sc = SkyCoord.from_name(name).galactic return sc.l.deg, sc.b.deg def sensitivity(l,b): """ function that returns the sensitivity map value at (l,b): Sensitivity is the phase-integrated (that is, point source) 95% Confidence Level integral energy flux, integrated above 100 MeV, in units of erg/cm2/s. """ return sensitivity_map[hp.ang2pix(hp.get_nside(sensitivity_map),l,b, lonlat=True)] # demonstrations name = 'Vela' print(f'The source {name} is at l,b = {get_lb(name)}') print( f'The sensitivity at {name}, which should be Nan, since no measurement there):', sensitivity(*get_lb(name)) ) print( f'Sensitvity at (10,10) is {sensitivity (10,10):.2e} erg s-1 cm-2')