Information taken from: "Evidence for Substructure in Rich Clusters of Galaxies from Radial-Velocity Measurements" by Alan Dressler and Stephen Shectman
For each galaxy with a measured radial velocity, we find the ten nearest neighbors with measured velocities. From this sample of 11, we calculate the local velocity mean and dispersion determined for the entire cluster sample. The global mean velocity $\bar{v}$ and dispersion $\sigma$ are calculated assuming a Gaussian distribution, but this assumption is ultimately removed by the normalizing procedure.
$$\delta^2 = (11/\sigma^2)[(\bar{v_{local}} - \bar{v})^2 + (\sigma_{local} - \sigma)^2]$$
In [1]:
import glob
from matplotlib import pyplot as plt
import numpy as np
import scipy.stats
import glob
from astropy.stats import biweight_location, biweight_midvariance, sigma_clip, bootstrap
from astropy import cosmology
from astropy.io import fits
from math import sqrt
%matplotlib inline
In [3]:
infile='biweight_center_scale.dat'
in1=np.recfromcsv(infile,delimiter='')
In [17]:
print in1.dtype.names
in1['name'] == 'NRGb155'
in1['center'][in1['name'] == 'NRGb155'][0]
Out[17]:
In [7]:
s='Abell2063'
print in1['name'] == s
In [22]:
def distance(ra,dec,all_ra,all_dec, n=10):
d = np.sqrt((ra-all_ra)**2+(dec-all_dec)**2)
t = np.argsort(d)
return t[0:n+1]
topfive=['NRGb155','NRGb226','Abell2063','NRGs117','NRGb247']
if __name__ == '__main__':
for cluster in topfive:
cl = fits.getdata(cluster + '_NSA.fits')
center_vel = in1['center'][in1['name'] == cluster][0]
sigma = in1['scale'][in1['name'] == cluster][0]
# cols = A2063_AGC[1].columns # Prints the titles of each column
# cols.info()
ds = np.zeros(len(cl.RA))
keepflag = abs(cl.ZDIST*3.e5 - center_vel) < 3.*sigma
ave_v=np.mean(cl.ZDIST[keepflag]*3.e5)
std_v=np.std(cl.ZDIST[keepflag]*3.e5)
#plt.figure()
#plt.plot(cl.RA,cl.DEC,'k.')
#for i in range(1):
for i in range(len(cl.RA)):
d_index = distance(cl.RA[i],cl.DEC[i],cl.RA,cl.DEC, n=10)
#print d_index
# calculate avg velocity
ave_local = np.mean(cl.ZDIST[d_index]*3.e5)
std_local = np.std(cl.ZDIST[d_index]*3.e5)
ds[i] = 11./std_v**2*((ave_local - ave_v)**2+ (std_local - std_v)**2)
# if i == 250:
# plt.plot(cl.RA[d_index],cl.DEC[d_index],'ro',markersize=10)
# plt.plot(cl.RA[i],cl.DEC[i],'bs')
# print ds[i]
# plt.scatter(cl.RA[i],cl.DEC[i],s=ds[i],color='k')
plt.figure()
plt.title(cluster)
plt.xlabel('RA')
plt.ylabel('DEC')
plt.scatter(cl.RA[keepflag],cl.DEC[keepflag],s=ds[keepflag],c=cl.ZDIST[keepflag]*3.e5)
plt.colorbar(label='Recession Velocity')
plt.figure()
plt.hist(cl.ZDIST*3e5, bins=np.arange(3000,10000,500),color='b')
plt.hist(cl.ZDIST[keepflag]*3e5, bins=np.arange(3000,10000,500),color='r')
plt.axvline(x=center_vel,color='k',ls='-')
plt.axvline(x=center_vel+3.*sigma,color='k',ls='--')
plt.axvline(x=center_vel-3.*sigma,color='k',ls='--')
s = '%5.1f, %5.1f'%(center_vel,sigma)
plt.text(0.1,0.9,s,transform=plt.gca().transAxes)
#plt.axis([230,232,7,10])
In [ ]:
In [ ]:
In [ ]: