%
% Missing Data Interpolation Function
%
%     [ yn, tn ] = missing( y, t, method, par )
%
% where  yn     = new variable time series
%        tn     = new series of times
%        y      = old variable time series
%        t      = old series of times (with missing
%                 points or uneven spacing)
%        method = 1  Fill missing points using average
%               = 2  Fill missing points using simple line
%               = 3  Fill missing points using polynomial
%               = 4  Fill missing points using spline
%               = 5  fixes uneven spacing using spline
%        par    = parameters for the methods
%                 methods 1 & 2 = [] (none required)
%				  method 3 = [ pts_before pts_after order ]
%				  method 4 = [ pts_before pts_after ]
%				  method 5 = [ pts_for_spline pts_overlap ]
%
% Bob Newell, February 1996
%
%------------------------------------------------------
%
function [ yn, tn ] = missing( y, t, method, par )
%
y = y(:);
t = t(:);
n = length( y );
%------------------------------------------------------
% first figure out the missing points
if method < 5,
  disp( 'Searching ...' )
  dt = min( diff( t ) );
  delt = round( diff(t) / dt );
  reg = [];
  tn = t(1);
  yn = y(1);
  count = 0;
  for i = 1:n-1,
    if delt(i) == 1,
	  tn = [ tn; t(i+1) ];
	  yn = [ yn; y(i+1) ];
	elseif delt(i) > 1,
	  tn = [ tn; dt*( 1:delt(i)-1 )' + ...
	              tn( length( tn ) ); t(i+1) ];
	  yn = [ yn; zeros( delt(i)-1, 1 ); y(i+1) ];
	  reg = [ reg; [ count+i+1 count+i+delt(i)-1 ] ];
	  count = count + delt(i) - 1 ;
    end;
  end;
  [ nr, i ] = size( reg );
  if nr > 0, 
    disp( [ 'Found ' num2str( count ) ...
	   ' missing points in ' num2str( nr ) ' regions:' ] )
    disp( reg )
  else,
    disp( 'No missing points found' )
	return;
  end;
end;
%------------------------------------------------------
% replace missing points by the mean value
%
if method == 1,
  newvalue = mean( y );
  for i = 1:nr,
    yn( reg(i,1) : reg(i,2) ) = ...
	    newvalue * ones( reg(i,2)-reg(i,1)+1, 1 );
  end;
%-------------------------------------------------------
% simple linear interpolation
elseif method == 2,
  for i = 1:nr,
    yn( reg(i,1) : reg(i,2) ) = interp1( ...
	         [ tn( reg(i,1)-1 ) tn( reg(i,2)+1 ) ], ...
	         [ yn( reg(i,1)-1 ) yn( reg(i,2)+1 ) ], ...
	         tn( reg(i,1) : reg(i,2) ) );
  end;
%------------------------------------------------------
% polynomial interpolation
elseif method == 3,
  for i = 1:nr,
    ndeg = par(3);
    pc = polyfit( ...
	         [ tn( reg(i,1)-par(1) ) ...
			        tn( reg(i,2)+par(2) ) ], ...
	         [ yn( reg(i,1)-par(1) ) ...
			        yn( reg(i,2)+par(2) ) ], ...
             ndeg );
    yn( reg(i,1) : reg(i,2) ) = polyval( pc, ...
	         tn( reg(i,1) : reg(i,2) ) );
  end;
%------------------------------------------------------
% spline interpolation
elseif method == 4,
  for i = 1:nr,
    yn( reg(i,1) : reg(i,2) ) = spline( ...
	         [ tn( reg(i,1)-par(1) ) ...
			         tn( reg(i,2)+par(2) ) ], ...
	         [ yn( reg(i,1)-par(1) ) ...
			         yn( reg(i,2)+par(2) ) ], ...
	         tn( reg(i,1) : reg(i,2) ) );
  end;
%------------------------------------------------------
% spline correction of uneven data
elseif method == 5,
  dt = ( t(n) - t(1) ) / ( n - 1 );
  tn = t(1) + dt * ( 0:n-1 );
  tn = tn(:);
  yn = zeros( size( tn ) );
  ns = par(1);
  if par(2) > ns/3,
    no = fix( ns / 3 );
  else,
    no = par(2);
  end;
  yn(1:no) = spline( t(1:ns-no), y(1:ns-no), tn(1:no) );
  i = 1;
  while i+ns <= n,
    yn(i+no:i+ns-no-1) = spline( t(i:i+ns-1), ...
	                 y(i:i+ns-1), tn(i+no:i+ns-no-1) );
	i = i + ns - 2 * no;
  end;
  yn(i:n) = spline( t(i-no:n), y(i-no:n), tn(i:n) );    
%------------------------------------------------------
end;
return;
%------------------------------------------------------
% the end
