pro color_fillin1,colorsin,indexin,colorsout,indexout,nchg,nf=nf
; This routine reads a save file colorsin containing colors and corresponding
; values of teff, logg, and logz, and also a save file indexin containing
; arrays of unique values teffu, loggu, logzu and the index file that points
; to values in the colors table, given indices into teffu, loggu, and logzu.
; This index file must have been written by basel_index1, or by this routine.
; It then systematically steps through the 3D space defined by teffu, loggu,
; logzu, in an attempt to locate cells that do not contain valid colors,
; but that have enough valid neighboring (up to 2 away) cells that a sensible interpolation
; or extrapolation can be performed.  This is defined as at least 40 valid
; neighbors (out of a possible 124).  If enough are found, then for each
; color a 3D quadratic is fit to the valid data:
;   c = A + B*x + C*y + D*z + E*x^2 +F*y^2 + G*z^2 + H*x*y + I*x*z + J*y*z
; where x = alog10(teff)-alog10(teff_0)
;       y = logg - logg_0
;       z = logz - logz_0
; and teff_0, logg_0, logz_0 are the values at the central point.  Invalid
; points are given zero weight, and one round of iteration is done to reduce
; the weights of extreme outliers.
; Interpolated colors are written to a new colors table, and the corresponding
; values of teff, logg, logz are written to their own arrays.  At the
; end of the processing, all of these arrays are concatenated appropriately,
; and a new colors file and index file are written to colorsout, indexout.
; The number of new points supplied by interpolation is returned in nchg.
; If keyword nf is set, the number of functions fitted is set to nf <= 10. (4 is a good value)

; constants
nmin=30                           ; min acceptable no of good neighbors
rrmax=2.e4                        ; max acceptable ratio of singular values
if(keyword_set(nf)) then nf0=nf

; restore the input data
restore,colorsin
restore,indexin

; get sizes of things
ntu=n_elements(teffu)
nlgu=n_elements(loggu)
nlzu=n_elements(logzu)
logtef=alog10(teff)
logtefu=alog10(teffu)
sz=size(colors)
nc=sz(1)
npt=sz(2)
nctot=0

loopc=0
loop:

; make dummy arrays for added colors, teff, logg, logz
cola=fltarr(nc)
teffa=[0.]
logga=[0.]
logza=[0.]
indxo=indx         ; output indx array
nchg=0

; do the work
for i=0,ntu-1 do begin
  for j=0,nlgu-1 do begin
    for k=0,nlzu-1 do begin
      if(indxo(i,j,k) eq -1) then begin   ; if colors ok, leave them alone

; count up valid neighbors
        iran=[(i-2) > 0, (i+2) < (ntu-1)]
        jran=[(j-2) > 0, (j+2) < (nlgu-1)]
        kran=[(k-2) > 0, (k+2) < (nlzu-1)]
        indnei=indxo(iran(0):iran(1),jran(0):jran(1),kran(0):kran(1))
        s=where(indnei gt -1,ns)
        if(ns ge nmin) then begin

; if we get here, indnei(s) contains the indices of the colors, teff, etc
; entries that are needed for the fits.  Do the fits.
          indc=indnei(s)
          tef0=alog10(teffu(i))
          lgg0=loggu(j)
          lgz0=logzu(k)
          xx=alog10(teff(indc))-tef0
          yy=logg(indc)-lgg0
          zz=logz(indc)-lgz0
          cfit=fltarr(nc)

; make array of functions to fit to
          funs=fltarr(ns,10)
          funs(*,0)=1.
          funs(*,1)=xx
          funs(*,2)=yy
          funs(*,3)=zz
          funs(*,4)=xx^2
          funs(*,5)=yy^2
          funs(*,6)=zz^2
          funs(*,7)=xx*yy
          funs(*,8)=xx*zz
          funs(*,9)=yy*zz

; create array of cross-products
          aa=fltarr(10,10)
          for m=0,9 do begin
            for n=0,9 do begin
              aa(m,n)=total(funs(*,m)*funs(*,n))
            endfor
          endfor

; compute SVD to see if we have singularity problems.  If so, reduce
; number of functions.
          svdc,aa,w,u,v
          rr=max(w)/min(w)
          if(rr lt rrmax) then begin
            nfun=10
          endif else begin
            svdc,aa(0:6,0:6),w,u,v
            rr=max(w)/min(w)
            if(rr lt rrmax) then begin
              nfun=7
            endif else begin
              svdc,aa(0:3,0:3),w,u,v
              rr=max(w)/min(w)
              if(rr lt rrmax) then begin
                nfun=4
              endif else begin
                goto,nofit
              endelse
            endelse
          endelse
       
          if(keyword_set(nf)) then nfun=nf0

          ct=fltarr(nc)   
          wts=fltarr(ns)+1.
          for m=0,nc-1 do begin
            dat=reform(colors(m,indc(s)),ns)
            cc=lstsqr(dat,funs,wts,nfun,rms,chisq,outp,1)
            s1=where(abs(outp/rms) ge 4.,ns1)
            if(ns1 gt 0) then begin
              wts(s1)=(4./(abs(outp(s1)/rms)))^2
              cc=lstsqr(dat,funs,wts,nfun,rms,chisq,outp,1)
            endif
            ct(m)=cc(0)
          endfor
        cola=[cola,ct]
        teffa=[teffa,teffu(i)]
        logga=[logga,loggu(j)]
        logza=[logza,logzu(k)]
        indxo(i,j,k)=npt+nchg
        nchg=nchg+1  

        endif
        nofit:
      endif
    endfor
  endfor
endfor

if(nchg gt 0) then begin
  cola=cola(nc:*)
  teffa=teffa(1:*)
  logga=logga(1:*)
  logza=logza(1:*)

  colors=reform([reform(colors,nc*long(npt)),cola],nc,npt+nchg)
  teff=[teff,teffa]
  logg=[logg,logga]
  logz=[logz,logza]
  ext=[ext,fltarr(nchg)]
  rv=[rv,fltarr(nchg)]
  indx=indxo
  npt=npt+nchg
endif

print,nchg
nctot=nctot+nchg
if(nchg gt 0) then begin
  loopc=loopc+1
  goto,loop
endif

nchg=nctot
save,colors,teff,logg,logz,ext,rv,names,file=colorsout
save,indx,teffu,loggu,logzu,file=indexout

end