% progressADVdeadrecrep see if wake from leg gets to ADV, repeating through a burst
%
% use progressive vectors to dead recon the end point of a particle N
% diameters from a leg (or other obstruction)
%
% function [badBursts badCount] = progressADVdeadrecrep(settings)
% settings:
%     settings.ADVFile = {'file1'....}; % cell array of ADV burst data
%     settings.ADVStats = {'file1'...}; % cell array of ADV statistics files
%           from one ADV that matches time span of settings.ADVFile input
%     settings.posleg{3}; [x, y] combinations
%     settings.posADV; [x, y] combinations
%     settings.verbose; %generate plots, good for making movies
%     settings.nd; % number of diameters
%     settings.diameter = 0.085 cm% diameter of obstructions
%     settings.wakewidth = 20; % wakes must be within +/- this angle to ADV
%     settings.ADVarea = 0.5 m; % wakes must get within this distance of ADV
%     settings.repinterval = 10; %sec
%
% badbursts = a cell array {ADV1, ADV2, ADVN}
%   where ADVN = [nburst, nleg]
%       if the ADV was shadowed by a leg during that burst, the ADV burst
%       number is recorded in the column for that leg.
%       for example, if ADVN was shadowed by leg #2 in posleg{3} at ADV
%       file index i, ADVN(i,2) = ADV{'burst'}(i) or the burst number at
%       file index i
%   thus, i = find(badbursts{ADVN}(:,2)) will return the indeces into the file of all
%   bursts from ADVN shadowed by leg #2 and badbursts{ADVN}(i,2) will
%   return the Sontek burst numbers of those bursts.
%

% Written by Marinna Martini for the U.S. Geological Survey
% Coastal and Marine Program Woods Hole Science Center, Woods Hole, MA

function [badBursts badCount] = progressADVdeadrecrep(settings)

%prog_ver = 'SVN $Revision:  $';

if ~exist('settings','var'), settings = []; end

if ~isstruct(settings),
    ADVFile = settings;
    clear settings
    settings.ADVFile = ADVFile;
end
if ~isfield(settings,'ADVFile'), return; end
if ~isfield(settings,'posleg'),
    disp('need leg position');
    return;
else
    nlegs = length(settings.posleg);
    legx = zeros(1,nlegs);
    legy=legx;
    for ileg = 1:nlegs,
        legx(ileg) = settings.posleg{ileg}(1);
        legy(ileg) = settings.posleg{ileg}(2);
    end
end
if ~isfield(settings,'posADV'), disp('need ADV position'); return; end
if ~isfield(settings,'verbose'), settings.verbose = 1; end
if ~isfield(settings,'nd'), settings.nd = 40; end
if ~isfield(settings,'diameter'), settings.diameter = 0.085; end
if ~isfield(settings,'wakewidth'), settings.wakewidth = 20; end
if ~isfield(settings,'ADVarea'), settings.ADVarea = 0.5; end
if ~isfield(settings,'repinterval'), settings.repinterval = 10; end %sec

nADVs = length(settings.ADVFile);

cdfin = cell(2);
for iADV = 1:nADVs,
    cdfin{iADV} = netcdf(settings.ADVFile{iADV});
    disp(sprintf('%s start time %s end time %s',cdfin{iADV}.start_time(:),...
        cdfin{iADV}.stop_time(:)))
    % make sure files are aligned in time
end

% look for a stats file
if isfield(settings,'ADVstats'),
    for iADV = 1:length(settings.ADVstats)
        cdfstat{iADV} = netcdf(settings.ADVstats{iADV});
    end
else
    cdfstat{iADV} = [];
end

% figure out some important things based on the blocks and samples
% calculate the harmonics
%nSamples = length(cdfin('sample')); % samples per ADV burst
nBursts = length(cdfin{1}('time'));
nSamples = length(cdfin{1}('sample'));
verbose = settings.verbose;
dt = 1/cdfin{1}.ADVDeploymentSetupSampleRate(1); % ADV sampling rate, Hz
pipedia = settings.diameter; %0.084; % 8.4 cm diameter leg size
wakewidth = settings.wakewidth; % wakes must be within +/- this angle to ADV
ADVarea = settings.ADVarea; % wakes must get within this distance of ADV
dist = settings.nd*pipedia;
badBursts = cell(nADVs,1);
badCount = badBursts;
u=badBursts;
v=badBursts;
vdir=badBursts;
vspeed=badBursts;
for iADV = 1:nADVs
    badBursts{iADV} = ones(nBursts,nlegs).*NaN;
    badCount{iADV} = badBursts{iADV};
    u{iADV} = ones(nSamples,1).*NaN;
    v{iADV} = u{iADV};
    vdir{iADV} = u{iADV};
    vspeed{iADV} = u{iADV};
