pro plot_refvsx,blocks,rf,zdb=zdb,ps=ps,log=log

; This routine makes postscript plots of the observed magnitudes of reference
; stars, uncorrected for atmospheric extinction, as a function of airmass.
; It does this for every filter for which observations exist, for every
; block listed in input array blocks.  
; If rf=0, it plots all reference stars that are observed at every opportunity
;          in the given block and filter
; If rf=1, it plots only stars in the R19 field
; If rf=2, it plots only stars in the NGC6811 field.
; If keyword ps is set, ps and gif files are created 
; If keyword log is set, values relating to the fits are written to a log
; file residing in a standard directory

print, '*****current version as of 071210*****'

; Trunk directories where PS and GIF files are created.  Make sure they exist!
; Likewise for logfile directory!
psdir='/mnt/server/tbrown/dbase/diagnostics/plots_ps/'
gifdir='/mnt/server/tbrown/dbase/diagnostics/plots_gif/'
;psdir='/knobos/e/kolinski/kepler_diagnostics/plots_ps/'
;gifdir='/knobos/e/kolinski/kepler_diagnostics/plots_gif/'
logdir='/mnt/server/tbrown/dbase/diagnostics/log/'
psflag = keyword_set(ps)

; constants
zdbase=getenv('ZDBASE')
if(keyword_set(zdb)) then zdbase=zdb
blkdbase=zdbase+'/survey/blocks/blocks'
obspath=zdbase+'/survey/observations/'
stardbase=zdbase+'/survey/stars/tstds'
starpath=zdbase+'/survey/stars'
imdbase=zdbase+'/survey/images/images'
filtdir=['u','g','r','i','z','gr','d5']
nfilt=n_elements(filtdir)
filtnam=['u        ','g        ','r        ','i        ','z        ',$
   'Gred     ','D51      ']
filtstar=['u','g','r','i','z','Gred','D51']   ; filter names in stars dbase
ral=[0.,19.5208,19.6367,8.8561]
decl=[0.,36.1167,46.5667,11.8333]
oname=['/tostds','/tostds','/tostds','/to132+011']
sname=['/tstds','/tstds','/tstds','/t132+011']

nblk=n_elements(blocks)

; loop over blocks, get range of images in block

if(keyword_set(log)) then begin
  openw,iunl,logdir+'refvsxlog_test',/get_lun,/append
  sst=systime()
  printf,iunl,'Created: '+sst
endif

for i=0,nblk-1 do begin

  ib=blocks(i)

print,'block = ',ib

  dbopen,blkdbase
  dbext,[ib],'iseqmin,iseqmax,jdmin',ismin,ismax,jdmin
  dbclose

  ;Form date name and create date-named directory for plot repository
  ;Dave Latham requested plot date change (subtracted 1 from jdmin)
  
  CALDAT, (jdmin-1.0D), mon, day, yr
  yr = STRTRIM(STRING(yr[0]),2)
  mon = STRTRIM(STRING(mon[0]),2)
  IF (STRLEN(mon) LT 2) THEN mon='0'+mon
  day = STRTRIM(STRING(day[0]),2)
  IF (STRLEN(day) LT 2) THEN day='0'+day
  sdatenam=yr+mon+day  ;for log file

  IF (psflag) THEN BEGIN
    pspath = psdir+yr+mon+day+'/'
    file0 = FINDFILE(pspath, count=ct)
    IF (ct EQ 0) THEN BEGIN
      str0 = 'mkdir '+pspath
      SPAWN, str0
    ENDIF
    gifpath = gifdir+yr+mon+day+'/'
    file1 = FINDFILE(gifpath, count=ct)
    IF (ct EQ 0) THEN BEGIN
      str1 = 'mkdir '+gifpath
      SPAWN, str1
    ENDIF

    psfile = pspath+yr+mon+day+'_rvx_b'+strtrim(string(ib),2)+'_test.ps'
    set_plot, 'ps'
    DEVICE, /helvetica, bits_per_pixel=8, filename=psfile, /landscape,  $
      /inches, xsize=10, ysize=8, xoffset=.5, yoffset=10.5
    set_plot, 'x'
  ENDIF

; loop over filters, find observations in tostds in each filter, in
; right range of iseq numbers, in right RA, Dec range if rf ne 0.
  for j=0,nfilt-1 do begin
    
