function get_allobs,ra,dec,ids,filtnams,cosmic,$
   xglob,ism,isx,extparms,chpparms,jd,mag,cflg,zdb=zdb
; this routine makes a time series of all observations of a given star, corrected
; for extinction.  On input,
;  ra = RA of desired star, in decimal degrees (used only to get correct tile)
;  dec = Dec of desired star, ditto
;  ids = starid of desired star (used to identify star)
;  filtnams = name of desired filter
;  cosmic,xglob = global cosmic color and zero airmass parms
;  ism(nblk) = min sequence number for each extinction block
;  isx(nblk) = max sequence number for each extinction block
;  extparms(4,7,nblk) = extinction a,b,c,k values for each filter 
;  chpparms(4,7,nblk) = chip parms for each filter (uniformly zero as this is written).
; On return:
;  jd(nobs) = dates of observations
;  mag(nobs) = extinction-corrected magnitudes.
;  cflg(nobs) = cloud flag for obs (0 = okay, not zero = cloudy)
; Returned value is the number of time samples found.

; constants
zdbase=getenv('ZDBASE')
if(keyword_set(zdb)) then zdbase=zdb
path=zdbase+'/survey/'
imdbase=path+'images/images'
blockdbase=path+'/blocks/blocks'
starpath=path+'stars/'
obspath=path+'observations/'
filtdir=['u','g','r','i','z','gr','d5']
colstar=[['u','g'],['g','r'],['r','i'],['i','z'],['i','z'],['g','r'],$
         ['g','r']]
filtnam=['u        ','g        ','r        ','i        ','z        ',$
   'Gred     ','D51      ']
;astrolib
nfilt=n_elements(filtdir)
rad=180.d00/!pi         ; degrees per radian

; get tile required, make db pathnames
tilenames,0,1,fix(15.*ra),fix(dec),stile
tilenames,0,1,fix(15.*ra),fix(dec),otile,/obj
stardb=starpath+stile
for i=0,6 do begin
  if(strtrim(filtnams,2) eq strtrim(filtnam(i),2)) then begin
    ifilt=i
  endif
endfor
obsdbase=obspath+filtdir(ifilt)+'/'+otile(0)

; open stars dbase, find star, make colors.
name=findfile(stardb+'.dbh',count=exists)
if(not exists) then begin
  nobs=0
  jd=[-1.d00]
  mag=[99.99]
  cflg=[1]
  goto,fini
endif

dbopen,stardb
ll=dbcircle(ra,dec,.05,dis,/silent)    ; see if this is fast because of indexed search
dbext,ll,'starid,u,g,r,i,z,gred,d51',sid,uu,gg,rr,ii,zz,gred,d51
s=where(sid eq ids,ns)
if(ns eq 1) then begin
  colors=[uu(s)-gg(s),gg(s)-rr(s),rr(s)-ii(s),ii(s)-zz(s),ii(s)-zz(s),gg(s)-rr(s),$
          gg(s)-rr(s)]
endif else begin
  nobs=0
  jd=[-1.d00]
  mag=[99.99]
  cflg=[1]
  goto,fini           ; no unambiguous match, so bail out.
endelse
dbclose

; look to see if observations dbase exists.  If not, set output parms and bail.
name=findfile(obsdbase+'.dbh',count=exists)
if(not exists) then begin
  nobs=0
  jd=[-1.d00]
  mag=[99.99]
  cflg=[1]
  goto,fini
endif

; open observations dbase, get all obs for star
dbopen,obsdbase
ll=dbcircle(ra,dec,.05,dis,/silent)    ; see if this is fast because of indexed search
dbext,ll,'starid,iseq,mag,x,chip',sido,iseqo,mago,xo,chipo
dbclose
s1=where(sido eq ids,ns1)
if(ns1 gt 0) then begin
  ll=ll(s1)
  sido=sido(s1)
  iseqo=iseqo(s1)
  mago=mago(s1)
  xo=xo(s1)
  chipo=chipo(s1)
  nobs=ns1
endif else begin
  nobs=0
  jd=[-1.d00]
  mag=[99.99]
  cflg=[1]
  goto,fini
endelse

; open images database, get jd values to go with sequence numbers
dbopen,imdbase
dbext,iseqo,'jd,cldflg',jdo,cfgo
dbclose
jd=jdo
cflg=cfgo

