function [xint, yint, xover] = interpolateRafosData(x, y, interval, gapsize, method)
%INTERPOLATERAFOSDATA  interpolate the rafos data over a specified gapsize
%
% SYNOPSIS:
%   function [xint,yint,xover]=interpolateRafosData(x,y,interval,gapsize,method)
%
% DESCRIPTION:
%   INTPRAF4 interpolates data using linear or spline interpolation.
%   If the gap over which was interpolated is greater than the gapsize
%   specified the interpolated points over this gap will be set to zero.
%   Only data points greater than zero will be interpolated.
%
% ARGUMENTS:
%     INPUT:    x,y             Data points to be interpolated over x
%               gapsize         Maximum gap over which to interpolate. The associated unit is that of the x-vector
%   I think this comment is incorrect: (Olaf)   gapsize         maximum number of data points over which interpolation shall be executed
%               interval        
%               method          String containing the method to use
%                               for interpolation. If method = 'spline'
%                               the points will be interpolated using
%                               spline interpolation. If method = 'linear'
%                               the points will be interpolated using
%                               linear interpolation.
%    OUTPUT:    yint            Interpolated y-data points
%               xint            Interpolated y-data points.
%               xover           x-values of interpolated data.
%   The function return the interpolated y-values in yint. All data points
%   which have been interpolated are given back by there indicies xint.
%   The interpolation coverd xover
%	This is the call
%   [itoa_date,itoa]=intpraf4(TOA_DATE(dx),toa(dx),interval,gapsize,imethod);
%  interval = time step of data
%
%  17JUL03, HHF:  This is not returning the correct values for 'xover', the flag which
%  tells whether point has been interpolated or not ... I think that it is returning a
%  '1' for interped only once every day (24hrs), no matter if all cycles < 24 hrs have been 
%  interpolated.
%  17 September 2012 HHF: this program is the problem - with infinite
%  division, e.g. 8./24, it returns time steps that do not match the
%  original msgdate.  I have performed '(round(().*100))./100' to clip to teh
%  hunreths digit all incoming time stamps.  For example, the 8-hourly
%  shedule with yields infinite divisiton is no clipped to 0.08, 0.42, 0.75
%  in day fraction, which include the 2-hr init offset fo rthe listening
%  window.  I will put fix at this function call, but not inside this
%  function.

% if input arguments are missing....
if nargin < 5
  method = 'spline';
end
if nargin < 4
  gapsize = 1;
end


if (strcmpi(method, 'none')) % if the interpolation method was none just return the input
    yint = y;
    xint = x;
    xover = [];
    return
end

if (strcmpi(method, 'cubic'))
    method = 'pchip';
end

[x, indx] = sort(x);   % sort the time axis, which seems to consist of at least 2 sebments
y = y(indx); 

gaps = diff(x);
% this finds the gaps larger than gapsize, over which no interpolation shall be perfomed
isgap = find(gaps > gapsize + 1);
gapstime = NaN(length(isgap), 2);
for ix = 1:length(isgap)
    gapstime(ix, :) = [x(isgap(ix)) x(isgap(ix) + 1)];
end

% this finds the gaps <=  gapsize, over which  interpolation shall be perfomed
isint = find(gaps <= gapsize + 1 & gaps > interval / 24);

% 17JULY03, HHF: Below finds times at which interp is preformed, BUT, decimates by 1 (i.e., 24 hours).
% This means that in a 6-hourly float, only one interp per day happens.
%%$for ix=1:length(isint)
%%$   int_index_array = [int_index_array, [x(isint(ix))+1:1:x(isint(ix)+1)-1]]; % these are the times at which interp is performed
%%$end

% 17JUL03, HHF: This is correct, accounts for vacriable cycle length, and will flag interpolated
% points correctly.
int_index_array = [];
for ix = 1:length(isint)
    int_index_array = [ ...
        int_index_array, ...
        [x(isint(ix))+(interval/24):(interval/24):x(isint(ix)+1)-(interval/24)] ...
        ]; % these are the times at which interp is performed
end

xint = [min(x):interval / 24:max(x)];           % interpolate time axis
yint = interp1(x, y, xint, lower(char(method))); % interpolate the data

%generate array with 1 where interpolated, 0 otherwise
xover = ismember(xint, int_index_array); % HHF, 17JUL03: Ahhhh! This is the problem!

if isempty(isgap)
    xint = toColumn(xint);
    yint = toColumn(yint);
    xover = toColumn(xover);
    return
end

% cut out segments which shall not be interpolated:
for ix = 1:length(isgap)
  gap = find((xint > gapstime(ix, 1)) & (xint < gapstime(ix, 2)));
  xint(gap) = ones(size(gap)) * (-1);   %set xint to -1 for gaps
end
% and now cut them out:
jx = find(xint ~= (-1));
xint = toColumn(xint(jx));
yint = toColumn(yint(jx));
xover = toColumn(xover(jx));


    function [columnVector] = toColumn(vector)
        if ~iscolumn(vector)
            columnVector = vector';
        else 
            columnVector = vector;
        end
    end


end