In [1]:
import healpix_util as hu
import astropy as ap
import math, sys, time
import numpy as np
from astropy.io import fits
from astropy.table import Table
from astropy.table import vstack
import astropy.io.ascii as ascii
from astropy.constants import c
import matplotlib.pyplot as plt
import math
import scipy.special as sp
from scipy import integrate
from astropy.stats import bayesian_blocks
from astropy.stats import histogram
from astropy.visualization import hist
%matplotlib inline

In [2]:
hmqdat1=ascii.read("./output/hmqbin1.dat")

In [3]:
hmqdat1


Out[3]:
<Table length=37885>
ZRADEC
float64float64float64
0.00111.887875-25.288639
0.001188.91611112.556111
0.001190.71833313.257222
0.001192.72055641.120278
0.001194.182521.683056
0.001197.94599-0.289915
0.001201.365583-43.018806
0.001202.46833347.194722
0.001202.58968244.688984
0.001204.983292-31.640083
.........
0.25258.5534242.34074
0.25258.86727954.419169
0.25259.53698331.494097
0.25267.562564.248889
0.25324.969722-43.252778
0.25325.47855810.752208
0.25339.553083-39.672028
0.25343.056389-8.963333
0.25345.75952-0.283827
0.25346.553822-0.986816

In [23]:
randra,randdec=hu.randsphere(num=len(hmqdat1['Z']),system='eq',ra_range=[min(hmqdat1['RA']),max(hmqdat1['RA'])],de_range=[min(hmqdat1['DEC']),max(hmqdat1['DEC'])])
randz=np.random.poisson(0.25,len(hmqdat1['Z']))

In [22]:
#hu.randsphere(num=400,system='eq',ra_range=[0,20],de_range=[0,360])
#[min(hmqdat1['DEC']),max(hmqdat1['DEC'])]=
hu.randsphere(num=len(hmqdat1['Z']),system='eq',ra_range=[0,20],de_range=[0,360])


Out[22]:
(array([  1.91980977,   7.59542376,   7.90522408, ...,  15.10234717,
         12.42486054,   2.36954371]),
 array([-28.04183662,  42.74128277,  -3.50169464, ...,  53.65479395,
         33.45348289,  77.33760118]))

In [4]:
help(np.random.poisson)


Help on built-in function poisson:

poisson(...)
    poisson(lam=1.0, size=None)
    
    Draw samples from a Poisson distribution.
    
    The Poisson distribution is the limit of the binomial distribution
    for large N.
    
    Parameters
    ----------
    lam : float or sequence of float
        Expectation of interval, should be >= 0. A sequence of expectation
        intervals must be broadcastable over the requested size.
    size : int or tuple of ints, optional
        Output shape.  If the given shape is, e.g., ``(m, n, k)``, then
        ``m * n * k`` samples are drawn.  Default is None, in which case a
        single value is returned.
    
    Returns
    -------
    samples : ndarray or scalar
        The drawn samples, of shape *size*, if it was provided.
    
    Notes
    -----
    The Poisson distribution
    
    .. math:: f(k; \lambda)=\frac{\lambda^k e^{-\lambda}}{k!}
    
    For events with an expected separation :math:`\lambda` the Poisson
    distribution :math:`f(k; \lambda)` describes the probability of
    :math:`k` events occurring within the observed
    interval :math:`\lambda`.
    
    Because the output is limited to the range of the C long type, a
    ValueError is raised when `lam` is within 10 sigma of the maximum
    representable value.
    
    References
    ----------
    .. [1] Weisstein, Eric W. "Poisson Distribution."
           From MathWorld--A Wolfram Web Resource.
           http://mathworld.wolfram.com/PoissonDistribution.html
    .. [2] Wikipedia, "Poisson distribution",
           http://en.wikipedia.org/wiki/Poisson_distribution
    
    Examples
    --------
    Draw samples from the distribution:
    
    >>> import numpy as np
    >>> s = np.random.poisson(5, 10000)
    
    Display histogram of the sample:
    
    >>> import matplotlib.pyplot as plt
    >>> count, bins, ignored = plt.hist(s, 14, normed=True)
    >>> plt.show()
    
    Draw each 100 values for lambda 100 and 500:
    
    >>> s = np.random.poisson(lam=(100., 500.), size=(100, 2))


In [ ]:


In [ ]:


In [8]:
rdat=ascii.read("./input/rdatf.dat")

In [5]:
hmqdat1


Out[5]:
<Table length=37885>
ZRADEC
float64float64float64
0.00111.887875-25.288639
0.001188.91611112.556111
0.001190.71833313.257222
0.001192.72055641.120278
0.001194.182521.683056
0.001197.94599-0.289915
0.001201.365583-43.018806
0.001202.46833347.194722
0.001202.58968244.688984
0.001204.983292-31.640083
.........
0.25258.5534242.34074
0.25258.86727954.419169
0.25259.53698331.494097
0.25267.562564.248889
0.25324.969722-43.252778
0.25325.47855810.752208
0.25339.553083-39.672028
0.25343.056389-8.963333
0.25345.75952-0.283827
0.25346.553822-0.986816

In [6]:
rdat


Out[6]:
<Table length=37885>
ZRADEC
float64float64float64
4.79121475247e-0672.377373403520.0512320354
5.16724118113e-06308.925925565-8.31255335055
1.89522769875e-05173.013645184-29.9838208208
1.92138087725e-0518.997778548420.678293975
2.57349485492e-05152.521111935.2451365534
2.63532549121e-0573.515340263624.7532901221
3.06301354178e-05120.893419183-13.7621959652
3.32600533562e-05322.55041759360.4296693751
4.31786092999e-05146.01388324-35.8467792924
4.67327098299e-05237.51372663352.0480118728
.........
0.249960006208116.379850093-53.0303303762
0.24996322218331.930717171-32.5528434175
0.249967582637259.68649147422.3017647038
0.249969846884333.555452555-35.6429474841
0.24997289713352.9619175-35.8820385318
0.2499758984879.717929980521.0700356321
0.249987342502265.701246713-18.6488828952
0.249988319448324.7062260675.7962410009
0.249992068468281.763179587-11.5775508903
0.2499935679169.416849142752.2534973327

