pro kic_qc2,density,gmr,imz,jmk,reso=reso,ps=ps,restart=restart
; This routine makes tv plots of star density and mean (weighted by star)
; colors for the region of the Kepler field,with resolution that by default
; is 0.2 degree, but that can be modified with the res keyword.
; stars are counted in bins with constant width res in ra and dec,
; but the density assigned is scaled as 1/cos(dec).
; Density is given in stars per square degree.
; Only stars with r mag le 17 and with valid g,r,i,z,J are counted

; constants
raran=[277.,302.]
decran=[32.,52.]
radian=180./!pi
respath='~/kdbase/kqc2_restrt.sav'

if(not keyword_set(restart)) then begin

; get the data
zdbase=getenv('ZDBASE')
tilelist=zdbase+'/survey/tilenames_all'
getallphot,tilelist,starid,ra,dec,u,g,r,ii,z,gred,d51,j2m,h2m,k2m,mags,names
ra=ra*15.
print,'tp1'

; select out the good stars
sg=where(ra gt 270. and g lt 40 and r lt 17 and ii lt 40 and z lt 40 and $
         j2m lt 40,nsg)
ra=ra(sg)
dec=dec(sg)
g=g(sg)
r=r(sg)
ii=ii(sg)
z=z(sg)
jj=j2m(sg)
k=k2m(sg)
print,'tp2'

; make the output arrays
if(keyword_set(res)) then dr=res else dr=0.2
nra=(raran(1)-raran(0))/dr+1
ndec=(decran(1)-decran(0))/dr+1
density=fltarr(nra,ndec)
gmr=fltarr(nra,ndec)
imz=fltarr(nra,ndec)
jmk=fltarr(nra,ndec)
raout=raran(0)+dr*findgen(nra)
decout=decran(0)+dr*findgen(ndec)

; do the work
print,'nra,ndec=',nra,ndec
for i=0,nra-1 do begin
  if(i mod 10 eq 0) then print,i
  for j=0,ndec-1 do begin
    s=where(ra ge raout(i) and ra lt raout(i)+dr and dec ge decout(j) and $
            dec lt decout(j)+dr,ns)
    if(ns gt 0) then begin
      cdec=cos((decout(j)+dr/2.)/radian)
      density(i,j)=ns/cdec
      gmr(i,j)=mean(g(s))-mean(r(s))
      imz(i,j)=mean(ii(s))-mean(z(s))
      jmk(i,j)=mean(jj(s))-mean(k(s))
    endif
  endfor
endfor

endif else begin
  restore,respath
endelse

save,ra,dec,g,r,ii,z,jj,k,nra,ndec,density,gmr,imz,jmk,raout,decout,file=respath


b1=rotate(bytscl(density,min=0.,max=650.),5) ; flip L/R because of RA
b2=rotate(bytscl(gmr,min=.3,max=1.),5)       ; appearance in sky
b3=rotate(bytscl(imz,min=.0,max=.3),5)
b4=rotate(bytscl(jmk,min=.3,max=.7),5)

npr=8*nra
npd=8*ndec
b1=rebin(b1,npr,npd,/samp)
b2=rebin(b2,npr,npd,/samp)
b3=rebin(b3,npr,npd,/samp)
b4=rebin(b4,npr,npd,/samp)
      
; make plots, with coordinates overlain
sz=size(b1)
xr=[max(raout),min(raout)]
yr=[min(decout),max(decout)]
xtit='RA (degree)'
ytit='Dec (degree)'
tit1='Star Density'
tit2='g-r Color'
tit3='i-z Color'
tit4='J-K Color'
rap=reverse(raout)
decp=reverse(decout)
cc=''

; do plots
if(keyword_set(ps)) then pflg=1 else pflg=0 
  
if(pflg) then begin
  psl,name=ps+'1.ps'
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit1,ticklen=1.,$
    /nodata,xsty=1,ysty=1,$
    charsiz=1.0
  psend
endif else begin
  window,0,xsiz=npr+100,ysiz=npd+100
  erase,0
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit1,ticklen=1.,/noerase,$
    /nodata,xsty=1,ysty=1,/device,charsiz=1.5
  erase,0                    ; get reasonable values into !x.window, etc.
endelse

px=!x.window*!d.x_vsize
py=!y.window*!d.y_vsize
xs=!x.window(1)-!x.window(0)
ys=!y.window(1)-!y.window(0)
p0x=!x.window(0)
p0y=!y.window(0)

if(pflg) then begin
  psl,name=ps+'1.ps'
  tv,b1,p0x,p0y,xsiz=xs,ysiz=ys,/norm
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit1,ticklen=1.,$
  /nodata,xsty=1,ysty=1,/noerase,$
  charsiz=1.0
  psend
endif else begin
  tv,b1,px(0),py(0)
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit1,ticklen=1.,/noerase,$
    /nodata,xsty=1,ysty=1,position=[px(0),py(0),px(0)+sz(1)-1,py(0)+sz(2)-1],/device,$
    charsiz=1.5
endelse

if(not pflg) then read,cc
if(pflg) then begin
  psl,name=ps+'2.ps'
  tv,b2,p0x,p0y,xsiz=xs,ysiz=ys,/norm
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit2,ticklen=1.,$
  /nodata,xsty=1,ysty=1,/noerase,$
  charsiz=1.0
  psend
endif else begin
  erase,0
  tv,b2,px(0),py(0)
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit2,ticklen=1.,/noerase,$
    /nodata,xsty=1,ysty=1,position=[px(0),py(0),px(0)+sz(1)-1,py(0)+sz(2)-1],/device,$
    charsiz=1.5
endelse

if(not pflg) then read,cc
if(pflg) then begin
  psl,name=ps+'3.ps'
  tv,b3,p0x,p0y,xsiz=xs,ysiz=ys,/norm
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit3,ticklen=1.,$
  /nodata,xsty=1,ysty=1,/noerase,$
  charsiz=1.0
  psend
endif else begin
  erase,0
  tv,b3,px(0),py(0)
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit3,ticklen=1.,/noerase,$
    /nodata,xsty=1,ysty=1,position=[px(0),py(0),px(0)+sz(1)-1,py(0)+sz(2)-1],/device,$
    charsiz=1.5
endelse

if(not pflg) then read,cc
if(pflg) then begin
  psl,name=ps+'4.ps'
  tv,b4,p0x,p0y,xsiz=xs,ysiz=ys,/norm
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit4,ticklen=1.,$
  /nodata,xsty=1,ysty=1,/noerase,$
  charsiz=1.0
  psend
endif else begin
  erase,0
  tv,b4,px(0),py(0)
  plot,rap,decp,xran=xr,yran=yr,xtit=xtit,ytit=ytit,tit=tit4,ticklen=1.,/noerase,$
    /nodata,xsty=1,ysty=1,position=[px(0),py(0),px(0)+sz(1)-1,py(0)+sz(2)-1],/device,$
    charsiz=1.5
endelse

if(pflg) then sun
if(not pflg) then read,cc

stop

end