pro histo1_polynear,NNB=NNB,p,Nmin,Nmax,d,e

loadct,39

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Reading the catalogs and defining variables;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
readcol2,'./polynear/output.etbl',w_id,w_poly
readcol2, '../test_carlos/results_2011/results_NNB'+strcompress(strtrim(NNB))+'/wideal-wcarlos',Id,w_ideal,wcarlos
readcol2,'../cata/CMC_i24_test_polynear.cat',v1,v2,v3,v4,v5,v6,v7,v8

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Tests;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
wideal_wpoly=[transpose(w_ideal),transpose(w_poly)]
openw,5,'./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/wideal_wpoly'
printf,5,wideal_wpoly
close,5

help,w_poly
help,w_ideal
a=n_elements(w_poly)
b=n_elements(w_ideal)
c=n_elements(w_id)
IF (a EQ b) THEN BEGIN
print, 'w_ideal and w_poly have the same amount of elements, (cool!)',a
ENDIF ELSE BEGIN
print,'w_ideal and w_poly do not have the same number of elements (bug!)'
stop
ENDELSE

IF (b EQ c) THEN BEGIN
print, 'w_ideal and w_id have the same amount of elements,and then have a good chance to be equal (cool!)'
ENDIF ELSE BEGIN
print,'w_ideal and w_id do not have the same amount of elements,there is no way they are equal (bug!)'
stop
ENDELSE
d=WHERE(w_id eq w_ideal)
IF n_elements(d) EQ c THEN BEGIN
print,'w_id is w_ideal (cool!)'
ENDIF ELSE BEGIN
g=WHERE(w_id ne w_ideal)
help,g
print,'there is a number of cases in which w_ideal is not w_id',n_elements(g)
print,'w_id is not w_ideal, if the number above is much bigger than 1, you should not continue!'
;stop
ENDELSE


e=WHERE(v1 eq w_id)
IF n_elements(e) EQ a THEN BEGIN
print,'The first column of the CMC_i24_test_polynear.cat file and w_id (first column of the output file) are the same (cool!)'
ENDIF ELSE BEGIN
print,'The first column of the CMC_i24_test_polynear.cat file and w_id (first column of the output file) are not the same (bug!)'
stop
ENDELSE


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;histograms;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;'info' files;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

;'info_poly_1': file containing the values of w_poly for w_ideal=1

w_ideal1=WHERE (w_ideal EQ 1.)
d=n_elements(w_ideal1)
print,'the number of galaxies for which w_ideal=1 is',d
poly1=w_poly[w_ideal1]
W1=[transpose(Id[w_ideal1]),transpose(w_ideal[w_ideal1]),transpose(poly1)]
help,W1
help,poly1
openw,10,'./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/info_poly_1'
printf,10,W1
close,10

Nmax_1=max(poly1)
print,'The max of w_poly[w_ideal=1] is',Nmax_1
FOR i=0L,d-1 DO BEGIN
IF poly1[i] EQ max(poly1) THEN print,'the rank of the max is',i
ENDFOR

Nmin_1=min(poly1)
print,'The min of w_poly[w_ideal=1] is',Nmin_1
FOR i=0L,d-1 DO BEGIN
IF poly1[i] EQ min(poly1) THEN print,'the rank of the min is',i
ENDFOR


;'info_poly_0': file containing the values of w_poly for w_ideal=0

w_ideal0=WHERE (w_ideal EQ 0.)
poly0=w_poly[w_ideal0]
nan_arr=finite(poly0)
good=where(nan_arr eq 1)
e=n_elements(w_ideal0)
print,'the number of galaxies for which w_ideal=0 is',e
W3=[transpose(Id[w_ideal0]),transpose(w_ideal[w_ideal0]),transpose(poly0)]
help,W3
openw,30,'./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/info_poly_0'
printf,30,W3
close,30

Nmax_0=max(poly0)
Nmin_0=min(poly0)
print,'The max of w_poly[w_ideal=0] is',Nmax_0
print,'the min of w_poly[w_ideal=0] is',Nmin_0
FOR i=0L,e-1 DO BEGIN
IF poly0[i] EQ min(poly0) THEN print,'the rank of the min is',i
ENDFOR

print,'for w_ideal=1, the values of w_poly go from '+strcompress(strtrim(Nmin_1,2))+' to '+strcompress(strtrim(Nmax_1,2))
print,'for w_ideal=0, the values of w_poly go from '+strcompress(strtrim(Nmin_0,2))+' to '+strcompress(strtrim(Nmax_0,2))

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;histo files;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

Nmax=max(Nmax_1,Nmax_0)
IF(min(poly0) LE Nmin_1)THEN Nmin=min(poly0) ELSE Nmin=Nmin_1
help,Nmin
print,'Thus, the maximum value of  w_poly will be '+strcompress(strtrim(Nmax,2))+' and the minimum will be '+strcompress(strtrim(Nmin,2))

;'histo_poly_1':file containig, for each intervalle of w_poly, the edges of the intervall and the number of galaxies within this intervalle


c=(Nmax-Nmin)/p
print,'the path is',c

wx=FLTARR(p+1)
n1=FLTARR(p+1)
n2=FLTARR(p+1)
n3=FLTARR(p+1)
n4=FLTARR(p+1)

FOR i=0L,p DO BEGIN
        wx1=WHERE (poly1 GE Nmin+i*c AND poly1 LT Nmin+c*(i+1))
	wx0=WHERE (poly0 GE Nmin+i*c AND poly0 LT Nmin+c*(i+1))
        n1[i]=Nmin+i*c
	n3[i]=Nmin+(i+1)*c
	IF (n_elements(wx1) EQ 1) THEN BEGIN
	IF (wx1 EQ -1.) THEN n2[i]=0 ELSE n2[i]=1
        ENDIF ELSE n2[i]=n_elements(wx1)
	IF (n_elements(wx0) EQ 1) THEN BEGIN
	IF (wx0 EQ -1.) THEN n4[i]=0 ELSE n4[i]=1
        ENDIF ELSE n4[i]=n_elements(wx0)
ENDFOR

W2=[transpose(n1),transpose(n3),transpose(n2)]
help,W2
openw,20,'./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/histo_poly1'
printf,20,W2
close,20
help,W2

W4=[transpose(n1),transpose(n3),transpose(n4)]
help, W4
openw,40,'./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/histo_poly0'
printf,40,W4
close,40
help,W4


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;'.ps' files;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

poly0=poly0(good)

;'poly_wideal_1.ps': file containing the histogram of w_poly for w_ideal=1
set_plot,'ps'
device,color=16
device,filename='./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/poly_wideal_1.ps'
plothist, poly1,bin=c,/halfbin,title='histogram of w_poly for w_ideal=1',ytitle='number of galaxies',xtitle='values of w_poly',color=0
device,/close

device,filename='./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/poly_wideal.ps'
plothist,poly1,bin=c,/halfbin,title='combined histogram of w_poly for w_ideal=1 and w_ideal=0',ytitle='number of galaxies',xtitle='values of w_poly',color=0,xrange=[-1,3],yrange=[0,3.5e4],xstyle=1,ystyle=1
plothist,poly1,bin=c,/halfbin,/overplot,color=230
x0=[-0.2] & y0=[5e3]
oplot,x0,y0,color=255
xyouts,-0.2,5e3,'w_ideal=0',color=60
x1=[1.2] & y1=[2.5e4]
oplot,x1,y1,color=255
xyouts,1.2,2.5e4,'w_ideal=1',color=230

;'poly_wideal_0.ps': file containing the histogram of w_poly for w_ideal=0

plothist,poly0,bin=c,/halfbin,/overplot,color=60
device,/close

set_plot,'ps'
device,filename='./histo_results2011/results_NNB'+strcompress(strtrim(NNB))+'/poly_wideal_0.ps'
plothist, poly0,bin=c,/halfbin,title='histogram of w_poly for w_ideal=0',ytitle='number of galaxies',xtitle='values of w_poly'
;plot,n4,n5,psym=2,symsize=3,xtitle='poly',ytitle='number of galaxies',title='histogram of poly values for w_ideal=0'
device,/close

set_plot,'x'
window,0
wset,0
plothist, poly1,bin=c,/halfbin,title='combined histogram of w_poly for w_ideal=1 and w_ideal=0',ytitle='number of galaxies',xtitle='values of w_poly',color=255,xrange=[-1,3],yrange=[0,3.5e4],xstyle=1,ystyle=1
plothist,poly1,bin=c,/halfbin,/overplot,color=230
plothist, poly0,bin=c,/halfbin,/overplot,color=60
x0=[-0.2]& y0=[5e3]
oplot,x0,y0,color=255;,/overplot
xyouts,-0.2,5e3,'w_ideal=0',color=60
x1=[1.2] & y1=[2.5e4]
oplot,x1,y1,color=255;,/overplot
xyouts,1.2,2.5e4,'w_ideal=1',color=230

END