function yf = fillgaps(y, med_hwidth, efoldt)
% yf = fillgaps(y, med_hwidth, efoldt)





% Turn a single row into a column; otherwise work on each
% column one at a time.
[nr, nc] = size(y);
if nr == 1,
   y = y(:);
   np = nc;
else
   np = nr;
end

yf = y;

n2 = round(med_hwidth/2);      % half-width of median window

for jj = 1:nc
   % Find the gaps:
   gap_mask = isnan(y(:,jj));
   good_mask = ~gap_mask;

   sgm = diff([0 ; gap_mask]) == 1;   % start_gap_mask
   egm = diff([gap_mask ; 0]) == -1;   % end_gap_mask
   istart = find(sgm);
   iend   = find(egm);
   ngaps = length(istart);


   gap_length = iend - istart + 1;

   for ii = 1:ngaps,
      dt = gap_length(ii);
      ig0 = istart(ii);
      ig1 = iend(ii);
      t = (1:dt)' - 0.5;

      % Points to the left of the gap:
      igood_l = find( good_mask( 1:(ig0-1) ) );
      ngl = length(igood_l);
      il0 = max( [ 1 (ngl-n2+1) ] );
      igood_l = igood_l(il0:ngl);  % this should be n2 good points.

      % Points to the right:
      igood_r = ig1 + find( good_mask( (ig1+1):np ) );
      ngr = length(igood_r);
      ir1 = min( [ ngr n2 ] );
      igood_r = igood_r(1:ir1);    % n2 good pts to the right.
      igood_lr = [igood_l ; igood_r];
      pc = polyfit(igood_lr, y(igood_lr,jj), 1);
      yp1 = polyval(pc, ig0-1);
      yp2 = polyval(pc, ig1+1);

      if ig0 == 1,
         y1 = 0;
      else
         y1 = y(ig0-1,jj) - yp1;
      end
      if ig1 == np,
         y2 = 0;
      else
         y2 = y(ig1+1,jj) - yp2;
      end

      q = exp(-dt/efoldt);

      b = (y1 - q * y2) / (1 - q^2);
      c = (y2 - q * y1) / (1 - q^2);

      yf(ig0:ig1,jj) = polyval(pc, (ig0:ig1)') + ...
                         b * exp(-t/efoldt) + c * exp( (t-dt)/efoldt ) ;

   end    % of loop through gaps in a given column

end  % of loop through columns

% If it was a row vector, turn the column back into a row.
if nr == 1, yf = yf.'; end