In [ ]:
#Creating output file 
rdpair = open("./output/rdbin1pair_f.dat",'w')
# We intend to obtain average redshift of the pair and deltaZ deltatheta Zdeltatheta values first essentially finding right pairs
rdpair.write("Z\t deltaZ\t deltatheta\t Zdeltatheta\n")
#Algorithm keeps removing rows after done finding quasars within deltaZ˜0.02 
start_time=time.time()
while len(hmqdat1['Z'])>1:
    j=0
    #for j in range(0,len(rdat['Z'])):
    while rdat[j]['Z']<=(hmqdat1[0]['Z']+0.05):
        deltaZ=abs(rdat[j]['Z']-hmqdat1[0]['Z'])
        Z=(rdat[j]['Z']+hmqdat1[0]['Z'])/2.0
        RA1=hmqdat1[0]['RA']*math.pi/180.0
        RA2=rdat[j]['RA']*math.pi/180.0
        DEC1=hmqdat1[0]['DEC']*math.pi/180.0
        DEC2=rdat[j]['DEC']*math.pi/180.0
        deltatheta=abs(math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2)))
        Zdeltatheta=Z*deltatheta
        if Zdeltatheta<=0.05:
            rdpair.write("%f\t %f\t %f\t %f\n" %(Z, deltaZ, deltatheta, Zdeltatheta))
        j=j+1
    hmqdat1.remove_row(0)
rdpair.close()
print "Total time taken to run: "
print time.time()-start_time

In [ ]:
#Creating output file 
Ralpha = open("./output/Ralpha.dat",'w')
# We intend to obtain average redshift of the pair and deltaZ deltatheta Zdeltatheta R alpha values
Ralpha.write("Z\t deltaZ\t deltatheta\t Zdeltatheta\t R\t alpha\t \n")
#Algorithm keeps removing rows after done finding quasars within deltaZ˜0.02 
start_time=time.time()
while len(hmqdat['Z'])>0:
    j=1
    while hmqdat[j]['Z']<=(hmqdat[0]['Z']+0.02):
        deltaZ=abs(hmqdat[j]['Z']-hmqdat[0]['Z'])
        Z=(Z1+Z2)/2.0
        RA1=hmqdat[0]['RA']*math.pi/180.0
        RA2=hmqdat[j]['RA']*math.pi/180.0
        DEC1=hmqdat[0]['DEC']*math.pi/180.0
        DEC2=hmqdat[j]['DEC']*math.pi/180.0
        deltatheta=abs(math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2)))
        Zdeltatheta=Z*deltatheta
            if Zdeltatheta<=0.02:
                R=math.sqrt(deltaZ**2+Zdeltatheta**2)
                alpha=math.acos(deltaZ/R)
                Ralpha.write("%f\t %f\t %f\t %f\t %f\t %f\n" %(Z, deltaZ, deltatheta, Zdeltatheta, R, alpha))
        j=j+1
    hmqdat.remove_row(0)
Ralpha.close()
print "Total time taken to run: "
print time.time()-start_time

In [30]:
R=np.sqrt(pairdats['deltaZ']**2+(pairdats['Zdeltatheta'])**2)

In [31]:
R


Out[31]:
<Column name='deltaZ' dtype='float64' length=6523551>
0.00128606128869
0.00179004672537
0.0022292269767
0.0021520895157
0.002311934458
0.00208023630539
0.002285352415
0.0028272443826
0.00310334409638
0.0025454423944
0.00322552511666
0.003254806648
...
0.0198295971495
0.0141191705445
0.0134927469285
0.0206283616848
0.0117651921637
0.0185912557457
0.0173086205675
0.0176120554639
0.0144317748444
0.010431996
0.0227972339893
0.021606718955

In [34]:
Rdat=pairdats.Column(name='R',data=R)
pairdats.add_column(Rdat)

In [36]:
alpha=np.arccos(pairdats['deltaZ']/pairdats['R'])


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: invalid value encountered in divide
  if __name__ == '__main__':

In [38]:
alphadat=pairdats.Column(name='alpha',data=alpha)
pairdats.add_column(alphadat)

In [43]:
om=0.3
ol=0.7
xLCDM=1.0/np.sqrt(om*(1+pairdats['Z'])**3+ol)

In [42]:
xLCDMdat=pairdats.Column(name='xLCDM',data=xLCDM)
pairdats.add_column(xLCDMdat)
pairdats


Out[42]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDM
float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.00067467016
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.00067467016
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.00044985352
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.00044985352
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.00044985352
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.00044985352
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.00044985352
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.00044985352
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.00044985352
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.00044985352
.....................
3.57550.0050.0035050.012530.01349274692851.191174265940.184312767301
3.5780.0160.0036390.013020.02062836168480.6830818949810.184165396292
3.57950.0110.0011660.0041750.01176519216370.3626462156480.184077066095
3.5820.0140.0034150.0122320.01859125574570.7181228114610.183930002874
3.58550.0070.0044150.015830.01730862056751.154448749250.18372443672
3.5960.0160.0020470.007360.01761205546390.4311909431630.183109985927
3.5960.0060.003650.0131240.01443177484441.14202982280.183109985927
3.5960.00.0029010.0104340.0104319961.570796326790.183109985927
3.59850.0170.0042210.0151910.02279723398930.7292047403450.182964183113
3.60.0180.003320.0119520.0216067189550.5861541768170.182876792494

In [93]:
yLCDM0=np.sqrt(om*(1+pairdats['Z'])**3+ol)/pairdats['Z']


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:1: RuntimeWarning: divide by zero encountered in divide
  if __name__ == '__main__':

In [95]:
LCDMf = lambda x: 1.0/math.sqrt(om*(1+x)**3+ol)
np.vectorize(LCDMf)

def LCDMfint(z):
    return integrate.quad(LCDMf, 0, z)

LCDMfint=np.vectorize(LCDMfint)

yLCDM=yLCDM0*LCDMfint(pairdats['Z'])[0]


/Users/rohin/anaconda/lib/python2.7/site-packages/ipykernel/__main__.py:9: RuntimeWarning: invalid value encountered in multiply

In [99]:
yLCDMdat=pairdats.Column(name='yLCDM',data=yLCDM)
pairdats.add_column(yLCDMdat)
pairdats


