% fixbadadvvel - use the results from flagbadadv.m to remove bad blocks of data
% in the burst file and then recalculate the statistics.
% This program operated directly on the files provided and does not make copies.  
% fixbadadv(dirtybFile, dirtysFile, qualityFile)

% Written by Marinna Martini for the U.S. Geological Survey 
% Coastal & Marine Program Woods Hole Field Center, Woods Hole, MA
% http://woodshole.er.usgs.gov/ Please report bugs to mmartini@usgs.gov

% 21 oct 2005 -- V 2.0 make this work with the new Quality file
% 29 Dec 2004 -- replace nanmean and nanstd with gmean and gstd so that
% users don't need the stats toolbox

function fixbadadvvel(dirtybFile, dirtysFile, qFile)

prog_ver = 2.0;

% those pesky input checks....
if nargin < 3 || nargin > 3,
    disp('Input information is incorrect, use ''help fixbadadv'' for help')
    return
end

cdfq = netCDF(qFile,'write');
if isempty(cdfq),
    disp(sprintf('Could not open netCDF quality file %s for writing',qFile))
    return
end
% make sure there's something to do
mask = cdfq{'flagbadadv_mask'}(:,:,:);
if isempty(find(mask,1)), % there's nothing to do!
    disp('No bad data to replace, aborting')
    close(cdfq)
    return
end

cdfb = netCDF(dirtybFile,'write');
if isempty(cdfb),
    disp(sprintf('Could not open netCDF burst data file %s for writing',dirtybFile))
    return
elseif ~strcmp(cdfb.DESCRIPT(:),'Sontek ADV raw data burst file'),
    disp(sprintf('Could not open netCDF burst data file %s',dirtybFile))
    close(cdfb)
    return
end
cdfs = netCDF(dirtysFile,'write');
if isempty(cdfs),
    disp(sprintf('Could not open netCDF burst data file %s for writing',dirtysFile))
    return
elseif ~strcmp(cdfs.DESCRIPT(:),'Sontek ADV raw data statistics file'),
    disp(sprintf('Could not open netCDF statistics data file %s',dirtysFile))
    close(cdfs)
    return
end

% the shape is Qdata.allbads{iBeam}{subburst}[badburst1, badburst2, badburst3,...]
dims = size(mask);
nBeams = dims(2);
nSubbursts = dims(3);
nBursts = dims(1);
nSamples = cdfb.ADVDeploymentSetupSamplesPerBurst(:);
nSubsamp = floor(nSamples/nSubbursts)
indSubb(:,1) = [1:nSubsamp:nSamples]'; % starting indeces for each subburst
if rem(nSamples, nSubbursts), 
    indSubb(:,2) = [nSubsamp:nSubsamp:(nSubsamp*nSubbursts)+nSubsamp]'; % starting indeces for each subburst
else
    indSubb(:,2) = [nSubsamp:nSubsamp:nSamples]'; % starting indeces for each subburst
end
indSubb = indSubb(1:nSubbursts,:);
% we're going to replace the entire subburst with NaN, so make a convenient
% NaN holder to be efficient
theFill = ones(nSubsamp,1).*NaN;
% which variables are we editing?
theVars = {'Velx','Vely','Velz'};
% define a counter for the number of good samples in a burst
cdfq{'fixbadadvvel_NGoodVelSamples'} = ncfloat('burst','axis'); ncobj = cdfq{'fixbadadvvel_NGoodVelSamples'};
ncobj.long_name = ncchar('Count of Good Samples in Burst');
ncobj.units = ncchar('counts');
ncobj.valid_range = ncfloat([0 nSamples]);
ncobj.comment = ncchar('Calculated by fixbadadvvel.m');
% strat by assuming they are all good...
cdfq{'fixbadadvvel_NGoodVelSamples'}(:,:) = ones(nBursts,nBeams).*nSamples;
burstNum = cdfq{'burst'}(:);

% we want to remove all three beams' data if any one is bad.
% the main loop
for iBeam = 1:nBeams,
    for iSubburst = 1:nSubbursts,
        nbadBursts = length(find(mask(:,iBeam,iSubburst)));        
        if ~isempty(nbadBursts),
            disp(sprintf('There are %d bad bursts in subburst %d for beam %d', nbadBursts, iSubburst, iBeam))
            % note iBurst here is the index to the burst in the file, not the burst number
            iBurst = find(mask(:,iBeam,iSubburst));
            for idx=1:length(iBurst), % this is not going through every burst, rather the list of bad ones
                % we need the index
                for iVar = 1:length(theVars),
                    % replace the velocity with the fill value
                    % TODO there are the wrong indeces here... fix it...
                    cdfb{theVars{iVar}}(iBurst(idx),indSubb(iSubburst,1):indSubb(iSubburst,2)) = theFill;
                    % now get the whole, corrected burst
                    data = cdfb{theVars{iVar}}(iBurst(idx),:);
                    % recalculate Mean and STD of the variables for the burst
                    % store them in the stats file
                    theName = sprintf('Mean%s',theVars{iVar});
                    cdfs{theName}(iBurst(idx)) = gmean(data);
                    theName = sprintf('Std%s',theVars{iVar});
                    cdfs{theName}(iBurst(idx)) = gstd(data);
                    ngood = nSamples - length(find(isnan(data)));
                    cdfq{'fixbadadvvel_NGoodVelSamples'}(iBurst(idx),iBeam) = ngood;
                end
            end
        else disp(sprintf('There are no bad bursts in subburst %d', iSubburst))
        end
    end
end

hstring = cdfs.history(:);
hstring = sprintf('edited by %s %f; %s', mfilename, prog_ver, hstring);
cdfs.history = hstring;
hstring = cdfb.history(:);
hstring = sprintf('edited by %s %f; %s', mfilename, prog_ver, hstring);
cdfb.history = hstring;

close(cdfb)
close(cdfs)
close(cdfq)

disp('Fixbadadv done')