print,'filtno = ',j

    window, 5, retain=2, xsize=800, ysize=650 ;open plot window before loop
    obsdbase=obspath+filtdir(j)+oname(rf)
    ierr=findfile(obsdbase+'.dbd')
    berr=byte(ierr)
    if(berr(0) eq 0) then goto,endloop     ; bail if target dbase doesn't exist
    dbopen,obsdbase
    sisr=strtrim(string(ismin),2)+' < iseq < '+strtrim(string(ismax),2)
    ll=dbfind(sisr)
    if(ll(0) gt 0) then begin
      if(rf gt 0) then begin
        ll1=dbcircle(ral(rf),decl(rf),15.,dis,ll,/silent)
	ll=ll1
      endif
      if(ll(0) lt 0) then goto,nostars
      dbext,ll,'starid,iseq,ra,dec,mag,x,chip',sid,iseq,ra,dec,mag,x,chip
      dbclose

; make list of unique stars, images, get standard magnitudes for all stars
      su=uniq(sid,sort(sid)) & usid=sid(su)
      nsu=n_elements(usid)
      sq=uniq(iseq,sort(iseq)) & uiseq=iseq(sq)
      nui=n_elements(uiseq)
      stardbase=starpath+sname(rf)
      dbopen,stardbase
      ll=dbget('starid',usid)
      dbext,ll,filtstar(j),mag0
      dbclose

; get fwhm and aperture corrections for the images
      dbopen,imdbase
      dbext,uiseq,'apercor,fwhm,ra,jd',apercor,fwhm,raio,jdio
      dbclose

; count observations for each unique star
      nos=lonarr(nsu)
      for kk=0,nsu-1 do begin
        s=where(sid eq usid(kk),ns)
	nos(kk)=ns
      endfor
      nosmax=max(nos)

; pick out subset of stars that were observed on every possible image but one
;     sg=where(nos ge nosmax-1,nsg)
;     if(j eq 3) then sg=where(nos ge nosmax-3,nsg)   ; ***temporary hack***
;      !p.multi = [0,1,2]

; pick out subset of stars that were observed exactly the same number of
; times as many other stars.
      ho=histogram(nos,min=0)
      hom=max(ho,ixh)
      sh=where(ho ge 0.6*hom,nsh)    ; find big peaks in the histogram
      msh=max(sh)                    ; big peak with max no of obs
      ssh=where(sh ge 0.75*msh,nssh) ; big peaks with similar number of obs
      gno=sh(ssh)                   ; no of obs for each of these peaks
      if (nssh gt 4) then gno=gno(nssh-4:nssh-1)  ; take at most 4 biggest
      ngno=n_elements(gno)           ; how many peaks have we got?
      sg=[-1]
      for iz=0,ngno-1 do begin
        ssg=where(nos eq gno(iz),nssg)  ; take stars with exactly this no of obs
	if(nssg gt 0) then sg=[sg,ssg]
      endfor
      sg=sg(1:*)
      nsg=n_elements(sg)

      if(nsg gt 5) then begin

; make medians of (mags of all stars on each image) and (mags-mags0)
; of all stars on each image, likewise of airmass for these stars on each image
; but only for stars with reasonable mag (or mag, mag0 for difference).
        medm=fltarr(nui)
        meddif=fltarr(nui)
        errdif=fltarr(nui)
        medx=fltarr(nui)
        tmag=fltarr(nsg,nui)+99.9
        tmag0=fltarr(nsg,nui)+99.9
        tdif=fltarr(nsg,nui)+99.9
        tx=fltarr(nsg,nui)+99.9
