pro workshop3

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

readcol, 'random.dat', ra_r, dec_r

;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)
nd = double(n_elements(ral)) ;low 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
;print, DD: 195.00000       392.00000       1028.0000       2919.0000       8470.0000       23877.000   66917.000       188257.00       59511461.

;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
;print, RR:      151.00000       265.00000       825.00000       2341.0000       6915.0000       19405.000   54669.000       154986.00       49755443.

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_W3_RR.ps', encapsulated=0, times=1, isolatin1=1
plot, theta_begin[0:-2], RR[1:-2], psym=sym(1), xtitle='Theta', ytitle='RR', /xl, /yl
device, /close

;Normalize the DD and RR vectors
DDn = double(DD)/double(nd*(nd-1.0)/2.0)
RRn = double(RR)/double(nr*(nr-1.0)/2.0)

;Computing the correlation function (natural estimator)
omega_l = DDn/RRn - 1.0
;print, omega_l :    0.079586683      0.23663077     0.041691261     0.042395520     0.023979072     0.028646060     0.023281653     0.015449941  -9.2399322e-05

;plot
yr=[5e-3, 1.]
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_W3_ACF.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, /yl, xtitle='theta', ytitle='omega(theta)'


;Add error bars
normaliz = double(nd*(nd-1.0)/2.0)
err_omega_l = sqrt(DD)/normaliz/RRn

oploterror, theta_begin[0:-2], omega_l[1:-2], err_omega_l[1:-2], psym=sym(1), thick=4, color=0
device, /close

;Create DD2 for the high mass sample
DD2 = dblarr(bins+2)
nd2 = double(n_elements(rah)) ;high mass catalog

;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

DD2n = double(DD2)/double(nd2*(nd2-1.0)/2.0)
omega_h = DD2n/RRn - 1.0
;Add error bars
normaliz2 = double(nd2*(nd2-1.0)/2.0)
err_omega_h = sqrt(DD2)/normaliz2/RRn

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_W3_ACF_both.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, /yl, xtitle='theta', ytitle='omega(theta)'
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
device, /close





stop
END