# 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')