; make block number for each observation
blkno=lonarr(nobs)
for i=0,nobs-1 do begin
  sb=where(ism le iseqo(i) and isx ge iseqo(i),nsb)
  if(nsb ge 0) then begin
    blkno(i)=sb(nsb-1)
  endif else begin
    print,'sequence number with no associated block!', iseqo(i)
  endelse
endfor

; make arrays of a,b,c,k,chipcor for each observation
aa=reform(extparms(0,ifilt,blkno))
bb=reform(extparms(1,ifilt,blkno))
cc=reform(extparms(2,ifilt,blkno))
kk=reform(extparms(3,ifilt,blkno))
chipcor=reform(chpparms(chipo,ifilt,blkno))

; make extinction-corrected magnitudes
mag=mago-aa-bb*colors(ifilt)-kk*(xo-xglob)-cc*(colors(ifilt)-cosmic(ifilt))*$
      (xo-xglob)-chipcor

fini:
return,nobs
end

pro pointing_eval1,filin,pointid,npoint,nsamp,rmspnt,rmsall,good,zdb=zdb
; This routine attempts to assess the quality of data for each pointing
; described in the ascii input file filin.
; *** This version of the routine, modified at Dave Latham's behest, tries
;     to compute variability of the r magnitude and of colors g-r, r-i,
;     i-z, and g-d51, not merely the individual magnitudes.
; ***
; For each pointing, it computes the extinction-corrected magnitudes for
; a sample of 20 bright (but not too bright) stars near the center of the
; pointing.  For each filter it then computes the number of distinct
; times (separated by 1 hour or more) that the pointing has been observed,
; and the typical scatter among pointings (if more than 1) for each star.
; Because of long-short exposure pairs, for some filters there are
; more exposures than pointings.  Statistics for the number and scatter
; of all exposures are also returned.
; On return,
;  pointid(npt) = ID number of pointint, encoding RA,Dec eg 39019408
;  npoint(nfilt,npt) = number of distinct visits (> 1 hour separation)
;  nsamp(nfilt,npt) = number of different images at this pointing
;  rmsall(nfilt,npt) = typical scatter among all images taken at this pointing
;  rmspnt(nfilt,npt) = typical scatter when results for each visit are averaged.
;  good(npt) = a quality index for the pointing, defined as follows:
;     -1 = no data for this pointing
;      0 = data for 1 visit but not including all of g,r,i,z,D51 filters
;      1 = all filters present for visit; no 2nd visit in any filter
;      2 = at least 2 visits, but some filters have only 1 visit
;      3 = all filters present for 2 or more visits, but typical scatter
;          values for some or all filters exceed norms.
;      4 = all filters present for 2 or more visits, typical scatter values
;          are within norms

; constants
zdbase=getenv('ZDBASE')
if(keyword_set(zdb)) then zdbase=zdb
path=zdbase+'/survey/'
imdbase=path+'images/images'
blockdbase=path+'/blocks/blocks'
starpath=path+'stars/'
obspath=path+'observations/'
filtdir=['u','g','r','i','z','gr','d5']
colstar=[['u','g'],['g','r'],['r','i'],['i','z'],['i','z'],['g','r'],$
         ['g','r']]
filtnam=['u        ','g        ','r        ','i        ','z        ',$
   'Gred     ','D51      ']
astrolib
nfilt=n_elements(filtdir)
rad=180.d00/!pi         ; degrees per radian
dw=0.30                 ; edge length of box centered on pointings to search
nwant=20                ; want this many stars per pointing for robust estimates
jdnow=2454000.d00       ; JD to get global parms
dt=1./24.               ; images at least 1 hour apart constitute separate visits
rmsthrsh=.03            ; threshold between good/bad pointing repeatability

; get global parameters.  choose a JD that gives parms for KeplerCam, since the cosmic
; colors and xglob never change,
; chip corrections are all zero, and the extinction parms all come out of the
; blocks database anyway. 
rd_global,jdnow,gnfilt,gfilt_names,gparms,geparms,chipz,cosmic,xglob,zdb=zdb

; get extinction parms for all blocks (not just ones in 'blocks' input parm)
dbopen,blockdbase
dbext,-1,'iseqmin,iseqmax',ism,isx
nblka=n_elements(ism)
extparms=fltarr(4,nfilt,nblka)
chpparms=fltarr(4,nfilt,nblka)
for i=0,nblka-1 do begin
  dbext,i+1,'a,b,c,k,chip2,chip3,chip4',a,b,c,k,c2,c3,c4
  extparms(0,*,i)=reform(a,1,nfilt)
  extparms(1,*,i)=reform(b,1,nfilt)
  extparms(2,*,i)=reform(c,1,nfilt)
  extparms(3,*,i)=reform(k,1,nfilt)
  chpparms(1,*,i)=reform( c2,1,nfilt)
  chpparms(2,*,i)=reform( c3,1,nfilt)
  chpparms(3,*,i)=reform( c4,1,nfilt)
endfor

; open input file, read the list of pointings and ra, dec values
openr,iun,filin,/get_lun
ss=''
nline=0
while(not eof(iun)) do begin
  readf,iun,ss
  nline=nline+1
endwhile
point_lun,iun,0
pid=strarr(nline)
rapt=dblarr(nline)
decpt=dblarr(nline)
f1='(i4,1x,a8,6x,a10,3x,a8)'
s1=''
s2=''
s3=''
for i=0,nline-1 do begin
  readf,iun,n1,s1,s2,s3,format=f1
  pid(i)=s1
  hh=double(strmid(s2,0,2))
  hm=double(strmid(s2,3,2))
  hs=double(strmid(s2,6,4))
  dd=double(strmid(s3,0,2))
  dm=double(strmid(s3,3,2))
  ds=double(strmid(s3,6,2))
  rapt(i)=15.*(hh+hm/60.+hs/3600.)   ; RA in decimal degrees
  decpt(i)=dd+dm/60.+ds/3600.        ; Dec in decimal degrees
endfor
close,iun
free_lun,iun

; loop over pointings
npoint=lonarr(nfilt,nline)
nsamp=lonarr(nfilt,nline)
rmsall=fltarr(nfilt,nline)
rmspnt=fltarr(nfilt,nline)
good=intarr(nline)
for ipt=0,nline-1 do begin

; determine tiles that are included in this pointing
  decptp=decpt(ipt)+dw/2.
  idecptp=fix(decptp)
  decptm=decpt(ipt)-dw/2.
  idecptm=fix(decptm)
  raptp=rapt(ipt)+dw/(2.*cos(decpt(ipt)/rad))
  iraptp=fix(raptp)
  raptm=rapt(ipt)-dw/(2.*cos(decpt(ipt)/rad))
  iraptm=fix(raptm)
  tilecr=[iraptm,iraptm,iraptp,iraptp]
  tilecd=[idecptm,idecptp,idecptm,idecptp]
  tilenames,0,4,tilecr,tilecd,tnames
  tilenames,0,4,tilecr,tilecd,tonames,/obj
  so=sort(tnames)
  tnames=tnames(so)
  tonames=tonames(so)
  tu=uniq(tnames)
  tunames=tnames(tu)
  tuonames=tonames(tu)
  ntile=n_elements(tu)

; find all stars in each of the target tiles that are in right brightness
; range.  Concatenate lists for different tiles.
  idlist=['']
  rmag=[0.]
  ras=[0.d00]
  decs=[0.d00]
  for i=0,ntile-1 do begin
    dbname=starpath+tunames(i)

    name=findfile(dbname+'.dbh',count=exists)
    if(not exists) then begin
      print,'notile'
      goto,notile
    endif

    dbopen,dbname
    ll=dbfind('10.5 < r < 15.',/silent)
    if(ll(0) ge 0) then begin
      dbext,ll,'starid,r,ra,dec',sid,r,ra,dec
      idlist=[idlist,sid]
      rmag=[rmag,r]
      ras=[ras,ra]
      decs=[decs,dec]
    endif
    dbclose
    notile:
  endfor

; extract only the stars within ra, dec ranges
  sc=where(15.*ras ge raptm and 15.*ras le raptp and decs ge decptm and $
           decs le decptp,nsc)
  if(nsc gt 0) then begin
    rmag=rmag(sc)
    idlist=idlist(sc)
    ras=ras(sc)
    decs=decs(sc)
  endif else begin
; stop
    goto,nostars
  endelse

; sort this list according to brightness.  Choose top nwant, or fail.  
  nstar=n_elements(rmag)-1

;********testing**********
if(nstar lt nwant) then print,'nostars!  ',rapt(ipt),decpt(ipt)
;**************************

  if(nstar lt nwant) then goto,nostars       ; bail if too few stars
  sob=sort(rmag)
  rmag=rmag(sob(0:nwant-1))
  idlist=idlist(sob(0:nwant-1))
  ras=ras(sob(0:nwant-1))
  decs=decs(sob(0:nwant-1))

; loop over stars
  nobsi=lonarr(nwant,nfilt)
  nvisi=lonarr(nwant,nfilt)
  rmsi=fltarr(nwant,nfilt)
  rmsp=fltarr(nwant,nfilt)

  for istar=0,nwant-1 do begin

; loop over filters, making sequential lists of ifilt, jd, mag
    jdv=[0.d0]
    magv=[0.]
    ifiltv=[0]
    for ifilt=0,nfilt-1 do begin

; get a time series for each star.
      rass=ras(istar)
      decss=decs(istar)
      idlists=idlist(istar)
      filtnams=filtnam(ifilt)
      nobs=get_allobs(rass,decss,idlists,filtnams,cosmic,xglob,$
        ism,isx,extparms,chpparms,jd,mag,cfg)

; reject obs for which cloud flag is set, redefine nobs
      scf=where(cfg eq 0,nobs)
      if(nobs gt 0) then begin
        jd=jd(scf)
        mag=mag(scf)
        cfg=cfg(scf)
      endif

      if(nobs gt 0) then begin
        jdv=[jdv,jd]
        magv=[magv,mag]
        ifiltv=[ifiltv,intarr(nobs)+ifilt]
      endif
      nobsi(istar,ifilt)=nobs
    endfor
    nts=n_elements(jdv)-1
    if(nts gt 0) then begin
      jdv=jdv(1:*)
      magv=magv(1:*)
      ifiltv=ifiltv(1:*)
    endif else begin
      goto,nodat
    endelse
      
;     soj=sort(jd)             ; sort into time order
;     jd=jd(soj)
;     mag=mag(soj)
;     if(nobs eq 0) then goto,endfilt

; first tag each observation with a visit number
    soj=sort(jdv)
    magv=magv(soj)
    ifiltv=ifiltv(soj)
    jdv=jdv(soj)
    nptt=n_elements(jdv)
    ivi=lonarr(nptt)-1

    ii=0
    jd0=jdv(0)
    for iv=0,nptt-1 do begin
      delt=jdv(iv)-jd0
      if(abs(delt) le dt) then begin
        ivi(iv)=ii
      endif else begin
        ii=ii+1
        jd0=jdv(iv)
        ivi(iv)=ii
      endelse
    endfor
    for ifilt=0,nfilt-1 do begin
      sf=where(ifiltv eq ifilt,nsf)
      if(nsf gt 0) then begin
        uvisi=uniq(ivi(sf))
        nvisi(istar,ifilt)=n_elements(uvisi)     ; num of distinct visits per star,filter
      endif else begin
        nvisi(istar,ifilt)=0
      endelse
    endfor
    nvismax=max(ivi)+1          ; num of visits to current star in all filters together

; make averages of magnitudes in separate visits
    magavg=fltarr(nfilt,nvismax)+99.9
    for ifilt=0,nfilt-1 do begin
      for ivis=0,nvismax-1 do begin
        s=where(ivi eq ivis and ifiltv eq ifilt,ns)
        if(ns gt 0) then magavg(ifilt,ivis)=total(magv(s))/ns
      endfor
    endfor

; make photometric indices for each visit.  Indices in order r, g-r, r-i, i-z, g-d51
    pindex=fltarr(nfilt,nvismax)+99.9
    pindex(0,*)=magavg(2,*)                 ; r mag
    difs=[[1,2],[2,3],[3,4],[1,6]]         ; filter indices of desired differences
    for iiv=0,nvismax-1 do begin
      for idd=0,3 do begin
        if(magavg(difs(0,idd),iiv) lt 40 and magavg(difs(1,idd),iiv) lt 40) then begin
          pindex(idd+1,iiv)=magavg(difs(0,idd),iiv)-magavg(difs(1,idd),iiv)
        endif else begin
          pindex(idd+1,iiv)=99.9
        endelse
      endfor
    endfor 

; estimate data dispersion among observations.  Measure used depends on nobs. nobs=0 => 0.
    for ifilt=0,nfilt-1 do begin
      sobs=where(ifiltv eq ifilt and magv lt 40,nobs)
      if(nobs gt 0) then mag=magv(sobs) 
      if(nobs eq 1) then rmsi(istar,ifilt)=0.
      if(nobs eq 2) then rmsi(istar,ifilt)=abs(mag(1)-mag(0))/sqrt(2.)
      if(nobs eq 3) then begin
        som=sort(mag)
        magso=mag(som)
        rmsi(istar,ifilt)=min([magso(1)-magso(0),magso(2)-magso(1)])
      endif
      if(nobs ge 3 and nobs le 6) then rmsi(istar,ifilt)=stddev(mag)
      if(nobs gt 6) then begin
        quartile,mag,med,q,dq
        rmsi(istar,ifilt)=dq/1.349
      endif  
    endfor

; estimate data dispersion among visits, as above
    for ifilt=0,4 do begin
      sii=where(pindex(ifilt,*) lt 40.,nii)
      if(nii gt 0) then magt=reform(pindex(ifilt,sii))
      if(nii eq 1) then rmsp(istar,ifilt)=0.
      if(nii eq 2) then rmsp(istar,ifilt)=abs(magt(1)-magt(0))/sqrt(2.)
      if(nii eq 3) then begin
        sov=sort(magt)
        magsv=magt(sov)
        rmsp(istar,ifilt)=min([magsv(1)-magsv(0),magsv(2)-magsv(1)])
      endif
      if(nii ge 4 and ii le 6) then rmsp(istar,ifilt)=stddev(magt)
      if(nii gt 6) then begin
        quartile,magt,med,q,dq
        rmsp(istar,ifilt)=dq/1.349
      endif
    endfor

    nodat:
  endfor           ; stars loop

;  pointid(npt) = ID number of pointint, encoding RA,Dec eg 39019408
;  npoint(nfilt,npt) = number of distinct visits (> 1 hour separation)
;  nsamp(nfilt,npt) = number of different images at this pointing
;  rmsall(nfilt,npt) = typical scatter among all images taken at this pointing
;  rmspnt(nfilt,npt) = typical scatter when results for each visit are averaged.
;  good(npt) = a quality index for the pointing, defined as follows:
;     -1 = no data for this pointing
;      0 = data for 1 visit but not including all of g,r,i,z,D51 filters
;      1 = all filters present for visit; no 2nd visit in any filter
;      2 = at least 2 visits, but some filters have only 1 visit
;      3 = all filters present for 2 or more visits, but typical scatter
;          values for some or all filters exceed norms.
;      4 = all filters present for 2 or more visits, typical scatter values
;          are within norms


; compute stats for this pointing
  for j=0,nfilt-1 do begin
    npoint(j,ipt)=median(nvisi(*,j))
    nsamp(j,ipt)=median(nobsi(*,j))
    s2=where(rmsi(*,j) gt 0.,ns2)
    if(ns2 ge 3) then rmsall(j,ipt)=median(rmsi(s2,j))
    s3=where(rmsp(*,j) gt 0.,ns3)
    if(ns3 ge 3) then rmspnt(j,ipt)=median(rmsp(s3,j))
  endfor
    
; messy logic to set values in good array
  if(max(npoint(*,ipt)) le 0) then good(ipt)=-1
  sgf=[1,2,3,4,6]               ; corresp [g,r,i,z,d51] filters
  sgf1=[0,1,2,3,4]              ; corresp [r,g-r,r-i,i-z,g-d51]
  if(max(npoint(sgf,ipt)) eq 1 and min(npoint(sgf,ipt)) le 0) then good(ipt)=0
  if(max(npoint(sgf,ipt)) eq 1 and min(npoint(sgf,ipt)) eq 1) then good(ipt)=1
  if(max(npoint(sgf,ipt)) ge 2 and min(npoint(sgf,ipt)) eq 1) then good(ipt)=2
  if(max(npoint(sgf,ipt)) ge 2 and min(npoint(sgf,ipt)) ge 2) then begin
    if(max(rmspnt(sgf1,ipt)) ge rmsthrsh and min(rmspnt(sgf1,ipt)) gt 0.) $
        then good(ipt)=3 else good(ipt)=4
  endif

nostars:
if(max(npoint(*,ipt)) le 0) then good(ipt)=-1
if((ipt mod 10) eq 0) then print,ipt
endfor             ; pointings loop

pointid=pid


end