In [1]:
import sys
sys.path.append('../../src/utils/')
from galenv import *
from astroquery.irsa import Irsa
Irsa.ROW_LIMIT = 10000
%matplotlib inline
In [2]:
def plot_cone(coord, theta, res, xSize=7.5, ySize=7.5, title='', show=True, savefig=False, imgname="plot.png"):
'''Only cone
coord = astropy coordinates
theta = Cone angle
res = result catalog
'''
ra = coord.ra.value
dec = coord.dec.value
fig = plt.figure(figsize=(xSize, ySize))
gs = gridspec.GridSpec(1, 1)
ax = plt.subplot(gs[0])
# ax.axis('equal')
limangle = 1.5*theta
ax.set_xlim((ra-limangle, ra+limangle))
ax.set_ylim((dec-limangle, dec+limangle))
# Central position/object
ax.plot(ra, dec, 'ro', alpha=0.5)
# Catalog object
ax.plot(res['ra'], res['dec'], 'k.')
plt.gca().invert_xaxis() # RA from E to W
ax.set_xlabel('RA (deg)')
ax.set_ylabel('DEC (deg)')
plt.title(title)
# Circle
# it is wrong if I draw a circle around (ra, dec) with radius theta
# due to small circle in celestial sphere for DEC
circle = plt.Circle((ra, dec), theta, fc='none', ec='black')
ax.add_artist(circle)
fig.tight_layout()
if savefig:
plt.savefig(imgname)
if show:
plt.show()
plt.close()
In [3]:
ga = Galenv()
In [4]:
def search_and_plot(objname, ra, dec, tangential_dist, cat='fp_psc'):
try:
print(objname)
z, v0, _ra, _dec = ga.queryobject_byname(objname)
print("NED (z, v, ra, dec): ", z, v0, _ra, _dec)
obj_coord = coordinates.SkyCoord(ra=ra, dec=dec, unit=(u.deg, u.deg))
dA, theta = ga.calc_dA_theta(z, tangential_dist)
print("From redshift & tangential_dist (dA, theta):", dA, theta)
result = Irsa.query_region(obj_coord, catalog=cat, spatial="Cone", radius= theta * u.deg)
plot_cone(obj_coord, theta, result, savefig=True, imgname=objname + '.png')
return result
print("----")
except:
print("error! maybe can not identify from name")
print("----")
In [5]:
objname = 'PKS J2253+1608'
ra_center = 343.49061
dec_center = 16.148211
radius = 7.5 # in Mpc
z = 0.859
data = search_and_plot(objname, ra_center, dec_center, radius, 'fp_psc')
In [6]:
data
Out[6]:
In [7]:
def check_overdensity(data, ra_center, dec_center, z, rmin, rmax):
dAmin, thetamin = ga.calc_dA_theta(z, rmin)
dAmax, thetamax = ga.calc_dA_theta(z, rmax)
inside = []
outside = []
for idata in data:
ra_obj = idata['ra']
dec_obj = idata['dec']
obj = coordinates.SkyCoord(ra=ra_obj, dec=dec_obj, unit=(u.deg, u.deg))
center = coordinates.SkyCoord(ra=ra_center, dec=dec_center, unit=(u.deg, u.deg))
dist = center.separation(obj)
dist = dist.value
if dist < thetamin:
inside.append(idata)
if dist > thetamin and dist < thetamax:
outside.append(idata)
return inside, outside, thetamin, thetamax
In [8]:
inside, outside, thetamin, thetamax = check_overdensity(data, ra_center, dec_center, z, 0.5, 1.0)
In [9]:
number_inside = len(inside)
number_outside = len(outside)
area_inside = np.pi * thetamin**2
area_outside = np.pi * (thetamax**2 - thetamin**2)
density_inside = number_inside/area_inside
density_outside = number_outside/area_outside
In [10]:
print(number_inside, number_outside, area_inside, area_outside)
In [11]:
print(density_inside, density_outside)
print(density_inside/density_outside)
In [ ]: