Written by Rose A. Finn, updated on January 5, 2015

PURPOSE: 
  This program calculates the biweight center and scale for the LCS clusters
  using existing programs from the astropy.stats package.  Errors on the center
  and scale are calculating using bootstrap resampling (1000 samples is the default).


CALLING SEQUENCE

   from within ipython

   % run  ~/svnRepository/pythonCode/LCSbiweight.py
   getbiweightall()

   to see the results plotted with velocity histogram for one cluster

   mkw11.plotvhist()

   to see a multipanel plot for all clusters, type

   plotall()

INPUT PARAMETERS
  none

OUTPUT PARAMETERS
  none

EXAMPLES
  see calling sequence

PROCEDURE

  The NSA catalog is cut to include only those galaxies with velocities within
  4000 km/s of the central velocity (from the literature) and within a projected
  radius of 1 degree.  The biweight location and scale are calculated, and the
  location is used as the median, and the biweight location and scale are recalculated.
  This is repeated until the scale changes by less than 1 km/s.  This typically requires
  2 iterations.

REQUIRED PYTHON MODULES
  numpy
  pylab
  astropy
  scipy

In [1]:
from matplotlib import pyplot as plt
import numpy as np
import scipy.stats
from astropy.stats import biweight_location, biweight_midvariance, sigma_clip, bootstrap
from astropy import cosmology
from astropy.io import fits
import glob

%matplotlib inline

In [2]:
def centralbi(x, M=None):
    x=np.array(x,'f')
    if M is None:
        M=np.median(x)
    
    MAD=np.median(abs(x-M))
    ui=((x-M)/(6*MAD))
    top=np.sum((x-M)*((1-ui**2)**2))
    bottom=np.sum((1-ui**2)**2)
    
    
    cbi=M + (top/bottom)
    #print self.clustername
    # print(cbi)
    
    #finds the biweight scale
    n=len(x)
    usbi=((x-M)/(9*MAD))
    upper= np.sum(((x-M)**2)*((1-usbi**2)**4))
    lower=np.sum((1-usbi**2)*(1-5*usbi**2))
    sbi=np.sqrt(n)*((sqrt(upper))/(abs(lower)))
    return cbi,sbi
def getbiweight(z):
    scale_cut=3.
    biweightscale=biweight_midvariance(z)
    biweightlocation=biweight_location(z)

    flag=np.abs(z-biweightlocation)/biweightscale < scale_cut
    #flag=np.ones(len(z),'bool')
    repeatflag=1
    nloop=0
    #print biweightlocation, biweightscale
    oldbiweightscale=biweightscale
    while repeatflag:
        newdata=z[flag]
        biweightscale=biweight_midvariance(newdata, M=biweightlocation)
        biweightlocation=biweight_location(newdata, M=biweightlocation)
        oldflag = flag
        #flag=abs(z-biweightlocation)/biweightscale < scale_cut
        nloop += 1
        #print nloop, biweightlocation, biweightscale, len(newdata), sum(flag)
        #if sum(flag == oldflag) == len(flag): 
        #    repeatflag=0
        #    print 'nloop = ', nloop
        if np.abs(biweightscale - oldbiweightscale) < 1.: 
            repeatflag=0
            #print 'nloop = ', nloop
        if nloop > 5:
            repeatflag = 0
        oldbiweightscale=biweightscale
        #print nloop, biweightlocation, biweightscale
    #flag=abs(z-biweightlocation)/biweightscale < 4.
    #biweightscale=biweight_midvariance(z[flag],M=biweightlocation)
    return biweightlocation, biweightscale

In [3]:
# read in sample.dat - has cluster, RA, Dec, recession velocity
names=[]
RA=[]
DEC=[]
vr=[]
infile1=open('sample.dat','r')
for line in infile1:
    #print line
    t=line.split()
    names.append(t[0])
    RA.append(t[1])
    DEC.append(t[2])
    vr.append(t[3])
infile1.close()
vr=np.array(vr,'f')
RA=np.array(RA,'f')
DEC=np.array(DEC,'f')
names=np.array(names)
print names


['NRGb004' 'NRGs027' 'NRGs038' 'NRGs076' 'NRGs090' 'NRGs110' 'NRGs117'
 'NRGb128' 'NRGb155' 'NRGb177' 'NRGb226' 'NRGb244' 'NRGb247' 'NRGs317'
 'Abell2063']

In [6]:
# read in fits AGC table
infiles=glob.glob('*_NSA.fits')
outfile=open('biweight_center_scale.dat','w')
outfile.write('# name center scale RA DEC \n')
for f in infiles:
    clustername=f.split('_')[0]
    clustervel=vr[names == clustername]
    clusterRA=RA[names == clustername]
    clusterDEC=DEC[names == clustername]
    dat = fits.getdata(f)
    # cut velocities to within +/- 4000 km/s
    keepflag = (np.abs(dat.ZDIST*3.e5 - clustervel) < 4000.) & (np.sqrt((dat.RA-clusterRA)**2+(dat.DEC-clusterDEC)**2) < .75)
    # calculate biweight center and scale
    a,b=getbiweight(dat.ZDIST[keepflag]*3.e5)
    print clustername,": center vel = %5.1f, scale = %5.1f, RA = %12.8f, DEC = %12.8f"%(a,b,clusterRA,clusterDEC)
    outfile.write(clustername+" %5.1f %5.1f %12.8f %12.8f\n"%(a,b,clusterRA,clusterDEC))
outfile.close()

#print infiles
#dat.columns


Abell2063 : center vel = 10408.3, scale = 890.9, RA = 230.75779724, DEC =   8.63939953
NRGb004 : center vel = 8408.0, scale = 358.4, RA = 129.54791260, DEC =  25.11666679
NRGb128 : center vel = 7601.3, scale = 490.3, RA = 170.58999634, DEC =  24.32805634
NRGb155 : center vel = 6421.4, scale = 756.4, RA = 176.18583679, DEC =  19.69972229
NRGb177 : center vel = 6981.4, scale = 413.3, RA = 181.07417297, DEC =  20.25499916
NRGb226 : center vel = 7004.1, scale = 1048.3, RA = 194.94708252, DEC =  27.92833328
NRGb244 : center vel = 6919.8, scale = 287.0, RA = 201.04499817, DEC =  13.97972202
NRGb247 : center vel = 6794.6, scale = 374.9, RA = 202.38000488, DEC =  11.78861141
NRGs027 : center vel = 8544.7, scale = 417.9, RA = 139.02583313, DEC =  17.60222244
NRGs038 : center vel = 9342.5, scale = 857.3, RA = 140.89582825, DEC =  22.33277702
NRGs076 : center vel = 8825.3, scale = 476.0, RA = 151.67416382, DEC =  14.43027782
NRGs090 : center vel = 9737.7, scale = 306.5, RA = 158.25874329, DEC =  12.09472179
NRGs110 : center vel = 10444.4, scale = 725.6, RA = 165.21208191, DEC =  10.55472183
NRGs117 : center vel = 9780.0, scale = 742.8, RA = 167.63082886, DEC =  28.72750092
NRGs317 : center vel = 8864.0, scale = 314.3, RA = 221.79083252, DEC =  13.70638847

In [4]:


In [4]: