pro workshop4

readcol, 'high.dat', rah, dech
readcol, 'low.dat', ral, decl

readcol, 'random.dat', ra_r, dec_r

;This was just to make my data to be a random catalog
;Creating the second random catalog
;ra_min = 0.0
;ra_max = 1.4
;dec_min= 0.0
;dec_max = 1.4
;ra_r2 = randomu(seed, 10000)* (ra_max-ra_min) + ra_min ; numpy.random.uniform()
;dec_r2 = randomu(seed, 10000)* (dec_max-dec_min) + dec_min
;openw, 1, 'random2.dat', width=1000
;for i=0, n_elements(ra_r2)-1 do printf, 1, ra_r2[i], dec_r2[i]
;close, 1
;ral =ra_r2
;decl = dec_r2

;Clustering parameters
theta_min= 5.d
theta_max= 200.0
bins = 7.d
nbins = bins + 1.0 
log_theta = alog10(theta_min) + (alog10(theta_max) -alog10(theta_min))*findgen(nbins)/float(nbins-1)
theta_begin=10.0^(log_theta)
bin_width1 = theta_begin - shift(theta_begin,1)
bin_width = bin_width1[1:-1]

;Create the DD
DD = dblarr(bins+2)
DD2 = dblarr(bins+2)
nd = double(n_elements(ral)) ;low mass catalog
nd2 = double(n_elements(rah)) ;high mass catalog

;DD
for ii = 0, nd-2 do begin
    GCIRC, 1, ral[ii]/15.0, decl[ii], ral[ii+1:-1]/15.0, decl[ii+1:-1], dis ;dist in arcsec
    pos =  floor(bins*(alog10(dis)-alog10(theta_begin[0]))/(alog10(theta_begin[-1])-alog10(theta_begin[0])) + 1.0) ; (log(x)-log(min))/binsize and binsize=(log(max)-log(min))/nbins. The +1 is because I am using first bin as additional one
    pos0 = where(pos lt 0)
    if (pos0[0] ne -1) then pos[pos0] = 0 ; if dist_proj <= Rmin then pos2=0
    pos1 = where(pos gt n_elements(theta_begin))
    if (pos1[0] ne -1) then pos[pos1] = n_elements(theta_begin) ; if dist_proj >Rmax then pos2=n_elements(R_begin) 
    for j=0, n_elements(pos)-1 do DD[pos[j]] = DD[pos[j]] + double(1.0) 
endfor

;DD2
for ii = 0, nd2-2 do begin
    GCIRC, 1, rah[ii]/15.0, dech[ii], rah[ii+1:-1]/15.0, dech[ii+1:-1], dis ;dist in arcsec
    pos =  floor(bins*(alog10(dis)-alog10(theta_begin[0]))/(alog10(theta_begin[-1])-alog10(theta_begin[0])) + 1.0) ; (log(x)-log(min))/binsize and binsize=(log(max)-log(min))/nbins. The +1 is because I am using first bin as additional one
    pos0 = where(pos lt 0)
    if (pos0[0] ne -1) then pos[pos0] = 0 ; if dist_proj <= Rmin then pos2=0
    pos1 = where(pos gt n_elements(theta_begin))
    if (pos1[0] ne -1) then pos[pos1] = n_elements(theta_begin) ; if dist_proj >Rmax then pos2=n_elements(R_begin) 
    for j=0, n_elements(pos)-1 do DD2[pos[j]] = DD2[pos[j]] + double(1.0) 
endfor

;Create the DR
DR = dblarr(bins+2)
DR2 = dblarr(bins+2)

;DR
for ii = 0, nd-2 do begin
    GCIRC, 1, ral[ii]/15.0, decl[ii], ra_r[0:-1]/15.0, dec_r[0:-1], dis ;dist in arcsec
    pos =  floor(bins*(alog10(dis)-alog10(theta_begin[0]))/(alog10(theta_begin[-1])-alog10(theta_begin[0])) + 1.0) ; (log(x)-log(min))/binsize and binsize=(log(max)-log(min))/nbins. The +1 is because I am using first bin as additional one
    pos0 = where(pos lt 0)
    if (pos0[0] ne -1) then pos[pos0] = 0 ; if dist_proj <= Rmin then pos2=0
    pos1 = where(pos gt n_elements(theta_begin))
    if (pos1[0] ne -1) then pos[pos1] = n_elements(theta_begin) ; if dist_proj >Rmax then pos2=n_elements(R_begin) 
    for j=0, n_elements(pos)-1 do DR[pos[j]] = DR[pos[j]] + double(1.0) 
endfor

;DR2
for ii = 0, nd2-2 do begin
    GCIRC, 1, rah[ii]/15.0, dech[ii], ra_r[0:-1]/15.0, dec_r[0:-1], dis ;dist in arcsec
    pos =  floor(bins*(alog10(dis)-alog10(theta_begin[0]))/(alog10(theta_begin[-1])-alog10(theta_begin[0])) + 1.0) ; (log(x)-log(min))/binsize and binsize=(log(max)-log(min))/nbins. The +1 is because I am using first bin as additional one
    pos0 = where(pos lt 0)
    if (pos0[0] ne -1) then pos[pos0] = 0 ; if dist_proj <= Rmin then pos2=0
    pos1 = where(pos gt n_elements(theta_begin))
    if (pos1[0] ne -1) then pos[pos1] = n_elements(theta_begin) ; if dist_proj >Rmax then pos2=n_elements(R_begin) 
    for j=0, n_elements(pos)-1 do DR2[pos[j]] = DR2[pos[j]] + double(1.0) 
endfor
;print, DR   308.00000       636.00000       1828.0000       5098.0000       14624.000       41970.000  120115.00       338762.00   1.0883666e+08
;print, DR2  360.00000       590.00000       1696.0000       5110.0000       14410.000       40833.000  117017.00       329266.00   1.0605072e+08

;Computation of RR
RR = dblarr(bins+2)
nr = double(n_elements(ra_r)) ;random catalog

;RR
for ii = 0, nr-2 do begin
    GCIRC, 1, ra_r[ii]/15.0, dec_r[ii], ra_r[ii+1:-1]/15.0, dec_r[ii+1:-1], dis ;dist in arcsec
    pos =  floor(bins*(alog10(dis)-alog10(theta_begin[0]))/(alog10(theta_begin[-1])-alog10(theta_begin[0])) + 1.0) ; (log(x)-log(min))/binsize and binsize=(log(max)-log(min))/nbins. The +1 is because I am using first bin as additional one
    pos0 = where(pos lt 0)
    if (pos0[0] ne -1) then pos[pos0] = 0 ; if dist_proj <= Rmin then pos2=0
    pos1 = where(pos gt n_elements(theta_begin))
    if (pos1[0] ne -1) then pos[pos1] = n_elements(theta_begin) ; if dist_proj >Rmax then pos2=n_elements(R_begin) 
    for j=0, n_elements(pos)-1 do RR[pos[j]] = RR[pos[j]] + double(1.0) 
endfor

;Normalize the DD and RR vectors
DDn = double(DD)/double(nd*(nd-1.0)/2.0)
DD2n = double(DD2)/double(nd2*(nd2-1.0)/2.0)
RRn = double(RR)/double(nr*(nr-1.0)/2.0)
DRn = double(DR)/double(nd*nr)
DR2n = double(DR2)/double(nd2*nr)

;Computing the correlation function (natural estimator)
omega_l = DDn/RRn - 1.0
omega_h = DD2n/RRn - 1.0
omega_l_LS = DDn/RRn -2.0*DRn/RRn + double(1.0)
omega_h_LS = DD2n/RRn -2.0*DR2n/RRn + double(1.0)


;Add error bars
normaliz = double(nd*(nd-1.0)/2.0)
normaliz2 = double(nd2*(nd2-1.0)/2.0)
err_omega_l = sqrt(DD)/normaliz/RRn
err_omega_h = sqrt(DD2)/normaliz2/RRn
err_omega_l_LS = (double(1.0)+omega_l_LS)/sqrt(double(DD))
err_omega_h_LS =(double(1.0)+omega_h_LS)/sqrt(double(DD2))

;plot DR
orig_device=!d.name 
@symbols_ps_kc
!p.font=0
set_plot,'ps'
Device, /color
TVLCT,255,0,0,100 ;red
TVLCT,0,255,0,200
TVLCT,0,0,255,300 ;blue
device,filename= 'plots_W4_DR.ps', encapsulated=0, times=1, isolatin1=1
plot, theta_begin[0:-2], DR[1:-2], psym=sym(1), xtitle='theta', ytitle='DR', /xl, /yl, title='Low mass sample'
device, /close

;plot
yr=[5e-3, 1.]
;yr=[-0.2, 0.2] ; for random-random ACF
xr= [2.0,theta_max]
orig_device=!d.name 
@symbols_ps_kc
!p.font=0
set_plot,'ps'
Device, /color
TVLCT,255,0,0,100 ;red
TVLCT,0,255,0,200
TVLCT,0,0,255,300 ;blue
;device,filename= 'plots_W4_ACF_random.ps', encapsulated=0, times=1, isolatin1=1
device,filename= 'plots_W4_ACF_estimators.ps', encapsulated=0, times=1, isolatin1=1
;plot, theta_begin[0:-2], omega_l[1:-2], psym=sym(1), xrange=xr, xstyle=1, yrange=yr, ystyle=1, /xl,  xtitle='theta', ytitle='omega(theta)'  ; for random-random ACF
plot, theta_begin[0:-2], omega_l[1:-2], psym=sym(1), xrange=xr, xstyle=1, yrange=yr, ystyle=1, /xl, /yl, xtitle='theta', ytitle='omega(theta)', title='Natural estimator'
oploterror, theta_begin[0:-2], omega_l[1:-2], err_omega_l[1:-2], psym=sym(1), thick=4, color=0
oploterror, theta_begin[0:-2], omega_h[1:-2], err_omega_h[1:-2], psym=sym(1), thick=4, color=100
plot, theta_begin[0:-2], omega_l_LS[1:-2], psym=sym(1), xrange=xr, xstyle=1, yrange=yr, ystyle=1, /xl, /yl, xtitle='theta', ytitle='omega(theta)', title='LS estimator'
oploterror, theta_begin[0:-2], omega_l_LS[1:-2], err_omega_l_LS[1:-2], psym=sym(1), thick=4, color=0
oploterror, theta_begin[0:-2], omega_h_LS[1:-2], err_omega_h_LS[1:-2], psym=sym(1), thick=4, color=100
device, /close

stop
END