;+
; NAME:
;     BGND_INTERP
; PURPOSE:
;     Given from one to three spectra representing the background at
;     different intervals, return a flat, linear, or quadratic
;     dynamic spectrum interpolated to the appropriate times.
; CATEGORY:
;     OVRO APC ANALYSIS
; CALLING SEQUENCE:
;     bgdynspec = bgnd_interp(bgspec,tspec,times)
; INPUTS:
;     bgspec    The spectrum or spectra representing the background in
;                 one to three intervals, of size [NF,N], where
;                   N=1 => a flat background,
;                   N=2 => a linear background,
;                   N=3 => a quadratic background,
;                   N>3 => a cubic spline background.
;                 If N=1, a 1-D array of size [NF] may be used as input.
;                 Any missing or bad data should be flagged with NaN.
;     tspec     The center times of the intervals used to generate the
;                 spectra in BGSPEC, in msec.
;     times     An array of length [NTIMES], containing the times [msec]
;                 at which to interpolate/extrapolate the background
;                 spectrum.
; OPTIONAL (KEYWORD) INPUT PARAMETERS:
; ROUTINES CALLED:
; OUTPUTS:
;     bgdynspec The background spectrum, of size [NF,NTIMES], representing
;                 the time-dependent background sampled at times given in
;                 the TIMES array.  If any frequency of the input DATA had
;                 all NaN, then NaN will be returned in that element.
;                 If input data are not as expected, -1 is returned.
; COMMENTS:
; SIDE EFFECTS:
; RESTRICTIONS:
; MODIFICATION HISTORY:
;     Written 23-Nov-2001 by Dale E. Gary
;-
function bgnd_interp,bgspec,tspec,times,accept_ng=accept_ng

   ; Squash any unit dimensions
   datin = reform(bgspec)

   ; Take care of the case where there are still too many dimensions
   if (size(datin,/n_dim) gt 2) then begin
      print,'BGND_INTERP: Error-unexpected number of dimensions in input array.'
      return,-1
   endif

   ; Determine length of TIMES array
   ntimes = n_elements(times)
   nspec = n_elements(tspec)

   ; Take care of the single background spectrum case immediately
   if (size(datin,/n_dim) eq 1) then begin
      if (n_elements(tspec) gt 1) then begin
         print,'BGND_INTERP: Error-unexpected number of dimensions in input array.'
         return,-1
      endif else begin
         ; Make an array of length NTIMES full of 1's
         ones = replicate(1.,ntimes)
         bgdynspec = datin#ones    ; This just makes NTIMES copies of datin
         return,bgdynspec
      endelse
   endif

   ; If there is more than a single background spectrum, verify that the
   ; input data agree.
   if (size(datin,/n_dim) eq 2 and $
      n_elements(datin[0,*]) ne nspec) then begin
      print,'BGND_INTERP: Error-unexpected number of dimensions in input array.'
      return,-1
   endif

   ; Determine the number of frequencies and create output array.
   nf = n_elements(datin[*,0])
   bgdynspec = fltarr(nf,ntimes)

   ; Loop over frequencies
   for i = 0, nf-1 do begin
      ; For this frequency, ensure that there are at least some good data
      good = where(finite(datin[i,*]),ngood)
      if (ngood eq 0) then begin
         bgdynspec[i,*] = !values.f_nan   ; Fill in entire frequency with NaN
         goto,skip
      endif
      ; There are at least some good data for this frequency
      if (ngood lt nspec) then begin
         ; Some of the values for this frequency contain bad data
         if (not keyword_set(accept_ng)) then begin
            ; If the user did not set the keyword to accept bad data,
            ; we have to wipe out the entire frequency
            bgdynspec[i,*] = !values.f_nan   ; Fill in entire frequency with NaN
            goto, skip
         endif
         if (keyword_set(accept_ng) and ngood eq 1) then begin
            ; Take care of case of only one good data value
            bgdynspec[i,*] = datin[i,good]
            goto,skip
         endif
      endif
      ; Now key on number of good frequencies to determine the type of
      ; interpolation to do:
      ;    NGOOD = 2 => Linear fit
      ;    NGOOD = 3 => Quadratic fit
      ;    NGOOD > 3 => Cubic spline
      X = double(tspec[good])
      Y = double(datin[i,good])
      sortx = sort(x)
      x = x(sortx)
      y = y(sortx)
      case NGOOD of
         2: begin
               out = poly_fit(x,y,1)
               bgdynspec[i,*] = out[0] + out[1]*double(times)
            end
         3: begin
               out = poly_fit(x,y,2)
               bgdynspec[i,*] = out[0] + out[1]*double(times) + out[2]*(double(times)^2)
            end
         ELSE: begin
               bgdynspec[i,*] = spline(x,y,double(times),1)
            end
      endcase

skip:
   endfor

   ; All done, so return

return,float(bgdynspec)
end