end
dirADV = cell(nADVs, nlegs);
distADV = dirADV;
repinterval = settings.repinterval; %sec
repsamples = repinterval/dt;
nreps = nSamples/repsamples;
disp(sprintf('Progressing in each burst by %d samples = %d seconds',repsamples,repinterval))

% relation of ADVs from legs, polar for Marinna's method
for iADV = 1:nADVs,
    for ileg = 1:nlegs,
        [dirADV{iADV, ileg}, distADV{iADV, ileg}] = ...
            uv2polar(settings.posADV{iADV}(1)-legx(ileg),...
            settings.posADV{iADV}(2)-legy(ileg));
        disp(sprintf('The %s ADV is %f degrees and %f m from the %s leg',...
            settings.ADVcolor{iADV},dirADV{iADV, ileg},distADV{iADV, ileg},settings.legname{ileg}))
    end
    badBursts{iADV} = zeros(nBursts,nlegs);
end

dummydata = zeros(1,cdfin{1}.ADVDeploymentSetupSamplesPerBurst(1));

mkdir('frames')

figure('position',[ 183         327        1031         595])
ax1 = subplot(3,1,3);
if ~isempty(cdfstat{1}),
    tm = datenum(gregorian(cdfstat{1}{'time'}(:)+cdfstat{1}{'time2'}(:)./(1000*24*3600)));
    plot(tm,cdfstat{1}{'SDP_850'}(:));
    hold on
    for iADV = 2:length(cdfstat),
        tj = cdfstat{iADV}{'time'}(:)+cdfstat{iADV}{'time2'}(:)./(1000*24*3600);
        plot(datenum(gregorian(tj)),cdfstat{iADV}{'SDP_850'}(:));
        tm = [tm; datenum(gregorian(tj))];
    end
    datetick('x',6,'keepticks','keeplimits');
    set(gca,'ylim',[0 50])
    ylabel(cdfstat{1}{'SDP_850'}.units(:))
    xlabel(cdfstat{1}{'SDP_850'}.long_name(:))
    if ischar(cdfstat{1}.MOORING(:)),
        text(0.1, -0.2,sprintf('Tripod %s',cdfstat{1}.MOORING(:)),'units','normalized');
    else
        text(0.1, -0.2,sprintf('Tripod %d',cdfstat{1}.MOORING(1)),'units','normalized');
    end
    y = get(gca,'ylim');
    bar = plot([tm(1) tm(1)],y,'w'); % hidden priming bar
    hold off
else
    text(0.5,0.5,'No stats data','units','normalized')
end

