;
; read Cosmos star list and plot
;
clearwins
dir = 'C:\Documents and Settings\Don Gavel\My Documents\Adaptive Optics\Keck\Keck NGAO\Keck NGAO Extragalactic Science Modeling\'
file = 'usno-search-COSMOS.txt'
u = read_ascii(dir+file,comment=';')
u = u.field01
ra = u[8,*] ; this is actually RA converted to degrees
dec = u[9,*] ; degrees
Bmag = u[18,*] ; B band magnitude
Rmag = u[19,*] ; R band magnitude
Imag = u[20,*] ; I band magnitude
;window,/free,xsize=700,ysize=700
;plot,ra,dec,xrange=minmax(ra),yrange=minmax(dec),psym=3,title='COSMOS field',xtitle='RA',ytitle='Dec'
;
n = 600 ; we'll divide the field into 0.1 arcmin (6 arcsec) patches
ra0 = min(ra)
dra = (max(ra)-min(ra))/float(n)
dec0 = min(dec)
ddec = (max(dec)-min(dec))/float(n)
ira = (ra-ra0)/dra < (n-1)
idec = (dec-dec0)/ddec < (n-1)
a = fltarr(n,n,3)
aa = fltarr(n,n)
bg = 25
aa -= bg
aa[ira,idec] = -(Imag<bg)
a[*,*,0] = aa
aa[ira,idec] = -(Rmag<bg)
a[*,*,1] = aa
aa[ira,idec] = -(Bmag<bg)
a[*,*,2] = aa
disp,a,x0 = [ra0,dec0], dx = [dra,ddec]
b = transpose(bytscl((a+bg)^(.2)),[2,0,1])
window,/free,xsize=n,ysize=n
tv,b,/true ; pretty color picture
;
; statistics: star/sq degree vs magnitude
dmag = .25
minmag = 8
maxmag = 21
nmags = (maxmag-minmag)/dmag + 1
mags = findgen(nmags)*dmag+minmag

maghist = fltarr(nmags,3)
for color = 0,2 do maghist[*,color] = histogram(-a[*,*,color],min=minmag,max=maxmag,nbins=nmags)
band = ['I magnitude','R magnitude','B magnitude']
ytitle = 'stars/sq degree'
window,/free,title='COSMOS field - Star distribution'

!p.multi = [0,2,2,0,0]
mx = reform(transpose([[mags],[mags]]),2*nmags)
mx = shift(mx,-1)
mx[2*nmags-1]=nmags*dmag+minmag
for color = 0,2 do begin
	my = maghist[*,color]
	my = reform(transpose([[my],[my]]),2*nmags)
	plot,mx,my>1,/ylog,xrange=[minmag,maxmag],xstyle=1,xtitle=band[color],ytitle=ytitle,charsize=1.5
endfor
!p.multi = 0
;
; statistics: cumulative star counts vs magnitude
magcum = fltarr(nmags,3)
ytitle = 'stars/sq deg < mag'
for color = 0,2 do magcum[*,color] = total(maghist[*,color],/cumulative)
window,/free,title='COSMOS field - Cumulative star counts'
!p.multi = [0,2,2,0,0]
mx = reform(transpose([[mags],[mags]]),2*nmags)
mx = shift(mx,-1)
mx[2*nmags-1]=nmags*dmag+minmag
for color = 0,2 do begin
	my = magcum[*,color]
	my = reform(transpose([[my],[my]]),2*nmags)
	plot,mx,my>1,/ylog,xrange=[minmag,maxmag],xstyle=1,xtitle=band[color],ytitle=ytitle,charsize=1.5
endfor
!p.multi = 0
;
end