pro workshop2

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

;Creating the random catalog
ra_min = 0.0
ra_max = 1.4
dec_min= 0.0
dec_max = 1.4

ra_r = randomu(seed, 10000)* (ra_max-ra_min) + ra_min ; numpy.random.uniform()
dec_r = randomu(seed, 10000)* (dec_max-dec_min) + dec_min

openw, 1, 'random.dat', width=1000
for i=0, n_elements(ra_r)-1 do printf, 1, ra_r[i], dec_r[i]
close, 1

;Plot random catalog
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_W2_random.ps', encapsulated=0, times=1, isolatin1=1
plot, ra_r, dec_r, psym=sym(1), symsize=0.3, xtitle='RA[deg]', ytitle='DEC[deg]', title='random catalog'
device, /close

;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)
print, 'theta_begin', theta_begin
;theta_begin       5.0000000       8.4690697       14.345028       24.297809       41.155967       69.710551118.07670       199.99997
bin_width1 = theta_begin - shift(theta_begin,1)
bin_width = bin_width1[1:-1]
print, 'bin width', bin_width
;bin width       3.4690697       5.8759586       9.9527806       16.858158       28.554584       48.366152 81.923262

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

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

;DD One option (inefficient)
;for i=0, nd-2 do begin
;    for j=i+1, nd-1 do begin
        ;compute the distance between the galaxy with positions ral[i], decl[i] and the galaxy with positions ral[j], decl[j]
        ;check what is the position in theta_begin in which we should store that distance, and call that position "pos"
        ;increase one in DD at that position: DD[pos] = DD[pos]+1.0
;    endfor
;endfor

;DD One better option
for ii=0, nd-2 do begin
    ;compute the distance between the galaxy with positions ral[i], decl[i] and all the rest of galaxies below that position: ral[i+1:-1], decl[i+1:-1]
    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
    
    ;This part is only becasue I am storing the number pairs at lower and higher values as theta_beg and theta_end
    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

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






END