function dout = degapint(din,flag,delt,tlims,tmax)
% Function to degap by interpolation a randomonly sampled time series with both small and large gaps
% Create February 9, 2001 for ISSIBook and ESS 204
% MODS: March 25, 2001, July 7, 2001
% din = Input array with with time as datenum in column 1 and a set of several variables in other cols
% flag = Flag in the input data
% delt = The desired sample rate in the output file (IN DATE NUMBER)
% tlims = Desired start and end datenumbers [tb, te] (IN DATE NUMBER)
% tmax = The largest duration of a gap that will be filled by interpolation (IN DATE NUMBER)
% dout = The output array with all gaps replaced by flagged records

% ELIMINATE ANY RECORDS WITH FLAGGED TIMES (These come from propagating the solar wind using vel measurments)
is = find(din(:,1)==flag); %Locs of flags in time column
din(is,:) = []; %Dump these records
fprintf(1,'Number of flagged time records = %d\n',length(is))

% SORT THE RECORDS IN ASCENDING TIME ORDER (my backwards time finder did not work)
[dn,id] = sort(din(:,1)); %Sort time column in ascending order
din = din(id,:); %Put all records in order of ascending time
% ELIMINATE NO ADVANCE TIME RECORDS TO GET MONOTONIC TIME VARIABLE
dn = din(:,1); %Extract the input time column
dt = diff(dn); %First difference of the time variable
is = find(dt==0); %No time change (next record at same time)
if(~isempty(is))
dn(is+1) = []; %Eliminate the bad record in the time column
din(is+1,:) = []; %Eliminate the bad record in the input arrayfor i =1:length(is)
end
fprintf(1,'Number of no time change records = %d\n',length(is))

noskip = 0;
if(noskip)
% FIND THE BEST POSSIBLE SAMPLE INTERVAL
dt = diff(din(:,1)); %Create the time between samples in the input data (in decimal days)
medt = 86400*median(dt); %Find the median sample interval (in seconds)
fprintf('median delt(sec) = %12.7f\n',medt) %This is likely to be very close to the actual sample rate
nuv = findnuv(86400*dt,flag); %Create a table of unique separations of samples, and their occ rate
pnuv = nuv';
fprintf(1,'%12.6f %7d\n',pnuv) %Print out the beginning of the table
disp('Compare your desired sample interval and maximum gap to this table')
disp(' Terminate the program if your choices seem wrong')
else
disp('Skipped the determination of unique delta time')
end

% CALCULATE THE NUMBER OF POINTS EXPECTED WITH INPUT SAMPLE RATE
[md, nd]= size(din); %Get the size of the input array
db = tlims(1); %The desired start time of output array
de = tlims(2); %The desired end time of output array
npts = round((de-db)/delt)+1; %Total number of points in requested interval
fprintf(1,'Input array has %7d rows\n',md)
fprintf(1,'Output array has %7d rows\n',npts)

% CREATE AN OUTPUT ARRAY TO HOLD THE DEGAPPED DATA
dout = repmat(flag,npts,nd); %Initialize output array to flags
dout(:,1) = linspace(db,de,npts)'; %Create the output time column

% LOCATE LARGE TIME GAPS IN THE INPUT DATA ARRAY
dt = diff(dn); %Recalculate first diff because of above changes in dn
is = find(dt>tmax); %Find all first diffs that are longer than max allowable interpolation
is =[0;is;length(dn)]; %Augment ends of is array to deal with initial and final data segments
fprintf(1,'Number of large gaps in time sequence = %d\n',length(is)-2)

% LOOP FOR EVERY SEGMENT BETWEEN GAPS
nmin = 4; %Need this for interpolation
mseg = round(tmax/delt); %Largest number of flags to interpolate
for i = 1:(length(is)-1)
jseg = is(i)+1:is(i+1); %This is segment to the right of the gap
if(length(jseg)>=nmin) %May be enough points to do interpolation
%Do an interpolation over the segment.
for j = 2:nd %Loop for every column except time
dat = din(jseg,j); %dat is a single column array of data with possible flag segments
tdn = dn(jseg); %tdn is a single column array of datenumber coresponding to dat
%Note that there are no time gaps longer than tmax in this segment
%However, there may be a combination of time gaps and flag intervals that are longer than tmax
%I am currently ignoring this (Aug 2, 2001)
% THIS SECTION IS TO INTERPOLATE TO NEW RESOLUTION, BUT TO AVOID INTERVALS OF FLAGS LONGER THAN TMAX
jbjejr = findinterp(dat,flag,mseg);%Find start and stop points of subsegs in seg dat to interpolate
[ml,nl] = size(jbjejr); %ml is the # of iterations in the interpolation loop for current seg
for k = 1:ml %Loop for each segment that can be interpolated
js = jbjejr(k,1):jbjejr(k,2); %Indices of subsegment in segment dat
jk = find(dat(js)==flag); %Find the flags in the subseg
js(jk) = []; %Discard flagged points in subsegment
if(~isempty(js)) %Can't do these if there are no points left
in = find(dout(:,1)>=(tdn(js(1))-delt) & dout(:,1)<=(tdn(js(end))+delt)); %Indicies of output points in subseg
kk = find(in<1 | in>length(dout(:,1)));
in(kk) =[]; %Null points outside of array
if(length(js)>=nmin) %Still enough points to do interpolation
dout(in,j) = interp1(tdn(js),dat(js),dout(in,1),'linear','extrap'); %Do a linear fit to non flag points
else %Give up and transfer
%Transfer the points for this component to closest time match in output
in = round((tdn(js)-db)/delt)+1; %Index in output array
dout(in,j) = din(js,j); %Map to closest time match
end
end
end
end
end
%%keyboard
if(length(jseg)<nmin) %There are not enough points for interpolation
%Transfer the points of all components to closest time match in output
in = round((dn(jseg)-db)/delt)+1; %Index in output array
dout(in,2:nd) = din(jseg,2:nd); %Map to closest time match
end
end
return