; also make aperture corrections, fwhm weighted by number of stars for ea image
        apca=fltarr(nui)
        fwhma=fltarr(nui)
        for mm=0,nui-1 do begin
          fwa=0.
          apa=0.
  	  ntot=0.
	  for cc=0,3 do begin
            si=where(iseq eq uiseq(mm) and chip eq cc,nsi)
	    apa=apa+nsi*apercor(cc,mm)
	    fwa=fwa+nsi*fwhm(cc,mm)
	    ntot=ntot+nsi
  	  endfor
  	  apca(mm)=apa/ntot
	  fwhma(mm)=fwa/ntot

	  for nn=0,nsg-1 do begin
	    st=where(sid eq usid(sg(nn)) and iseq eq uiseq(mm),nst)
	    if(nst eq 1) then begin
	      tmag(nn,mm)=mag(st)
	      tmag0(nn,mm)=mag0(sg(nn))
	      tdif(nn,mm)=mag(st)-mag0(sg(nn))
	      tx(nn,mm)=x(st)
	    endif
          endfor
	  sug=where(tmag(*,mm) lt 40,nsug)
  	  sue=where(tmag(*,mm) lt 40 and tmag0(*,mm) lt 40,nsue)
	  if(nsue ge 3) then begin
	    medm(mm)=median(tmag(sug,mm))
	    meddif(mm)=median(tmag(sue,mm)-tmag0(sue,mm))
	    quartile,tmag(sue,mm)-tmag0(sue,mm),medofd,qrt,dqrt
            errdif(mm)=dqrt
	    medx(mm)=median(tx(sug,mm))
	  endif
        endfor

; compute linear fits to meddif vs medx, residuals
        if(nui ge 2 and max(abs(medx)) gt 0.) then begin
          ss=where(medx ne 0.,nss)
	  if(nss gt 0) then begin
	    medx=medx(ss)
	    meddif=meddif(ss)
	    medm=medm(ss)
	    nui=nss
            rai=raio[ss]
            jdi=jdio[ss]
	  endif
          dd=poly_fit(medx,meddif,1)
	  azero=dd(0)
	  slope=dd(1)

; look for and reject outliers in the residuals
          yyf=poly(medx,dd)
	  resid=meddif-yyf
	  nresid=n_elements(resid)
	  if(nresid ge 5) then begin
	    quartile,resid,medrr,qrr,dqrr
	    dmed=abs(resid-medrr)/dqrr
	    sgdd=where(dmed le 3.*1.349,nsgdd)     ; data within 3-sigma 
	    if(nsgdd lt nresid and nsgdd ge 4) then begin
	      dd=poly_fit(medx(sgdd),meddif(sgdd),1)
	      azero=dd(0)
	      slope=dd(1)
	      print,nresid-nsgdd,'Bad data rejected, filter=',j
	    endif
	  endif
        endif else begin
          azero=meddif(0)
	  slope=0.
          dd=[azero,slope]
          if(max(abs(meddif)) eq 0.) then begin
            print,'no nonzero meddif -- bailing'
            goto,endloop
          endif
        endelse
        yyf=poly(medx,dd)
        resid=meddif-yyf
	if(n_elements(resid) ge 2) then medrms=stddev(resid) else medrms=99.
	minx=min(medx)                ; smallest image median airmass
	maxx=max(medx)                ; largest
	nxpts=n_elements(medx)        ; number of points to which fit is made
	msiga=medrms/sqrt((nxpts-2) > 1)
	mxran=(maxx-minx) > .01
	msigk=medrms/(mxran*sqrt((nxpts-2) > 1))   ; rough and ready estimates,
	                                           ; do better later
	f1='(i5,2x,1a8,2x,1a4,7(2x,f7.3),i6)'
	if(keyword_set(log)) then begin
	  fs=filtstar(j)
	  printf,iunl,ib,sdatenam,fs,azero,msiga,slope,msigk,medrms,$
	      minx,maxx,nsg,format=f1
	endif

        xx=.8+.5*findgen(4)
        yyp=poly(xx,dd)
        stra=string(azero+1.215*slope,format='(f6.3)')
        strs=string(slope,format='(f6.3)')
	strsiga=string(msiga,format='(f7.4)')
	strsigk=string(msigk,format='(f7.4)')
        strrms=string(medrms,format='(f7.4)')

; get hourangle
        get_hourangle, rai, jdi, hrai
        ha = INTARR(nui)
        sha = WHERE(hrai GT 0, nsha)
        IF (nsha GT 0) THEN ha[sha]=1

; make plots of useful stuff
        IF (psflag) THEN BEGIN ;plot ps and gif files first, then to screen
          set_plot, 'ps'    
          !p.multi = [0,1,2]     
          tit='Test: Block '+strtrim(string(blocks(i)),2)+'  filter '+filtstar(j)
          xtit='Airmass'
          ytit='Magnitudes'
          xr=[0.9,2.2]
          medmed=total(medm)/nui
          difmed=total(meddif)/nui
          yr1=[medmed-0.2,medmed+0.3]
          yr2=[difmed-.3,difmed+.4]
  
          plot,medx,meddif,xran=xr,yran=yr2,tit=tit,xtit=xtit,ytit=ytit, $
             /nodata,/xsty,/ysty
          for mm=0,nui-1 do begin
            IF (ha[mm] EQ 0) THEN symb=4 ELSE symb=5
            oplot,[medx(mm)],[meddif(mm)],psym=symb,thick=3
;           oplot,tx(*,mm),tdif(*,mm),psym=3
            oploterr,[medx(mm)],[meddif(mm)],[errdif(mm)]
          endfor
          oplot,xx,yyp,line=2
          xyouts,1.,yr2(1)-.05, 'Triangle symbol denotes positive hour angle'
          xyouts,1.,yr2(1)-.1,'a = '+stra+' +/- '+strsiga
          xyouts,1.,yr2(1)-.15,'k = '+strs+' +/- '+strsigk
          xyouts,1.,yr2(1)-.2, 'rms = '+strrms

          ytit2='Fit Residual (mag)'
          xtit1='FWHM (pix)'
          xr1=[0.,12.]
          yr1=[-.06,.06]
          plot,fwhma,resid,psym=2,tit=tit,xtit=xtit1,ytit=ytit2,xran=xr1,$
            yran=yr1,/xsty,/ysty  

          set_plot, 'x'  
        ENDIF

        !p.multi = [0,1,2] 
        tit='Test: Block '+strtrim(string(blocks(i)),2)+'  filter '+filtstar(j)
        xtit='Airmass'
        ytit='Magnitudes'
        xr=[0.9,2.2]
        medmed=total(medm)/nui
        difmed=total(meddif)/nui
        yr1=[medmed-0.2,medmed+0.3]
        yr2=[difmed-.3,difmed+.4]
      
        plot,medx,meddif,xran=xr,yran=yr2,tit=tit,xtit=xtit,ytit=ytit,/nodata,$
          /xsty,/ysty
        for mm=0,nui-1 do begin
          IF (ha[mm] EQ 0) THEN symb=4 ELSE symb=5
          oplot,[medx(mm)], [meddif(mm)], psym=symb, thick=3
;          oplot,tx(*,mm),tdif(*,mm),psym=3
          oploterr,[medx(mm)],[meddif(mm)],[errdif(mm)]
        endfor
        oplot,xx,yyp,line=2
        xyouts,1.,yr2(1)-.05, 'Triangle symbol denotes positive hour angle'
        xyouts,1.,yr2(1)-.1,'a = '+stra+' +/- '+strsiga
        xyouts,1.,yr2(1)-.15,'k = '+strs+' +/- '+strsigk
        xyouts,1.,yr2(1)-.2, 'rms = '+strrms

        ytit2='Fit Residual (mag)'
        xtit1='FWHM (pix)'
        xr1=[0.,12.]
        yr1=[-.06,.06]
        plot,fwhma,resid,psym=2,tit=tit,xtit=xtit1,ytit=ytit2,xran=xr1,$
          yran=yr1,/xsty,/ysty

        ;Create gif
        IF (psflag) THEN BEGIN
          giffname=gifpath+yr+mon+day+'_rvx_b'+ $
            strtrim(string(ib),2)+'_'+strtrim(filtnam(j),2)+'_test.gif'
          tvlct, r, g, b, /get
          r[0]=255 & g[0]=255 & b[0]=255
          r[255]=0 & g[255]=0 & b[255]=0
          gif_image = tvrd()
          write_gif, giffname, gif_image, r, g, b
        ENDIF
      endif  ;from: if(nsg gt 5) then begin
      nostars:
      wdelete, 5
    endif
    endloop:
  endfor        ; end of filters loop

  if (psflag) then begin
    set_plot, 'ps'
    device, /close
    set_plot, 'x'
  endif

endfor          ; end of blocks loop
!p.multi=[0,1,1]
if (keyword_set(log)) then begin
  close, iunl
  free_lun,iunl
endif

;stop
end