pro fitstds1,blockid,fit2=fit2,noedit=noedit,showraw=showraw, PS=PS,zdb=zdb,$
    noseg=noseg

; This routine fits standard-star data from the block identified by blockid
; to yield extinction coefficients for the block.  By default, only the
; zero points are fit.
; If keyword fit2 is set, then both zero points and 1st-order extinction
; coefficients are fit.  In any case,
; transformation and 2nd-order extinction coeffs are taken from the
; global_parms.asc file.
; If keyword noedit is set, the routine does not plot residuals en route,
; nor allow the user to edit possibly bad points.
; Results are stored in appropriate arrays in the 'blocks' database.
; During the fitting process, some automatic data-culling is done,
; and diagnostic plots are produced.
; Fits are attempted for all filters for which the block contains data,
; ie, all for which valid=1.
; If keyword showraw is set, then diagnostic plots show raw differences
; between obsd and std magnitudes for each star, not residuals after fit.
; If keyword zdb is set, its value overrides the ZDBASE environment var
; If keyword ps is set, one ps file for all filters and one jpeg file per 
; filter are created within subdirectories of psdir and gifdir named by date.
; If keyword noseg is set, standard obs come from a single file called tostds.
; Otherwise, from a time-segmented database with a name like tostds123.

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

; Trunk directories where PS and GIF files are created.  Make sure they exist!
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/'
;psdir='/knobos/e/kolinski/kepler_diagnostics/test/'
;gifdir='/knobos/e/kolinski/kepler_diagnostics/test/'

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'
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
colstar=[['u','g'],['g','r'],['r','i'],['i','z'],['i','z'],['g','r'],$
         ['g','r']]           ; star dbase filters needed for colors
sigstar=['sig_u','sig_g','sig_r','sig_i','sig_g','sig_Gred','sig_D51']
crng=[[-.5,3.0],[-.5,1.5],[-.3,.7],[-.3,.5],[-.3,.5],[-.5,1.5],[-.5,1.5]]
f1='(f13.5)'
f2='(f10.5)'
mt=40.           ; threshold indicating bad magnitude values
et=0.31          ; threshold for bad magnitude errors.
qthrsh=5.        ; threshold for bad data point rejection, in sigma
astrolib
!priv=2
ss=''
jd0=2453000.2d0   ; base JD for calculating obs segment name
nday=15L          ; number of days covered by an obs segment

; get existing block info
dbopen,blkdbase
dbext,[blockid],'jdmin,jdmax,iseqmin,iseqmax,valid,a,b,c,k',$
              jdmin,jdmax,iseqmin,iseqmax,valid,a0,b0,c0,k0
sjdmin=string(jdmin,format=f1)
sjdmax=string(jdmax,format=f1)
dbclose,blkdbase

; make the 3-digit extension for std obs segment name
if(keyword_set(noseg)) then begin
  ext3=''
endif else begin
  iseg=long(jdmin-jd0)/nday
  siseg=strtrim(string(iseg),2)
  nbiseg=n_elements(byte(siseg))
  if(nbiseg eq 1) then siseg='00'+siseg
  if(nbiseg eq 2) then siseg='0'+siseg
  ext3=siseg
endelse

;Form date name and create date-named directory for plot repository
IF (psflag) THEN BEGIN
  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
  pspath = psdir+yr+mon+day+'/'
  file0 = FINDFILE(pspath, count=ct)
  IF (ct EQ 0) THEN BEGIN
    str0 = 'mkdir '+pspath
    SPAWN, 'umask a+rwx ; ' + str0 + ' ; umask 22'  ; insure directory is world-writable
  ENDIF
  gifpath = gifdir+yr+mon+day+'/'
  file1 = FINDFILE(gifpath, count=ct)
  IF (ct EQ 0) THEN BEGIN
    str1 = 'mkdir '+gifpath
    SPAWN, 'umask a+rwx ; ' + str1 + ' ; umask 22'  ; insure directory is world-writable
  ENDIF

  psfile = pspath+yr+mon+day+'_fitstds_b'+strtrim(string(blockid),2)+'.ps'
  set_plot, 'ps'
  DEVICE, /helvetica, bits_per_pixel=8, filename=psfile, XSIZE=8.0, $
         YSIZE=9, XOFFSET=.25, YOFFSET=1, /INCHES, /PORTRAIT
  set_plot, 'x'