Out[99]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDMyLCDM
float64float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.000674670160.999662947356
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.000674670160.999662947356
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.000449853520.9997751988
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.000449853520.9997751988
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.000449853520.9997751988
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.000449853520.9997751988
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.000449853520.9997751988
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.000449853520.9997751988
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.000449853520.9997751988
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.000449853520.9997751988
........................
3.57550.0050.0035050.012530.01349274692851.191174265940.1843127673012.42932784764
3.5780.0160.0036390.013020.02062836168480.6830818949810.1841653962922.43027204961
3.57950.0110.0011660.0041750.01176519216370.3626462156480.1840770660952.43083847858
3.5820.0140.0034150.0122320.01859125574570.7181228114610.1839300028742.4317823732
3.58550.0070.0044150.015830.01730862056751.154448749250.183724436722.43310350314
3.5960.0160.0020470.007360.01761205546390.4311909431630.1831099859272.43706463723
3.5960.0060.003650.0131240.01443177484441.14202982280.1831099859272.43706463723
3.5960.00.0029010.0104340.0104319961.570796326790.1831099859272.43706463723
3.59850.0170.0042210.0151910.02279723398930.7292047403450.1829641831132.43800726613
3.60.0180.003320.0119520.0216067189550.5861541768170.1828767924942.43857275155

In [100]:
sLCDM=pairdats['xLCDM']*np.sqrt((pairdats['deltaZ'])**2+(pairdats['yLCDM']*pairdats['Zdeltatheta'])**2)

In [101]:
sLCDMdat=pairdats.Column(name='sLCDM',data=sLCDM)
pairdats.add_column(sLCDMdat)
pairdats


Out[101]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDMyLCDMsLCDM
float64float64float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.000674670160.9996629473560.0012869625095
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.000674670160.9996629473560.00179110666594
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.000449853520.99977519880.0022303051794
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.000449853520.99977519880.00215311611837
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.000449853520.99977519880.00231296456834
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.000449853520.99977519880.00208108886284
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.000449853520.99977519880.00228634775052
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.000449853520.99977519880.00282796744178
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.000449853520.99977519880.00310439404409
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.000449853520.99977519880.00254663322885
...........................
3.57550.0050.0035050.012530.01349274692851.191174265940.1843127673012.429327847640.00568556885124
3.5780.0160.0036390.013020.02062836168480.6830818949810.1841653962922.430272049610.00653002077534
3.57950.0110.0011660.0041750.01176519216370.3626462156480.1840770660952.430838478580.00275499566675
3.5820.0140.0034150.0122320.01859125574570.7181228114610.1839300028742.43178237320.00604679069142
3.58550.0070.0044150.015830.01730862056751.154448749250.183724436722.433103503140.00719225310651
3.5960.0160.0020470.007360.01761205546390.4311909431630.1831099859272.437064637230.00440122912702
3.5960.0060.003650.0131240.01443177484441.14202982280.1831099859272.437064637230.00595875619834
3.5960.00.0029010.0104340.0104319961.570796326790.1831099859272.437064637230.00465618159246
3.59850.0170.0042210.0151910.02279723398930.7292047403450.1829641831132.438007266130.00745598271651
3.60.0180.003320.0119520.0216067189550.5861541768170.1828767924942.438572751550.00626464175486

In [104]:
ascii.write(pairdats,"./output/sLCDM_correl.dat")

In [105]:
pairdats


Out[105]:
<Table length=6523551>
ZdeltaZdeltathetaZdeltathetaRalphaxLCDMyLCDMsLCDM
float64float64float64float64float64float64float64float64float64
-0.00150.0010.539116-0.0008090.001286061288690.6800076331861.000674670160.9996629473560.0012869625095
-0.00150.0010.989785-0.0014850.001790046725370.9780455646811.000674670160.9996629473560.00179110666594
-0.0010.0020.984608-0.0009850.00222922697670.4574718758971.000449853520.99977519880.0022303051794
-0.0010.0020.794663-0.0007950.00215208951570.3782038314041.000449853520.99977519880.00215311611837
-0.0010.0021.159759-0.001160.0023119344580.5254936212221.000449853520.99977519880.00231296456834
-0.0010.0020.572174-0.0005720.002080236305390.2786442099751.000449853520.99977519880.00208108886284
-0.0010.0021.105819-0.0011060.0022853524150.5050742474061.000449853520.99977519880.00228634775052
-0.0010.0021.998327-0.0019980.00282724438260.7849797384161.000449853520.99977519880.00282796744178
-0.0010.0022.372919-0.0023730.003103344096380.8704715209171.000449853520.99977519880.00310439404409
-0.0010.0021.574572-0.0015750.00254544239440.6669402811261.000449853520.99977519880.00254663322885
...........................
3.57550.0050.0035050.012530.01349274692851.191174265940.1843127673012.429327847640.00568556885124
3.5780.0160.0036390.013020.02062836168480.6830818949810.1841653962922.430272049610.00653002077534
3.57950.0110.0011660.0041750.01176519216370.3626462156480.1840770660952.430838478580.00275499566675
3.5820.0140.0034150.0122320.01859125574570.7181228114610.1839300028742.43178237320.00604679069142
3.58550.0070.0044150.015830.01730862056751.154448749250.183724436722.433103503140.00719225310651
3.5960.0160.0020470.007360.01761205546390.4311909431630.1831099859272.437064637230.00440122912702
3.5960.0060.003650.0131240.01443177484441.14202982280.1831099859272.437064637230.00595875619834
3.5960.00.0029010.0104340.0104319961.570796326790.1831099859272.437064637230.00465618159246
3.59850.0170.0042210.0151910.02279723398930.7292047403450.1829641831132.438007266130.00745598271651
3.60.0180.003320.0119520.0216067189550.5861541768170.1828767924942.438572751550.00626464175486

In [ ]:
spairs_knuth=hist(pairdats['sLCDM'][2253:len(pairdats['sLCDM']-1)],bins='knuth')
spairs_blocks=hist(pairdats['sLCDM'][2253:len(pairdats['sLCDM']-1)],bins='blocks')
DDsk=AutoCorrelation(spairs_knuth[0])
DDsb=AutoCorrelation(spairs_blocks[0])
plt.plot(spairs_knuth[1][1:len(spairs_knuth[1])],DDsk)
plt.savefig("./images/DD_sLCDM_knuth.pdf")
plt.savefig("./images/DD_sLCDM_knuth.jpeg")
plt.plot(spairs_blocks[1][1:len(spairs_blocks[1])],DDsb)
plt.savefig("./images/DD_sLCDM_blocks.pdf")
plt.savefig("./images/DD_sLCDM_blocks.jpeg")

In [129]:
spairs_1000=hist(sLCDM[2253:len(sLCDM)-1],bins=1000)
plt.savefig("./images/hist_sLCDM_1000.pdf")
plt.savefig("./images/hist_sLCDM_1000.jpeg")
DDs1000=AutoCorrelation(spairs_1000[0])
plt.figure()
plt.plot(spairs_1000[1][1:len(spairs_1000[1])],DDs1000)
plt.savefig("./images/DD_sLCDM_1000.pdf")
plt.savefig("./images/DD_sLCDM_1000.jpeg")



In [112]:
max(sLCDM)


Out[112]:
0.028187966602223606

In [114]:
np.all(np.isfinite(sLCDM))


Out[114]:
False

In [115]:
min(sLCDM)


Out[115]:
9.6817160036876236e-07

In [116]:
help(np.isfinite)


Help on ufunc object:

isfinite = class ufunc(__builtin__.object)
 |  Functions that operate element by element on whole arrays.
 |  
 |  To see the documentation for a specific ufunc, use np.info().  For
 |  example, np.info(np.sin).  Because ufuncs are written in C
 |  (for speed) and linked into Python with NumPy's ufunc facility,
 |  Python's help() function finds this page whenever help() is called
 |  on a ufunc.
 |  
 |  A detailed explanation of ufuncs can be found in the "ufuncs.rst"
 |  file in the NumPy reference guide.
 |  
 |  Unary ufuncs:
 |  =============
 |  
 |  op(X, out=None)
 |  Apply op to X elementwise
 |  
 |  Parameters
 |  ----------
 |  X : array_like
 |      Input array.
 |  out : array_like
 |      An array to store the output. Must be the same shape as `X`.
 |  
 |  Returns
 |  -------
 |  r : array_like
 |      `r` will have the same shape as `X`; if out is provided, `r`
 |      will be equal to out.
 |  
 |  Binary ufuncs:
 |  ==============
 |  
 |  op(X, Y, out=None)
 |  Apply `op` to `X` and `Y` elementwise. May "broadcast" to make
 |  the shapes of `X` and `Y` congruent.
 |  
 |  The broadcasting rules are:
 |  
 |  * Dimensions of length 1 may be prepended to either array.
 |  * Arrays may be repeated along dimensions of length 1.
 |  
 |  Parameters
 |  ----------
 |  X : array_like
 |      First input array.
 |  Y : array_like
 |      Second input array.
 |  out : array_like
 |      An array to store the output. Must be the same shape as the
 |      output would have.
 |  
 |  Returns
 |  -------
 |  r : array_like
 |      The return value; if out is provided, `r` will be equal to out.
 |  
 |  Methods defined here:
 |  
 |  __call__(...)
 |      x.__call__(...) <==> x(...)
 |  
 |  __repr__(...)
 |      x.__repr__() <==> repr(x)
 |  
 |  __str__(...)
 |      x.__str__() <==> str(x)
 |  
 |  accumulate(...)
 |      accumulate(array, axis=0, dtype=None, out=None)
 |      
 |      Accumulate the result of applying the operator to all elements.
 |      
 |      For a one-dimensional array, accumulate produces results equivalent to::
 |      
 |        r = np.empty(len(A))
 |        t = op.identity        # op = the ufunc being applied to A's  elements
 |        for i in range(len(A)):
 |            t = op(t, A[i])
 |            r[i] = t
 |        return r
 |      
 |      For example, add.accumulate() is equivalent to np.cumsum().
 |      
 |      For a multi-dimensional array, accumulate is applied along only one
 |      axis (axis zero by default; see Examples below) so repeated use is
 |      necessary if one wants to accumulate over multiple axes.
 |      
 |      Parameters
 |      ----------
 |      array : array_like
 |          The array to act on.
 |      axis : int, optional
 |          The axis along which to apply the accumulation; default is zero.
 |      dtype : data-type code, optional
 |          The data-type used to represent the intermediate results. Defaults
 |          to the data-type of the output array if such is provided, or the
 |          the data-type of the input array if no output array is provided.
 |      out : ndarray, optional
 |          A location into which the result is stored. If not provided a
 |          freshly-allocated array is returned.
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          The accumulated values. If `out` was supplied, `r` is a reference to
 |          `out`.
 |      
 |      Examples
 |      --------
 |      1-D array examples:
 |      
 |      >>> np.add.accumulate([2, 3, 5])
 |      array([ 2,  5, 10])
 |      >>> np.multiply.accumulate([2, 3, 5])
 |      array([ 2,  6, 30])
 |      
 |      2-D array examples:
 |      
 |      >>> I = np.eye(2)
 |      >>> I
 |      array([[ 1.,  0.],
 |             [ 0.,  1.]])
 |      
 |      Accumulate along axis 0 (rows), down columns:
 |      
 |      >>> np.add.accumulate(I, 0)
 |      array([[ 1.,  0.],
 |             [ 1.,  1.]])
 |      >>> np.add.accumulate(I) # no axis specified = axis zero
 |      array([[ 1.,  0.],
 |             [ 1.,  1.]])
 |      
 |      Accumulate along axis 1 (columns), through rows:
 |      
 |      >>> np.add.accumulate(I, 1)
 |      array([[ 1.,  1.],
 |             [ 0.,  1.]])
 |  
 |  at(...)
 |      at(a, indices, b=None)
 |      
 |      Performs unbuffered in place operation on operand 'a' for elements
 |      specified by 'indices'. For addition ufunc, this method is equivalent to
 |      `a[indices] += b`, except that results are accumulated for elements that
 |      are indexed more than once. For example, `a[[0,0]] += 1` will only
 |      increment the first element once because of buffering, whereas
 |      `add.at(a, [0,0], 1)` will increment the first element twice.
 |      
 |      .. versionadded:: 1.8.0
 |      
 |      Parameters
 |      ----------
 |      a : array_like
 |          The array to perform in place operation on.
 |      indices : array_like or tuple
 |          Array like index object or slice object for indexing into first
 |          operand. If first operand has multiple dimensions, indices can be a
 |          tuple of array like index objects or slice objects.
 |      b : array_like
 |          Second operand for ufuncs requiring two operands. Operand must be
 |          broadcastable over first operand after indexing or slicing.
 |      
 |      Examples
 |      --------
 |      Set items 0 and 1 to their negative values:
 |      
 |      >>> a = np.array([1, 2, 3, 4])
 |      >>> np.negative.at(a, [0, 1])
 |      >>> print(a)
 |      array([-1, -2, 3, 4])
 |      
 |      ::
 |      
 |      Increment items 0 and 1, and increment item 2 twice:
 |      
 |      >>> a = np.array([1, 2, 3, 4])
 |      >>> np.add.at(a, [0, 1, 2, 2], 1)
 |      >>> print(a)
 |      array([2, 3, 5, 4])
 |      
 |      ::
 |      
 |      Add items 0 and 1 in first array to second array,
 |      and store results in first array:
 |      
 |      >>> a = np.array([1, 2, 3, 4])
 |      >>> b = np.array([1, 2])
 |      >>> np.add.at(a, [0, 1], b)
 |      >>> print(a)
 |      array([2, 4, 3, 4])
 |  
 |  outer(...)
 |      outer(A, B)
 |      
 |      Apply the ufunc `op` to all pairs (a, b) with a in `A` and b in `B`.
 |      
 |      Let ``M = A.ndim``, ``N = B.ndim``. Then the result, `C`, of
 |      ``op.outer(A, B)`` is an array of dimension M + N such that:
 |      
 |      .. math:: C[i_0, ..., i_{M-1}, j_0, ..., j_{N-1}] =
 |         op(A[i_0, ..., i_{M-1}], B[j_0, ..., j_{N-1}])
 |      
 |      For `A` and `B` one-dimensional, this is equivalent to::
 |      
 |        r = empty(len(A),len(B))
 |        for i in range(len(A)):
 |            for j in range(len(B)):
 |                r[i,j] = op(A[i], B[j]) # op = ufunc in question
 |      
 |      Parameters
 |      ----------
 |      A : array_like
 |          First array
 |      B : array_like
 |          Second array
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          Output array
 |      
 |      See Also
 |      --------
 |      numpy.outer
 |      
 |      Examples
 |      --------
 |      >>> np.multiply.outer([1, 2, 3], [4, 5, 6])
 |      array([[ 4,  5,  6],
 |             [ 8, 10, 12],
 |             [12, 15, 18]])
 |      
 |      A multi-dimensional example:
 |      
 |      >>> A = np.array([[1, 2, 3], [4, 5, 6]])
 |      >>> A.shape
 |      (2, 3)
 |      >>> B = np.array([[1, 2, 3, 4]])
 |      >>> B.shape
 |      (1, 4)
 |      >>> C = np.multiply.outer(A, B)
 |      >>> C.shape; C
 |      (2, 3, 1, 4)
 |      array([[[[ 1,  2,  3,  4]],
 |              [[ 2,  4,  6,  8]],
 |              [[ 3,  6,  9, 12]]],
 |             [[[ 4,  8, 12, 16]],
 |              [[ 5, 10, 15, 20]],
 |              [[ 6, 12, 18, 24]]]])
 |  
 |  reduce(...)
 |      reduce(a, axis=0, dtype=None, out=None, keepdims=False)
 |      
 |      Reduces `a`'s dimension by one, by applying ufunc along one axis.
 |      
 |      Let :math:`a.shape = (N_0, ..., N_i, ..., N_{M-1})`.  Then
 |      :math:`ufunc.reduce(a, axis=i)[k_0, ..,k_{i-1}, k_{i+1}, .., k_{M-1}]` =
 |      the result of iterating `j` over :math:`range(N_i)`, cumulatively applying
 |      ufunc to each :math:`a[k_0, ..,k_{i-1}, j, k_{i+1}, .., k_{M-1}]`.
 |      For a one-dimensional array, reduce produces results equivalent to:
 |      ::
 |      
 |       r = op.identity # op = ufunc
 |       for i in range(len(A)):
 |         r = op(r, A[i])
 |       return r
 |      
 |      For example, add.reduce() is equivalent to sum().
 |      
 |      Parameters
 |      ----------
 |      a : array_like
 |          The array to act on.
 |      axis : None or int or tuple of ints, optional
 |          Axis or axes along which a reduction is performed.
 |          The default (`axis` = 0) is perform a reduction over the first
 |          dimension of the input array. `axis` may be negative, in
 |          which case it counts from the last to the first axis.
 |      
 |          .. versionadded:: 1.7.0
 |      
 |          If this is `None`, a reduction is performed over all the axes.
 |          If this is a tuple of ints, a reduction is performed on multiple
 |          axes, instead of a single axis or all the axes as before.
 |      
 |          For operations which are either not commutative or not associative,
 |          doing a reduction over multiple axes is not well-defined. The
 |          ufuncs do not currently raise an exception in this case, but will
 |          likely do so in the future.
 |      dtype : data-type code, optional
 |          The type used to represent the intermediate results. Defaults
 |          to the data-type of the output array if this is provided, or
 |          the data-type of the input array if no output array is provided.
 |      out : ndarray, optional
 |          A location into which the result is stored. If not provided, a
 |          freshly-allocated array is returned.
 |      keepdims : bool, optional
 |          If this is set to True, the axes which are reduced are left
 |          in the result as dimensions with size one. With this option,
 |          the result will broadcast correctly against the original `arr`.
 |      
 |          .. versionadded:: 1.7.0
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          The reduced array. If `out` was supplied, `r` is a reference to it.
 |      
 |      Examples
 |      --------
 |      >>> np.multiply.reduce([2,3,5])
 |      30
 |      
 |      A multi-dimensional array example:
 |      
 |      >>> X = np.arange(8).reshape((2,2,2))
 |      >>> X
 |      array([[[0, 1],
 |              [2, 3]],
 |             [[4, 5],
 |              [6, 7]]])
 |      >>> np.add.reduce(X, 0)
 |      array([[ 4,  6],
 |             [ 8, 10]])
 |      >>> np.add.reduce(X) # confirm: default axis value is 0
 |      array([[ 4,  6],
 |             [ 8, 10]])
 |      >>> np.add.reduce(X, 1)
 |      array([[ 2,  4],
 |             [10, 12]])
 |      >>> np.add.reduce(X, 2)
 |      array([[ 1,  5],
 |             [ 9, 13]])
 |  
 |  reduceat(...)
 |      reduceat(a, indices, axis=0, dtype=None, out=None)
 |      
 |      Performs a (local) reduce with specified slices over a single axis.
 |      
 |      For i in ``range(len(indices))``, `reduceat` computes
 |      ``ufunc.reduce(a[indices[i]:indices[i+1]])``, which becomes the i-th
 |      generalized "row" parallel to `axis` in the final result (i.e., in a
 |      2-D array, for example, if `axis = 0`, it becomes the i-th row, but if
 |      `axis = 1`, it becomes the i-th column).  There are three exceptions to this:
 |      
 |      * when ``i = len(indices) - 1`` (so for the last index),
 |        ``indices[i+1] = a.shape[axis]``.
 |      * if ``indices[i] >= indices[i + 1]``, the i-th generalized "row" is
 |        simply ``a[indices[i]]``.
 |      * if ``indices[i] >= len(a)`` or ``indices[i] < 0``, an error is raised.
 |      
 |      The shape of the output depends on the size of `indices`, and may be
 |      larger than `a` (this happens if ``len(indices) > a.shape[axis]``).
 |      
 |      Parameters
 |      ----------
 |      a : array_like
 |          The array to act on.
 |      indices : array_like
 |          Paired indices, comma separated (not colon), specifying slices to
 |          reduce.
 |      axis : int, optional
 |          The axis along which to apply the reduceat.
 |      dtype : data-type code, optional
 |          The type used to represent the intermediate results. Defaults
 |          to the data type of the output array if this is provided, or
 |          the data type of the input array if no output array is provided.
 |      out : ndarray, optional
 |          A location into which the result is stored. If not provided a
 |          freshly-allocated array is returned.
 |      
 |      Returns
 |      -------
 |      r : ndarray
 |          The reduced values. If `out` was supplied, `r` is a reference to
 |          `out`.
 |      
 |      Notes
 |      -----
 |      A descriptive example:
 |      
 |      If `a` is 1-D, the function `ufunc.accumulate(a)` is the same as
 |      ``ufunc.reduceat(a, indices)[::2]`` where `indices` is
 |      ``range(len(array) - 1)`` with a zero placed
 |      in every other element:
 |      ``indices = zeros(2 * len(a) - 1)``, ``indices[1::2] = range(1, len(a))``.
 |      
 |      Don't be fooled by this attribute's name: `reduceat(a)` is not
 |      necessarily smaller than `a`.
 |      
 |      Examples
 |      --------
 |      To take the running sum of four successive values:
 |      
 |      >>> np.add.reduceat(np.arange(8),[0,4, 1,5, 2,6, 3,7])[::2]
 |      array([ 6, 10, 14, 18])
 |      
 |      A 2-D example:
 |      
 |      >>> x = np.linspace(0, 15, 16).reshape(4,4)
 |      >>> x
 |      array([[  0.,   1.,   2.,   3.],
 |             [  4.,   5.,   6.,   7.],
 |             [  8.,   9.,  10.,  11.],
 |             [ 12.,  13.,  14.,  15.]])
 |      
 |      ::
 |      
 |       # reduce such that the result has the following five rows:
 |       # [row1 + row2 + row3]
 |       # [row4]
 |       # [row2]
 |       # [row3]
 |       # [row1 + row2 + row3 + row4]
 |      
 |      >>> np.add.reduceat(x, [0, 3, 1, 2, 0])
 |      array([[ 12.,  15.,  18.,  21.],
 |             [ 12.,  13.,  14.,  15.],
 |             [  4.,   5.,   6.,   7.],
 |             [  8.,   9.,  10.,  11.],
 |             [ 24.,  28.,  32.,  36.]])
 |      
 |      ::
 |      
 |       # reduce such that result has the following two columns:
 |       # [col1 * col2 * col3, col4]
 |      
 |      >>> np.multiply.reduceat(x, [0, 3], 1)
 |      array([[    0.,     3.],
 |             [  120.,     7.],
 |             [  720.,    11.],
 |             [ 2184.,    15.]])
 |  
 |  ----------------------------------------------------------------------
 |  Data descriptors defined here:
 |  
 |  identity
 |      The identity value.
 |      
 |      Data attribute containing the identity element for the ufunc, if it has one.
 |      If it does not, the attribute value is None.
 |      
 |      Examples
 |      --------
 |      >>> np.add.identity
 |      0
 |      >>> np.multiply.identity
 |      1
 |      >>> np.power.identity
 |      1
 |      >>> print(np.exp.identity)
 |      None
 |  
 |  nargs
 |      The number of arguments.
 |      
 |      Data attribute containing the number of arguments the ufunc takes, including
 |      optional ones.
 |      
 |      Notes
 |      -----
 |      Typically this value will be one more than what you might expect because all
 |      ufuncs take  the optional "out" argument.
 |      
 |      Examples
 |      --------
 |      >>> np.add.nargs
 |      3
 |      >>> np.multiply.nargs
 |      3
 |      >>> np.power.nargs
 |      3
 |      >>> np.exp.nargs
 |      2
 |  
 |  nin
 |      The number of inputs.
 |      
 |      Data attribute containing the number of arguments the ufunc treats as input.
 |      
 |      Examples
 |      --------
 |      >>> np.add.nin
 |      2
 |      >>> np.multiply.nin
 |      2
 |      >>> np.power.nin
 |      2
 |      >>> np.exp.nin
 |      1
 |  
 |  nout
 |      The number of outputs.
 |      
 |      Data attribute containing the number of arguments the ufunc treats as output.
 |      
 |      Notes
 |      -----
 |      Since all ufuncs can take output arguments, this will always be (at least) 1.
 |      
 |      Examples
 |      --------
 |      >>> np.add.nout
 |      1
 |      >>> np.multiply.nout
 |      1
 |      >>> np.power.nout
 |      1
 |      >>> np.exp.nout
 |      1
 |  
 |  ntypes
 |      The number of types.
 |      
 |      The number of numerical NumPy types - of which there are 18 total - on which
 |      the ufunc can operate.
 |      
 |      See Also
 |      --------
 |      numpy.ufunc.types
 |      
 |      Examples
 |      --------
 |      >>> np.add.ntypes
 |      18
 |      >>> np.multiply.ntypes
 |      18
 |      >>> np.power.ntypes
 |      17
 |      >>> np.exp.ntypes
 |      7
 |      >>> np.remainder.ntypes
 |      14
 |  
 |  signature
 |  
 |  types
 |      Returns a list with types grouped input->output.
 |      
 |      Data attribute listing the data-type "Domain-Range" groupings the ufunc can
 |      deliver. The data-types are given using the character codes.
 |      
 |      See Also
 |      --------
 |      numpy.ufunc.ntypes
 |      
 |      Examples
 |      --------
 |      >>> np.add.types
 |      ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
 |      'LL->L', 'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D',
 |      'GG->G', 'OO->O']
 |      
 |      >>> np.multiply.types
 |      ['??->?', 'bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l',
 |      'LL->L', 'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D',
 |      'GG->G', 'OO->O']
 |      
 |      >>> np.power.types
 |      ['bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L',
 |      'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'FF->F', 'DD->D', 'GG->G',
 |      'OO->O']
 |      
 |      >>> np.exp.types
 |      ['f->f', 'd->d', 'g->g', 'F->F', 'D->D', 'G->G', 'O->O']
 |      
 |      >>> np.remainder.types
 |      ['bb->b', 'BB->B', 'hh->h', 'HH->H', 'ii->i', 'II->I', 'll->l', 'LL->L',
 |      'qq->q', 'QQ->Q', 'ff->f', 'dd->d', 'gg->g', 'OO->O']


