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 % dout = degapint(din,flag,delt,tlims,tmax); % 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 input data either before or after the interval tlims ii = find(din(:,1)tlims(2)); din(ii,:) = []; % 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 if(length(jseg)