ENDIF

; get current global parameters B and C, and cosmic colors
rd_global,jdmin(0),gnfilt,gfilt_names,gparms,geparms,chipz,cosmic,xglob,zdb=zdb

; make output arrays
a=fltarr(nfilt)+99.9
b=fltarr(nfilt)+99.9
c=fltarr(nfilt)+99.9
k=fltarr(nfilt)+99.9
qual=fltarr(nfilt)+99.9
chipo=fltarr(3,nfilt)+99.9

; loop over filters
for i=0,nfilt-1 do begin
  print
  if(valid(i) gt 0) then begin
    IF (NOT KEYWORD_SET(noedit)) THEN BEGIN
      repeat begin
        print,'Skip filter '+strtrim(filtnam(i),2)+'? (y/n/end)'
        ss = GET_KBRD(1)
        ss=strmid(strtrim(ss,2),0,1)
      endrep until (ss eq 'y' or ss eq 'n' or ss eq 'e')
      if (ss eq 'y') then goto,filtloop
      if (ss eq 'e') then goto,endloop
      print, 'please wait...'
    ENDIF ELSE print, 'Working on filter '+strtrim(filtnam(i),2)

;Create gif filename
    IF (psflag) THEN giffname=gifpath+yr+mon+day+'_fitstds_b'+ $
              strtrim(string(blockid),2)+'_'+strtrim(filtnam(i),2)+'.gif'
  
; make strings to be used to get values out of stars dbase
    fc=filtstar(i)+','
    ec=sigstar(i)+','
    cc1=colstar(0,i)+','
    cc2=colstar(1,i)
    
; go to images dbase, seek images in this filter w/ correct
; time range, go to std obs database, find all obs corresp to these images   
    dbopen,imdbase
    st=dbfind(sjdmin+' < jd < '+sjdmax+',filter='+filtnam(i), /silent)
    dbext,st,'entry,ra,dec,jd,cldflg,astflg',imiseq,imra,imdec,imjd,cldf,astf  
           ; for tagging with reference field name
    dbclose,imdbase

; go to standard obs dbase for this color, find all star obs that come
; from any of these images.
    obsdbase=obspath+filtdir(i)+'/tostds'+ext3
; check that database exists.  If not, set valid entry to zero, bail
    ierrf=findfile(obsdbase+'.dbd')
    berrf=byte(ierrf(0))
    if(berrf(0) eq 0) then begin
      valid(i)=0
      print,'No standard star data for filter '+strtrim(filtnam(i),2)+' !!!'
      goto,filtloop
    endif
    dbopen,obsdbase
    cx=0                             ; to prevent mysterious failures
    sto=dbget('iseq',st,-1,count=cx, /silent)
    nobs=n_elements(sto)
    if(sto(0) gt 0 and nobs ge 3) then begin
      dbext,sto,'starid,iseq,ra,dec,mag,err,x,chip,xcen,ycen,sky',$
        sido,obsseq,ra,dec,mag,err,x,chip,xcen,ycen,sky
      dbclose,obsdbase
      sra=string(ra,format=f2)
      sdec=string(dec,format=f2)
      ch2=fltarr(nobs)      ; arrays =1 for chip=2,3,4 resp
      ch3=fltarr(nobs)
      ch4=fltarr(nobs)
      s2=where(chip eq 2,ns2)
      s3=where(chip eq 3,ns3)
      s4=where(chip eq 4,ns4)
      if(ns2 gt 0) then ch2(s2)=1
      if(ns3 gt 0) then ch3(s3)=1
      if(ns4 gt 0) then ch4(s4)=1