In [117]:
np.argwhere(np.isnan(sLCDM))


Out[117]:
array([[ 206],
       [ 207],
       [ 208],
       ..., 
       [2250],
       [2251],
       [2252]])

In [119]:
sLCDM[2252]


Out[119]:
nan

In [121]:
pairdats['Z'][2252]


Out[121]:
0.0

In [ ]:
spairs1000=spairs_knuth

In [134]:
plt.figure()
plt.loglog(4.1e3*spairs_1000[1][1:len(spairs_1000[1])],DDs1000)
plt.savefig("./images/DD_sLCDMlog_1000.pdf")
plt.savefig("./images/DD_sLCDMlog_1000.jpeg")



In [135]:
plt.figure()
plt.loglog(4.1e3*spairs_1000[1][1:len(spairs_1000[1])],spairs_1000[1][1:len(spairs_1000[1])]**2*DDs1000)
plt.savefig("./images/DD_sLCDMlog_1000.pdf")
plt.savefig("./images/DD_sLCDMlog_1000.jpeg")



In [54]:
np.where(hmqdat['Z']<=0.25)


Out[54]:
(array([    0,     1,     2, ..., 37950, 37951, 37952]),)

In [55]:
np.where(hmqdat['Z']>0)


Out[55]:
(array([    67,     68,     69, ..., 510352, 510353, 510354]),)

In [56]:
hmqbin1=hmqdat[67:37952]

In [57]:
hmqbin1


Out[57]:
<Table length=37885>
ZRADEC
float64float64float64
0.00111.887875-25.288639
0.001188.91611112.556111
0.001190.71833313.257222
0.001192.72055641.120278
0.001194.182521.683056
0.001197.94599-0.289915
0.001201.365583-43.018806
0.001202.46833347.194722
0.001202.58968244.688984
0.001204.983292-31.640083
.........
0.25258.5534242.34074
0.25258.86727954.419169
0.25259.53698331.494097
0.25267.562564.248889
0.25324.969722-43.252778
0.25325.47855810.752208
0.25339.553083-39.672028
0.25343.056389-8.963333
0.25345.75952-0.283827
0.25346.553822-0.986816

In [58]:
ascii.write(hmqbin1,"./output/hmqbin1.dat")

In [59]:
len(hmqbin1)


Out[59]:
37885

In [ ]:


In [60]:
randra,randdec=hu.randsphere(37885)
randz=0.25*np.random.random_sample(37885)

In [61]:
randz


Out[61]:
array([ 0.09055629,  0.17966635,  0.19842346, ...,  0.00283044,
        0.14986337,  0.1608593 ])

In [62]:
randra


Out[62]:
array([  82.83397398,   16.6656554 ,  165.51338515, ...,  115.29707449,
        191.65893917,  258.64839955])

In [63]:
randdec


Out[63]:
array([-31.4041433 , -32.79569862,  22.5723584 , ..., -28.97524559,
        77.49756571, -28.649747  ])

In [64]:
rdat = Table([randz, randra, randdec], names=('Z', 'RA', 'DEC'))

In [65]:
ascii.write(rdat,"./input/rdatf.dat")

In [66]:
rdat


Out[66]:
<Table length=37885>
ZRADEC
float64float64float64
0.090556294765882.8339739799-31.4041432973
0.17966635020116.6656554042-32.7956986173
0.198423456837165.51338515422.5723583964
0.0269945747649.8786430626-47.9838130232
0.249567027275189.926592265-14.1726643577
0.20805363853410.1352114864-63.9908692299
0.123940268516347.577168388-46.005290915
0.0208993081512321.130083571-49.7373804336
0.0597156114309330.07736732818.4620279676
0.20782996978255.8054223472-32.5626952012
.........
0.22630613076272.347986703-32.6878124709
0.17417361379236.947505738-15.7671372385
0.0508709212651281.209876081-34.3132613188
0.199287293881329.683055775-80.2396491292
0.0213769230744208.34418634242.9564082458
0.0902636684931286.865628061-32.0319052474
0.127189969569314.31914424829.6982196894
0.00283043855431115.297074488-28.9752455884
0.149863371893191.65893916777.4975657083
0.160859297555258.648399553-28.6497470047

In [67]:
rdat.sort('Z')

In [69]:
ascii.write(rdat,"./input/rdatf.dat",overwrite=True)

In [68]:
rdat


Out[68]:
<Table length=37885>
ZRADEC
float64float64float64
4.79121475247e-0672.377373403520.0512320354
5.16724118113e-06308.925925565-8.31255335055
1.89522769875e-05173.013645184-29.9838208208
1.92138087725e-0518.997778548420.678293975
2.57349485492e-05152.521111935.2451365534
2.63532549121e-0573.515340263624.7532901221
3.06301354178e-05120.893419183-13.7621959652
3.32600533562e-05322.55041759360.4296693751
4.31786092999e-05146.01388324-35.8467792924
4.67327098299e-05237.51372663352.0480118728
.........
0.249960006208116.379850093-53.0303303762
0.24996322218331.930717171-32.5528434175
0.249967582637259.68649147422.3017647038
0.249969846884333.555452555-35.6429474841
0.24997289713352.9619175-35.8820385318
0.2499758984879.717929980521.0700356321
0.249987342502265.701246713-18.6488828952
0.249988319448324.7062260675.7962410009
0.249992068468281.763179587-11.5775508903
0.2499935679169.416849142752.2534973327

In [78]:
#Creating output file 
ddpair = open("./output/ddbin1pair.dat",'w')
# We intend to obtain average redshift of the pair and deltaZ deltatheta Zdeltatheta values first essentially finding right pairs
ddpair.write("Z\t deltaZ\t deltatheta\t Zdeltatheta\n")
#Algorithm keeps removing rows after done finding quasars within deltaZ˜0.02 
start_time=time.time()
while len(hmqbin1['Z'])>1:
    j=1
    while hmqbin1[j]['Z']<=(hmqbin1[0]['Z']+0.05):
        deltaZ=abs(hmqbin1[j]['Z']-hmqbin1[0]['Z'])
        Z=(hmqbin1[j]['Z']+hmqbin1[0]['Z'])/2.0
        RA1=hmqbin1[0]['RA']*math.pi/180.0
        RA2=hmqbin1[j]['RA']*math.pi/180.0
        DEC1=hmqbin1[0]['DEC']*math.pi/180.0
        DEC2=hmqbin1[j]['DEC']*math.pi/180.0
        deltatheta=abs(math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2)))
        Zdeltatheta=Z*deltatheta
        if Zdeltatheta<=0.05:
            ddpair.write("%f\t %f\t %f\t %f\n" %(Z, deltaZ, deltatheta, Zdeltatheta))
        j=j+1
    hmqbin1.remove_row(0)
ddpair.close()
print "Total time taken to run: "
print time.time()-start_time


---------------------------------------------------------------------------
IndexError                                Traceback (most recent call last)
<ipython-input-78-a60bd2d973f0> in <module>()
      7 while len(hmqbin1['Z'])>1:
      8     j=1
----> 9     while hmqbin1[j]['Z']<=(hmqbin1[0]['Z']+0.05):
     10         deltaZ=abs(hmqbin1[j]['Z']-hmqbin1[0]['Z'])
     11         Z=(hmqbin1[j]['Z']+hmqbin1[0]['Z'])/2.0

/Users/rohin/anaconda/lib/python2.7/site-packages/astropy/table/table.pyc in __getitem__(self, item)
   1196             return self.columns[item]
   1197         elif isinstance(item, (int, np.integer)):
-> 1198             return self.Row(self, item)
   1199         elif (isinstance(item, (tuple, list)) and item and
   1200               all(isinstance(x, six.string_types) for x in item)):

/Users/rohin/anaconda/lib/python2.7/site-packages/astropy/table/row.pyc in __init__(self, table, index)
     41         if index < -n or index >= n:
     42             raise IndexError('index {0} out of range for table with length {1}'
---> 43                              .format(index, len(table)))
     44 
     45     def __getitem__(self, item):

IndexError: index 5744 out of range for table with length 5744

In [70]:
from astropy.table import vstack

In [74]:
rdraw=vstack([hmqbin1,rdat])

In [76]:
rdraw.sort('Z')

In [77]:
rdraw


Out[77]:
<Table length=75770>
ZRADEC
float64float64float64
4.79121475247e-0672.377373403520.0512320354
5.16724118113e-06308.925925565-8.31255335055
1.89522769875e-05173.013645184-29.9838208208
1.92138087725e-0518.997778548420.678293975
2.57349485492e-05152.521111935.2451365534
2.63532549121e-0573.515340263624.7532901221
3.06301354178e-05120.893419183-13.7621959652
3.32600533562e-05322.55041759360.4296693751
4.31786092999e-05146.01388324-35.8467792924
4.67327098299e-05237.51372663352.0480118728
.........
0.25203.1799137.63222
0.25204.40517131.270606
0.25204.70186445.77337
0.25205.06027833.4125
0.25207.2225567.946605
0.25208.7147729.425905
0.25210.4144442.95594
0.25210.8421172.375816
0.25198.65487712.805779
0.25346.553822-0.986816

In [80]:
hmqdat1=ascii.read("./output/hmqbin1.dat")

In [ ]:
hu.

In [ ]:
#Creating output file 
rrpair = open("./output/rrbin1pair.dat",'w')
# We intend to obtain average redshift of the pair and deltaZ deltatheta Zdeltatheta values first essentially finding right pairs
rrpair.write("Z\t deltaZ\t deltatheta\t Zdeltatheta\n")
#Algorithm keeps removing rows after done finding quasars within deltaZ˜0.02 
start_time=time.time()
while len(rdat['Z'])>1:
    j=1
    while rdat[j]['Z']<=(rdat[0]['Z']+0.05):
        deltaZ=abs(rdat[j]['Z']-rdat[0]['Z'])
        Z=(rdat[j]['Z']+rdat[0]['Z'])/2.0
        RA1=rdat[0]['RA']*math.pi/180.0
        RA2=rdat[j]['RA']*math.pi/180.0
        DEC1=rdat[0]['DEC']*math.pi/180.0
        DEC2=rdat[j]['DEC']*math.pi/180.0
        deltatheta=abs(math.acos(math.sin(DEC1)*math.sin(DEC2)+math.cos(DEC1)*math.cos(DEC2)*math.cos(RA1-RA2)))
        Zdeltatheta=Z*deltatheta
        if Zdeltatheta<=0.05:
            rrpair.write("%f\t %f\t %f\t %f\n" %(Z, deltaZ, deltatheta, Zdeltatheta))
        j=j+1
    rdat.remove_row(0)
rrpair.close()
print "Total time taken to run: "
print time.time()-start_time