;+
; NAME:
;     TPINTERP
; PURPOSE:
;     Interpolate missing TP calibration factors, based on adjacent good measurements
;     and the ND increment at the missing frequency, for 27-m total power calibration.
; CATEGORY:
;     OVRO APC CALIBRATION
; CALLING SEQUENCE:
;     tpinterp,fghz,tpfac,tprms,ndavg,tpnew
; INPUTS:
;     fghz    The array of frequencies at which the factors were measured
;     tpfac   The array of total power factors, as determined from TPSOLVE,
;               as an array of size (2,2,85), where the first index is the
;               feed (R-position or L-position), the second is the antenna
;               (antennas 1 and 2), and the third is the frequency.
;     tprms   The corresponding array of RMS variations (not used)
;     ndavg   The noise diode increments corresponding to the input data,
;               of identical size and meaning to TPFAC.
; OPTIONAL (KEYWORD) INPUT PARAMETERS
; ROUTINES CALLED:
; OUTPUTS:
;     tpnew   The output array, same as TPFAC but with missing frequencies
;               interpolated.
; COMMENTS:
; SIDE EFFECTS:
;     A total power calibration file (.TPC file) is written)
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 02-Oct-1998 by Dale E. Gary
;     04-Oct-1998  DG
;       Improve results by normalizing to the noise diode (which should
;       take out the frequency-to-frequency variations due to the receiver,
;       standing waves, etc.), eliminating obviously bad measurements,
;       fitting a smooth gaussian, interpolating missing measurements,
;       then de-normalizing to obtain final factors.
;     22-Jul-1999  DG
;       Changed MAXDIF from 2 to 10, due to the fact that there really is
;       a systematic difference between (R-L) and (L-R) as defined here.
;       Typically it can reach a value of 6 and still be perfectly good
;       data.
;-
pro tpinterp,fghz,tpfac,tprms,ndavg,tpnew

   tpnew = tpfac
   for ifeed = 0, 1 do begin
      for iant = 0, 1 do begin
         ; Remove frequencies with missing ND value from further consideration
         igood = where(ndavg(ifeed,iant,*) gt 0,ngood)
         f = fghz(igood)
         tp = reform(tpfac(ifeed,iant,igood))
         nd = reform(ndavg(ifeed,iant,igood))
         ; Calculate the total power normalized to the ND
         ratio = tp/nd

         ; Find those values with missing total power measurement (indicated by negative RATIO)
         bad = where(ratio lt 0, nbad)
         if (nbad gt 0) then begin
            ; There are some missing TP values
            good = where(ratio gt 0, ngood)
            if (ngood gt 2) then begin
               ; And there are some good values.  Get the relevant values of F and RATIO
               tpgood = tp(good)
               rgood = ratio(good)
               fgood = f(good)
            endif else begin
               print, 'Not enough good data for Feed ',ifeed,' Ant ',iant+1
               goto,skip
            endelse
            ; Loop over the bad TP values, and replace them with estimates based on
            ; the RATIO at nearby good frequencies
            for i = 0, nbad-1 do begin
               ; Find the indices of the nearest two good frequencies.
               if (f[bad[i]] lt fgood[0]) then begin
                  ; We have to extrapolate to lower frequencies
                  k1 = 0
                  k2 = 1
               endif else if (f[bad[i]] gt fgood[ngood-1]) then begin
                  ; We have to extrapolate to higher frequencies
                  k1 = ngood-2
                  k2 = ngood-1
               endif else begin
                  ; Our bad frequency falls somewhere in the middle of the good
                  ; frequencies, so we can interpolate.  Find the nearest neighbors.
                  k1 = (where(f[bad[i]] gt fgood,nlt))(nlt-1)
                  k2 = (where(f[bad[i]] lt fgood,ngt))(0)
               endelse
               tp[bad[i]] = (rgood[k1] + (rgood[k2]-rgood[k1])* $
                            ((f[bad[i]]-fgood[k1])/(fgood[k2]-fgood[k1])))*nd[bad[i]]
            endfor
            ; Restore TP, including corrections, into output array
         endif
skip:
         tpnew[ifeed,iant,igood] = tp
      endfor
   endfor

return
end