; make arrays for cloud and astrometry flags that go with each obs
      cldfo=intarr(nobs)
      astfo=intarr(nobs)

; go to stars dbase and get standard mags, colors, errors
      stdmag=fltarr(nobs)
      stdcol=fltarr(nobs)
      wts=fltarr(nobs)
      mc1a=fltarr(nobs)
      mc2a=fltarr(nobs)
      serra=fltarr(nobs)
      dbopen,stardbase
      raerr=fltarr(nobs)           ; diagnostic
      decerr=fltarr(nobs)          ; diagnostic
      time=dblarr(nobs)            ; local time in fractional days
      for j=0L,nobs-1 do begin
;       so=dbfind('ra='+sra(j)+'(.00018),dec='+sdec(j)+'(.0024)', /silent)   
        so=dbfind('ra='+sra(j)+'(.00006),dec='+sdec(j)+'(.0008)', /silent)   
        nso=n_elements(so)
        if(so(0) le 0 or nso ne 1) then begin
          so=dbfind('starid='+sido(j), /silent) ; do only if ra search fails
          nso=n_elements(so)
          if(so(0) le 0 or nso ne 1) then begin
          print,'Big Problems in stars database!  Obs for which are no stars!'
          stop
          endif
        endif
        dbext,so,fc+ec+cc1+cc2+',ra,dec',smag,serr,mc1,mc2,starra,stardec
        stdmag(j)=smag
        stdcol(j)=mc1-mc2
        mc1a(j)=mc1
        mc2a(j)=mc2
        serra(j)=serr
        wts(j)=1./(serr^2+err(j)^2)
        raerr(j)=ra(j)-starra(0)
        decerr(j)=dec(j)-stardec(0)
        stim=where(imiseq eq obsseq(j))
        time(j)=imjd(stim(0)) 
        sj=where(imiseq eq obsseq(j),nsj)
        cldfo(j)=cldf(sj(0))
        astfo(j)=astf(sj(0))
      endfor
      minjd=min(long(time))
      time=time-minjd-7./24.          ; constant gets us onto local time, AZ
      time=time*24.
      
; now set weights to zero for invalid mags, errors, cloud or astrometry flags
      s=where(stdmag gt mt or mc1a gt mt or mc2a gt mt or serra gt et $
          or stdcol lt crng(0,i) or stdcol gt crng(1,i) or cldfo gt 0 or $
          astfo gt 0,ns)
      if(ns gt 0) then wts(s)=0.
      sa=where(wts gt 0.,nsa)

; someday set weights to zero for bad columns, etc.

; do a fit, with fitting functions depending on fit2 keyword.
      sgf=where(gfilt_names eq filtnam(i))
      btmp=gparms(0,sgf) & btmp=btmp(0)
      ctmp=gparms(1,sgf) & ctmp=ctmp(0)
      costmp=cosmic(sgf) & costmp=costmp(0)
      ch2tmp=chipz(0,sgf) & ch2tmp=ch2tmp(0)
      ch3tmp=chipz(1,sgf) & ch3tmp=ch3tmp(0)
      ch4tmp=chipz(2,sgf) & ch4tmp=ch4tmp(0)
      aglob=geparms(0,sgf) & aglob=aglob(0)
      kglob=geparms(1,sgf) & kglob=kglob(0)

; test for all zeros in chip-related functions
      if(total(ch2*wts) eq 0. or total(ch3*wts) eq 0. or total(ch4*wts) eq 0)$
             then begin
        print,'No valid data for some chip, filter = '+filtnam(i)
;stop
;       if(keyword_set(global)) then goto,filtloop
      endif

      if(not keyword_set(fit2)) then begin
	nfuns=1
        funs=fltarr(nobs,nfuns)
        dif=mag-stdmag-(x-xglob)*kglob
        funs(*,0)=1.
      endif else begin
        nfuns=2
        funs=fltarr(nobs,nfuns)
        dif=mag-stdmag - stdcol*btmp - (stdcol-costmp)*(x-xglob)*ctmp - $
           ch2tmp*ch2 - ch3tmp*ch3 - ch4tmp*ch4
        funs(*,0)=1.
        funs(*,1)=x-xglob
      endelse
      if(nsa gt nfuns) then begin
	if(nfuns eq 1) then begin
	  cc=total(dif*wts*funs)/total(wts*funs)
	  resid=dif-cc*funs
	endif else begin
          cc=lstsqr(dif,funs,wts,nfuns,rms,chi2,resid,1)
	endelse
      endif else begin
        print,'Too few valid observations, filter ='+filtnam(i)
        goto,filtloop
      endelse

; identify discrepant points, set their weights to zero, redo fit
      igood=fltarr(nobs)+1.
      quartile,resid,med,q,dq
      s=where(abs(resid)/dq ge qthrsh/1.39 or cldfo gt 0 or astfo gt 0,ns)  
;                 factor of 1.39 => dq to sigma
      if(ns gt 0) then igood(s)=0.
fitloop:
      if(nfuns eq 1) then begin
        cc=total(dif*wts*igood*funs)/total(wts*igood*funs)
        resid=dif-cc*funs
	cc=[cc,kglob]
	chi2=total(resid^2*wts*igood)/(total(igood)-1)
      endif else begin
        cc=lstsqr(dif,funs,wts*igood,nfuns,rms,chi2,resid,1)
      endelse
      line=[(aglob-cc(0))-(kglob-cc(1))*xglob,kglob-cc(1)]

; plot results, solicit user interaction
      if(not keyword_set(noedit)) then begin
        if(keyword_set(showraw)) then rr=dif else rr=resid
        resid_plot2,x,stdcol,time,rr,wts,igood,filtstar(i),blockid,psflag, $
          giffname,/edit,crn=crng(*,i),line=line
        print,'Satisfied? (y/n)'
        ss='' &  ss = GET_KBRD(1)
        if(strmid(strtrim(ss,2),0,1) ne 'y') then begin
          psflag=0    ; set to zero so no double plots are produced
          goto,fitloop
        endif
        psflag=keyword_set(ps)  ; set back to original value
      endif ELSE BEGIN
        if(keyword_set(showraw)) then rr=dif else rr=resid
        resid_plot2,x,stdcol,time,rr,wts,igood,filtstar(i),blockid,psflag, $
          giffname,crn=crng(*,i),line=line
        psflag=keyword_set(ps)  ; set back to original value
      ENDELSE
; save results of fit, do next color
      a(i)=cc(0)
      k(i)=cc(1)
      qual(i)=chi2
      if(nfuns ge 6) then begin
        c(i)=cc(2)
        chipo(0,i)=cc(3)
        chipo(1,i)=cc(4)
        chipo(2,i)=cc(5)
        if(nfuns eq 7) then b(i)=cc(6) else b(i)=0.
      endif else begin
        b(i)=btmp
        c(i)=ctmp
        chipo(0,i)=ch2tmp
        chipo(1,i)=ch3tmp
        chipo(2,i)=ch4tmp
      endelse

    endif else begin
      print,'Only ',nobs,' found for filter '+filtnam(i)+' -- no fit'
      stop
    endelse
    
; overwrite original coeff values with new ones
    a0[i]=a[i]
    b0[i]=b[i]
    c0[i]=c[i]
    k0[i]=k[i]
    
  endif
  filtloop:
  
endfor  ;loop over filters
endloop:

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

; update blocks dbase
chip2=reform(chipo(0,*))
chip3=reform(chipo(1,*))
chip4=reform(chipo(2,*))
dbopen,blkdbase,1
dbupdate,[blockid],'a,b,c,k,chip2,chip3,chip4,qual,valid',a0,b0,c0,k0, $
               chip2,chip3,chip4,qual,valid
dbclose,blkdbase

end