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, 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, z, _ra, _dec
print("----")
except:
print("error! maybe can not identify from name")
print("----")
In [5]:
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)
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
comparison = density_inside/density_outside
return inside, outside, thetamin, thetamax, [number_inside, number_outside, area_inside, area_outside, density_inside, density_outside, comparison]
In [6]:
objlist = ['WISE J161021.87-395858.4',
'[HB89] 1741-038',
'3C 454.3',
'PKS 0539-057',
'PKS 0601-70',
'SSTSL2 J113006.83-144912.6',
'NGC 4945',
'[HB89] 0333+321 ABS01',
'PKS 0003-066',
'PKS 1830-21',
'[HB89] 0234+285',
'[HB89] 0748+126',
'[HB89] 1749+096',
'WISE J094857.31+002225.6',
'[HB89] 1104-445',
'3C 279',
'WISE J183005.92+061915.7',
'MESSIER 084',
'[HB89] 1514+197',
'LQAC 069+030 001',
'[HB89] 2131-021',
'4C +00.81',
'LQAC 066+023 001',
'PKS 1622-29',
'PKS 1057-79']
res = []
for objname in objlist:
try:
data, z, ra_center, dec_center = search_and_plot(objname, 2.0, 'fp_psc')
if len(data) >= Irsa.ROW_LIMIT: print("Probably number of rows limit < number of data")
inside, outside, thetamin, thetamax, result_par = check_overdensity(data, ra_center, dec_center, z, 0.5, 1.0)
print("Overdensity: ", result_par[-1])
res.append([objname, data, z, ra_center, dec_center, inside, outside, thetamin, thetamax, result_par])
print("-----")
except:
print("Error\n")
In [7]:
number_of_overdensity = 0
number_of_lowdensity = 0
for ires in res:
od = ires[-1][-1]
print(od)
if od > 1.0:
number_of_overdensity += 1
else:
number_of_lowdensity +=1
print("Low:", number_of_lowdensity, " High:", number_of_overdensity)
In [8]:
objlist = ['WISE J161021.87-395858.4',
'[HB89] 1741-038',
'3C 454.3',
'PKS 0539-057',
'PKS 0601-70',
'SSTSL2 J113006.83-144912.6',
'NGC 4945',
'[HB89] 0333+321 ABS01',
'PKS 0003-066',
'PKS 1830-21',
'[HB89] 0234+285',
'[HB89] 0748+126',
'[HB89] 1749+096',
'WISE J094857.31+002225.6',
'[HB89] 1104-445',
'3C 279',
'WISE J183005.92+061915.7',
'MESSIER 084',
'[HB89] 1514+197',
'LQAC 069+030 001',
'[HB89] 2131-021',
'4C +00.81',
'LQAC 066+023 001',
'PKS 1622-29',
'PKS 1057-79']
res = []
for objname in objlist:
try:
data, z, ra_center, dec_center = search_and_plot(objname, 2.0, 'fp_psc')
if len(data) >= Irsa.ROW_LIMIT: print("Probably number of rows limit < number of data")
inside, outside, thetamin, thetamax, result_par = check_overdensity(data, ra_center, dec_center, z, 1, 2.0)
print("Overdensity: ", result_par[-1])
res.append([objname, data, z, ra_center, dec_center, inside, outside, thetamin, thetamax, result_par])
print("-----")
except:
print("Error\n")
In [9]:
number_of_overdensity = 0
number_of_lowdensity = 0
for ires in res:
od = ires[-1][-1]
print(od)
if od > 1.0:
number_of_overdensity += 1
else:
number_of_lowdensity +=1
print("Low:", number_of_lowdensity, " High:", number_of_overdensity)
In [11]:
objlist = ['WISE J161021.87-395858.4',
'[HB89] 1741-038',
'3C 454.3',
'PKS 0539-057',
'PKS 0601-70',
'SSTSL2 J113006.83-144912.6',
'NGC 4945',
'[HB89] 0333+321 ABS01',
'PKS 0003-066',
'PKS 1830-21',
'[HB89] 0234+285',
'[HB89] 0748+126',
'[HB89] 1749+096',
'WISE J094857.31+002225.6',
'[HB89] 1104-445',
'3C 279',
'WISE J183005.92+061915.7',
'MESSIER 084',
'[HB89] 1514+197',
'LQAC 069+030 001',
'[HB89] 2131-021',
'4C +00.81',
'LQAC 066+023 001',
'PKS 1622-29',
'PKS 1057-79']
res = []
for objname in objlist:
try:
data, z, ra_center, dec_center = search_and_plot(objname, 2.0, 'fp_psc')
if len(data) >= Irsa.ROW_LIMIT: print("Probably number of rows limit < number of data")
inside, outside, thetamin, thetamax, result_par = check_overdensity(data, ra_center, dec_center, z, 0.5, 2.0)
print("Overdensity: ", result_par[-1])
res.append([objname, data, z, ra_center, dec_center, inside, outside, thetamin, thetamax, result_par])
print("-----")
except:
print("Error\n")
In [12]:
number_of_overdensity = 0
number_of_lowdensity = 0
for ires in res:
od = ires[-1][-1]
print(od)
if od > 1.0:
number_of_overdensity += 1
else:
number_of_lowdensity +=1
print("Low:", number_of_lowdensity, " High:", number_of_overdensity)
In [13]:
0.07030455605796143*3600
Out[13]:
In [ ]: