In [ ]:

    

from math import *
from numpy import *
from scipy import integrate
from scipy.interpolate import InterpolatedUnivariateSpline
fg = open("lrg.ascii", "r")
fd = open("random.ascii", "r")
sd = open("bao.ascii", "w")


    

In [ ]:

    

# Values for the cosmic parameters
c = 3.0*10**5
Om = 0.3
H0 = 72.0
Ol = 0.7
rty = 2.0


    

In [ ]:

    

# Formula for radial diameter distance
func = lambda z: 1.0/sqrt(Om*(1.0 + z)**3 +0.7)
listyy = []
listss = []
ss = 0.01


    

In [ ]:

    

# Integrating the radial distance and forming the angular diameter distance with a spline of the radial distance
while ss <= rty:
      y, err = integrate.quad(func, 0.0, ss)
      listss.append(ss)
      listyy.append(y)
      ss = ss + 0.01
Hz = InterpolatedUnivariateSpline(listss, listyy)
def angdist(zz):
      value = c*(1.0 + zz)/H0 * Hz(zz)
      return(value)
listxD = []
listyD = []
listzD = []
weightD = []
listxR = []
listyR = []
listzR = []
weightR = []
a=1
b=2
n=0
k=0
for line in fg:
  red = float(line.split()[2])
  angle = float(line.split()[0])
  dec = float(line.split()[1])
  weight = float(line.split()[6])
  dist = angdist(red)
  # Converting declination into polar angle
  dd = (pi/2.0) - dec
  # Converting into spherical coordinates
  xx = dist*cos(angle)*sin(dd)
  yy = dist*sin(angle)*sin(dd)
  zz = dist*cos(dd)
  listxD.append(xx)
  listyD.append(yy)
  listzD.append(zz)
  weightD.append(weight)
  n=n+1
# As above but for the random catalogue
for line in fd:
  red = float(line.split()[2])
  angle = float(line.split()[0])
  dec = float(line.split()[1])
  weight = float(line.split()[6])
  dist = angdist(red)
  dd = (pi/2.0) - dec
  xx = dist*cos(angle)*sin(dd)
  yy = dist*sin(angle)*sin(dd)
  zz = dist*cos(dd)
  listxR.append(xx)
  listyR.append(yy)
  listzR.append(zz)
  weightR.append(weight)
  k=k+1

fuzzy = 0
compare = 100000
bins = 201
# This is the size of the bins


    

In [ ]:

    

size = 1.0
counter = 0
listD = bins*[0]
listR = bins*[0]
listDR = bins*[0]


    

In [ ]:

    

# To reduce computation time, instead of checking if distances are less than 200 Mpc, instead check off against 200^2. This prevents the code from unnecessarily doing a square root for every calculation.
comp = 40000.0
m=0
while m < n:
    x0 = listxD[m]
    y0 = listyD[m]
    z0 = listzD[m]
    w = weightD[m]
    # Setting the parameter space around the origin galaxy
    u = x0 + 200.0
    i = x0 - 200.0
    o = y0 + 200.0
    p = y0 - 200.0
    r = z0 + 200.0
    t = z0 - 200.0
    listxt = []
    listyt = []
    listzt = []
    wwt = []
    oo = 0
for j in range(m + 1, n):
  x1 = listxD[j]
  y1 = listyD[j]
  z1 = listzD[j]
  ww1 = weightD[j]
#Checking to see which galaxies are within the volume demarcated around the origin galaxy
  if i < x1 < u and p < y1 < o and t < z1 < r:
    listxt.append(x1)
    listyt.append(y1)
    listzt.append(z1)
    wwt.append(ww1)
    oo = oo + 1
for e in range(0, oo):
  x2 = listxt[e]
  y2 = listyt[e]
  z2 = listzt[e]
  ww2 = wwt[e]
  # Calculating the distance from the x, y, z coordinates
  dd = (x0 - x2)**2 + (y0 - y2)**2 + (z0 - z2)**2
  if dd <= comp:
      # Checking to see which bin the distance is to be assigned to
      ds = int(sqrt(dd)/size)
      io = ww2 * w
      listD[ds] = listD[ds] + io
      counter = counter + 1
      fuzzy = fuzzy + 1
      if fuzzy == compare:
          print fuzzy
          compare = compare + 100000
m=m+1


    

In [ ]:

    

# As above but now for the DR correlation v=0
while v < n:
  x0 = listxD[v]
  y0 = listyD[v]
  z0 = listzD[v]
  w = weightD[v]
  u = x0 + 200.0
  i = x0 - 200.0
  o = y0 + 200.0
  p = y0 - 200.0
  r = z0 + 200.0
  t = z0 - 200.0
  listxt = []
  listyt = []
  listzt = []
  wwt = []
  oo = 0
  for j in range(0, k):
    x1 = listxR[j]
    y1 = listyR[j]
    z1 = listzR[j]
    ww1 = weightR[j]
    if i < x1 < u and p < y1 < o and t < z1 < r:
      listxt.append(x1)
      listyt.append(y1)
      listzt.append(z1)
      wwt.append(ww1)
      oo = oo + 1
  for e in range(0, oo):
    x2 = listxt[e]
    y2 = listyt[e]
    z2 = listzt[e]
    ww2 = wwt[e]
# Calculating the distance from the x, y, z coordinates
    dd = (x0 - x2)**2 + (y0 - y2)**2 + (z0 - z2)**2
    if dd <= comp:
    # Checking to see which bin the distance is to be assigned to
      ds = int(sqrt(dd)/size)
      io = ww2 * w
      listDR[ds] = listDR[ds] + io
      counter = counter + 1
      fuzzy = fuzzy + 1
      if fuzzy == compare:
          print fuzzy
          compare = compare + 100000
v=v+1


    

In [ ]:

    

# As above for the RR correlation q=0
while q < k:
  x0 = listxR[q]
  y0 = listyR[q]
  z0 = listzR[q]
  w = weightR[q]
  u = x0 + 200.0
  i = x0 - 200.0
  o = y0 + 200.0
  p = y0 - 200.0
  r = z0 + 200.0
  t = z0 - 200.0
  listxt = []
  listyt = []
  listzt = []
  wwt = []
  oo = 0
  for j in range(q + 1, k):
    x1 = listxR[j]
    y1 = listyR[j]
    z1 = listzR[j]
    ww1 = weightR[j]
    if i < x1 < u and p < y1 < o and t < z1 < r:
      listxt.append(x1)
      listyt.append(y1)
      listzt.append(z1)
      wwt.append(ww1)
      oo = oo + 1
  for e in range(0, oo):
    x2 = listxt[e]
    y2 = listyt[e]
    z2 = listzt[e]
    ww2 = wwt[e]
    dd = (x0 - x2)**2 + (y0 - y2)**2 + (z0 - z2)**2
    if dd <= comp:
      ds = int(sqrt(dd)/size)
      io = ww2 * w
      listR[ds] = listR[ds] + io
      counter = counter + 1
      fuzzy = fuzzy + 1
      if fuzzy == compare:
        print fuzzy
        compare = compare + 100000
q=q+1
# Writing the DD, DR and RR bins to file where they will then be used to calculate the 2-point correlation function
for l in range(0, bins):
    xl = listD[l]
    fr = listDR[l]
    op = listR[l]
    er = l * size
    sd.write("%f %f %f %f\n" % (er, xl, fr, op))
print counter
fg.close()
fd.close()
sd.close()