dx = cell(nADVs,1); dy = dx;
uc = dx; vc = dx;
dirc = dx; distc = dx;
wake(1:nADVs) = struct('speed',0,'time',0,'samples',0);
ax2 = subplot(3,1,[1 2]);
%ax3 = subplot(3,2,[2 4]);
h = zeros(nADVs,2);
tic
for iBurst = 1:nBursts,
    burstDate = datestr(datenum(gregorian(cdfin{1}{'time'}(iBurst)+cdfin{1}{'time2'}(iBurst)./(1000*24*3600))));
    burstNum = cdfin{1}{'burst'}(iBurst);
    project = cdfin{1}.PROJECT(:);
    disp(sprintf('Burst %d at file index #%d',burstNum,iBurst))

    %get the data
    for iADV = 1:nADVs,
        %if burstNum == 111, keyboard; end
        u{iADV} = cdfin{iADV}{'u_1205'}(iBurst,:);
        v{iADV} = cdfin{iADV}{'v_1206'}(iBurst,:);
        if any(isnan(u{iADV})) || any(isnan(v{iADV})) || all(~u{iADV}) || all(~v{iADV}),
            u{iADV} = dummydata;
            v{iADV} = dummydata;
            disp(sprintf('\tNo %s velocity data', settings.ADVcolor{iADV}))
        end

        [vdir{iADV}, vspeed{iADV}] = uv2polar(u{iADV},v{iADV});
        wake(iADV).speed = gmean(vspeed{iADV})./100;
        if wake(iADV).speed < 0.001, % got to be bogus, ADV can't measure this small
            % now keep it from screwing up things downstream
            wake(iADV).speed = 0;
            wake(iADV).time = 0;
            wake(iADV).samples = 1;
        else
            wake(iADV).time = dist/wake(iADV).speed;
            wake(iADV).samples = floor(wake(iADV).time/dt);
        end
    end

    % now for each ADV from each leg see if it shadows by each ADV
    for sADV = 1:nADVs, % this is the ADV we are checking for shadowing
        for ileg = 1:nlegs
            for iADV = 1:nADVs, % this is the ADV whose data we are using
                if verbose, disp(sprintf('\tChecking if %s is shadowed by %s data from %s leg',...
                    settings.ADVcolor{sADV}, settings.ADVcolor{iADV}, settings.legname{ileg})); end
                if 1,
                    for irep = 0:nreps-1, % do it repreatedly throught the burst
                        s1 = irep*(repinterval./dt)+1;
                        s2 = s1+wake(iADV).samples-1;
                        % this here is dumping lots of potential shadowing
                        if s2<length(u{iADV}),
                            dx{iADV} = u{iADV}(s1:s2).*dt./100; % incremental displacement in X direction (cm/s)*s*(1/100 m/cm) = m
                            dy{iADV} = v{iADV}(s1:s2).*dt./100; % incremental displacement in Y direction

                            % do a cumulative distance of u and v in meters
                            uc{iADV} = cumsum(dx{iADV});
                            vc{iADV} = cumsum(dy{iADV});
                            [dirc{iADV}, distc{iADV}] = uv2polar(uc{iADV}(end)-uc{iADV}(1),vc{iADV}(end)-vc{iADV}(1));
                            % mathematically see if it's near the ADVs
                            % from the leg, see if it gets near each ADV
                            % if we are within +/- 10 degrees of the bearing of the
                            % ADV from the leg, check the range
                            dirdiff = abs(dirc{iADV}-dirADV{sADV, ileg});
                            if (dirdiff <= wakewidth) && (distc{iADV} > distADV{sADV, ileg}-ADVarea),
                                % we are impinging on the ADV, note it.
                                badBursts{sADV}(iBurst, ileg) = burstNum; % set to burstnum for now.
                                badCount{sADV}(iBurst, ileg) = badCount{sADV}(iBurst, ileg)+1;
                                if verbose, disp(sprintf('\t\t%s shadowed by %s leg in rep %d',settings.ADVcolor{sADV},settings.legname{ileg},irep)); end
                            end
                            if 0, %(burstNum == 52) & verbose,  diagnostic check
                                disp(sprintf('\t\trep %d: dirdiff = %f, wakewidth = %f, distc = %f, distADV = %f',...
                                    irep, dirdiff, wakewidth, distc{iADV}, distADV{sADV,ileg}))
                            end
                        end
                    end
                end
                if 0,
                    dx{iADV} = u{iADV}.*dt./100; % incremental displacement in X direction (cm/s)*s*(1/100 m/cm) = m
                    dy{iADV} = v{iADV}.*dt./100; % incremental displacement in Y direction

                    % do a cumulative distance of u and v in meters
                    uc{iADV} = cumsum(dx{iADV});
                    vc{iADV} = cumsum(dy{iADV});
                    [dirc{iADV}, distc{iADV}] = uv2polar(uc{iADV}(end)-uc{iADV}(1),vc{iADV}(end)-vc{iADV}(1));
                    % mathematically see if it's near the ADVs
                    % from the leg, see if it gets near each ADV
                    % if we are within +/- 10 degrees of the bearing of the
                    % ADV from the leg, check the range
                    dirdiff = abs(dirc{iADV}-dirADV{sADV, ileg});
                    if (dirdiff <= wakewidth) && (distc{iADV} > distADV{sADV, ileg}-ADVarea),
                        % we are impinging on the ADV, note it.
                        badBursts{sADV}(iBurst, ileg) = burstNum; % set to burstnum for now.
                        badCount{sADV}(iBurst, ileg) = badCount{sADV}(iBurst, ileg)+1;
                        disp(sprintf('\t\t%s shadowed by %s leg',settings.ADVcolor{iADV},settings.legname{ileg}))
                    end
                end
                if verbose, disp(sprintf('\t\tdirdiff = %f, distc = %f, distADV = %f',...
                    dirdiff(end), distc{iADV}(end), distADV{sADV,ileg})); end
            end
        end
    end

    %if burstNum == 52, keyboard; end
    axes(ax2)
    cla
    plot(legx,legy,'o');
    text(legx+0.1,legy,settings.legname)
    lim = 4;
    set(gca,'ylim',[-lim lim])
    set(gca,'xlim',[-lim lim])
    hold on
    for iADV = 1:nADVs,
        plot(settings.posADV{iADV}(1),settings.posADV{iADV}(2),...
            [settings.ADVcolor{iADV}(1) 'v'])
        text(settings.posADV{iADV}(1)+0.1,settings.posADV{iADV}(2),settings.ADVcolor{iADV})
        text(0.01, 0.6-0.05*iADV*2,sprintf('%s shadowed by',settings.ADVcolor{iADV}), 'units', 'normalized')
        for ileg = 1:length(legx),
            if badBursts{iADV}(iBurst,ileg),
                text(0.02+0.05*ileg, 0.55-0.05*iADV*2,sprintf('%s',settings.legname{ileg}), 'units', 'normalized')
            end
        end
        if any(u{iADV}) && any(v{iADV}),
            for ileg = 1:length(legx),
                h(ileg,1) = plot(settings.posleg{ileg}(1)+uc{iADV},settings.posleg{ileg}(2)+vc{iADV},...
                    [settings.ADVcolor{iADV}(1) '-']);
                h(ileg,2) = plot(settings.posleg{ileg}(1)+uc{iADV}(end),settings.posleg{ileg}(2)+vc{iADV}(end),...
                    [settings.ADVcolor{iADV}(1) 'x']);
            end
            text(0.01, 0.05*iADV,sprintf('%s Mean east %3.2f, north %3.2f, m/s, dist %3.2f m at %d degrees',...
                settings.ADVcolor{iADV}, gmean(u{iADV}(1:wake(iADV).samples))./100,gmean(v{iADV}(1:wake(iADV).samples))./100,...
                distc{iADV}(end), floor(dirc{iADV}(end))), 'units', 'normalized')
            text(0.01, 1-0.05*iADV,sprintf('%s %d diameters = %d samples in %f sec',...
                settings.ADVcolor{iADV}, settings.nd, wake(iADV).samples, wake(iADV).time), 'units', 'normalized')
        else
            text(0.01, 0.05*iADV,sprintf('%s no data',settings.ADVcolor{iADV}), 'units', 'normalized')
            text(0.01, 1-0.05*iADV,sprintf('%s no data',settings.ADVcolor{iADV}), 'units', 'normalized')
        end
    end
    text(0.01, 1-0.05*nADVs-0.05,sprintf('Testing for wake width of %d deg. and %2.1f m area around ADV',...
        wakewidth, ADVarea), 'units', 'normalized')
    text(0.01, 1-0.05*nADVs-0.1,sprintf('%d sec repeat interval',...
        settings.repinterval), 'units', 'normalized')
    text(0.01, 1-0.05*nADVs-0.15,sprintf('%d diameters = %4.1f m',...
        settings.nd, dist), 'units', 'normalized')
    ylabel('North-South, m')
    xlabel('East-West, m')
    title(sprintf('%s Progressive Vectors for burst %d on %s',project,burstNum,burstDate))
    hold off
    axes(ax1)
    hold on
    delete(bar)
    idx = find(tm >= datenum(burstDate),1,'first');
    bar = plot([tm(idx) tm(idx)],y,'r');
    hold off

    % save plot
    drawnow
    eval(sprintf('print -dpng .\\frames\\b%d.png',burstNum));
end
if ~mod(iBurst,10), disp(sprintf('Finished %d bursts %4.2f min elapsed',iBurst,toc./60)), end

for iADV = 1:nADVs
    close(cdfin{iADV})
end
for iADV = 1:length(cdfstat)
    close(cdfstat{iADV})
end

return