function [time_spl, data_spl, deriv_spl] = splineWithDummy(pTime, pData, pDummy, pInterval)
% function [time_spl,data_spl,deriv_spl] = spline_with_dummy(time,data,dummy,interval);
%
% Inputs:
% time: time basis [usually in RAFOS days]
% data: whatever data in whatever unit
% dummy: the dummy that is used in the input 'data'
% interval: the interval at which the resulting spline shall be subsampled, e.g. every 2h. Given in [days]. 
% if dummy is numerical and >0, e.g. 999, values smaller than dummy are assumed to be valid
% if dummy is numerical and <0, e.g. -1,  values greater than dummy are assumed to be valid
%
% Outputs:
% time_spl: subsampled time, steady without any dummy
% data_spl: splined and subsampled data; with NaNs at intervals where dummys are found in the input. 
% deriv_spl: splined and subsampled first derivative; with NaNs at intervals where dummys are found in the input. 
%            units are [data_spl]/[time_spl]
%
% Olaf Boebel, URI 22 August 2000

% For testing purposes
%[header,data] = loadrfc('C:\data\floats\rafos\final\kapex\','rfc195.rfc');
%time = data(:,3);
%data = data(:,5);
%dummy = 0;

% determine spline interval
timestep = pTime(2)-pTime(1); % this is the time step of the raw data
if (nargin == 4)
   si = pInterval;   % the time step of the splined and interpolated output data
else
   si = timestep; %: output at same intervals as input
end

% define the dummy value
if (nargin < 3)
   pDummy = 'NaN';
end

% find the dummies
if strcmp(pDummy,'NaN')
   index = find(isnan(pData));
elseif pDummy > 0
   index = find(pData >= pDummy);
elseif pDummy <= 0
      index = find(pData <= pDummy);
end
index_2 = find(diff(pTime(index))>timestep);    %find the transition points

if ~isempty(index) 
   if ( pTime(index(length(index))) == pTime(length(pTime)) && pTime(index(1)) == pTime(1) ) % the first and last data points are dummies
      mask_interval = [pTime([index(1); index(index_2+1)-1]) pTime([index((index_2))+1; index(length(index))])] ;  %  -si is needed to compensate the + si below
   elseif ( pTime(index(length(index))) == pTime(length(pTime)) )                         % the  last data point is a dummy
   	mask_interval = [pTime([index(1); index(index_2+1)]-1) pTime([index((index_2))+1; index(length(index))])] ;
   elseif ( pTime(index(1)) == pTime(1))                                               % the first data point is a dummy
      mask_interval = [pTime([index(1); index(index_2+1)-1]) pTime([index((index_2)); index(length(index))]+1)] ;
   else                                                                              % only interior dummies
      mask_interval = [pTime([index(1); index(index_2+1)]-1) pTime([index((index_2)); index(length(index))]+1)] ;
   end
   time_gap_array = [];
   
   % modified by T. reynaud and G. Roudaut 13/07/05
   %function rounddigit added to get ride of roundoff errors while
   %comparing two real numbers
   for indi = 1:size(mask_interval,1)
      if (rounddigit(pTime(index(indi)),5) == rounddigit(pTime(indi),5) && indi == 1) % the first data point is a dummy
            if (rounddigit(mask_interval(indi,1)+si,5) <= rounddigit(mask_interval(indi,2)-si,5)) % at least one element in between
                time_gap = linspace(mask_interval(indi,1),mask_interval(indi,2)-si,-(mask_interval(indi,1)-mask_interval(indi,2))*1/si)';
            elseif (rounddigit(mask_interval(indi,1)+si,5) == rounddigit(mask_interval(indi,2),5)) % no element in between
                time_gap = linspace(mask_interval(indi,1),mask_interval(indi,2),-(mask_interval(indi,1)-mask_interval(indi,2))*1/si+1)';
            end
      else
      	if (rounddigit(mask_interval(indi,1)+si,5) <= rounddigit(mask_interval(indi,2)-si,5)) % at least one element in between
      	   time_gap = linspace(mask_interval(indi,1)+si,mask_interval(indi,2)-si,-(mask_interval(indi,1)-mask_interval(indi,2))*1/si-1)';
      	elseif  (rounddigit(mask_interval(indi,1)+si,5) == rounddigit(mask_interval(indi,2),5)) % no element in between
      	   time_gap = linspace(mask_interval(indi,1)+si,mask_interval(indi,2),-(mask_interval(indi,1)-mask_interval(indi,2))*1/si)';
      	end
      end
      
      time_gap_array = [time_gap_array; time_gap];
   end
end

% find the valid data points

if strcmp(pDummy,'NaN')
   index = find(~isnan(pData));
elseif pDummy > 0
   index = find(pData < pDummy);
elseif pDummy <= 0
   index = find(pData > pDummy);
end

% here is where filling in all untrackable segments, and creating a longer vector than 
% either u- or v- velocities (or props in the end) - HHF 28 oct 2003.  
% should we cull after the fact? 
time_spl      = linspace(pTime(1),pTime(size(pTime,1)),(pTime(size(pTime,1))-pTime(1))*1/si+1);  % new time basis in rafos days

pData   = pData(index); pTime = pTime(index);

data_spl_f    = csape(pTime,pData,[2 1]); 
data_spl      = fnval(data_spl_f,time_spl);
deriv_spl_f   = fnder(data_spl_f);
deriv_spl     = fnval(deriv_spl_f,time_spl) ; % unit is [data]/[time] 

%plot(time_spl,data_spl); hold on

% fill in the NaN's
if exist('time_gap_array', 'var')
	indi = find(ismember(round(time_spl(:)*24*60),round(time_gap_array(:)*24*60)));
	data_spl(indi)  = NaN*ones(size(indi)); 
	deriv_spl(indi) = NaN*ones(size(indi)); 
end

%plot(time_spl,data_spl,'g');
%plot(time_spl,data_spl,'gx');
%plot(time,data,'rx');