function adcp_whisk(theYears, dataDirectory, plotDirectory, siteArchive, theSite)
%
%   adcp_plots(theYears, dataDirectory, plotDirectory, siteArchive, theSite)
%
%   Function to create a whisker plot of the low passed wind stress and current
%   observations from the profiling current meters for theYears.
%   Files are output as SITE_whisk_YYYY.

%Soupy Alexander, 10/29/03

%Load in the wind stress
[windTime, tauU, tauV] = mblt_windstress(dataDirectory);

currTime_lp = [];
cU5_lp = [];
cV5_lp = [];
cU10_lp = [];
cV10_lp = [];
cU1_lp = [];
cV1_lp = [];

%Load in the site A data
for indexFile = 1:length(siteArchive)
    fileVars = siteArchive(indexFile).theVars;
    fileStart = siteArchive(indexFile).startTime;
    gregStart = gregorian(fileStart);
    fileEnd = siteArchive(indexFile).endTime;
    gregEnd = gregorian(fileEnd);
    if ~any(strcmp(fileVars, 'Werr_1201')) | ~any(ismember([gregStart(1) gregEnd(1)], theYears))
        continue
    end
    currFileA = netcdf(fullfile(dataDirectory, siteArchive(indexFile).fileList), 'nowrite');
    fileDepth = currFileA{'depth'}(:);
    fileTime = singlejd(currFileA{'time'}(:), currFileA{'time2'}(:));
    if length(fileDepth) == 1
        continue
    end
    
    fileU = currFileA{'u_1205'}(:);
    fileV = currFileA{'v_1206'}(:);
    bads = find(fileU > 1e3 | fileV > 1e3);
    fileU(bads) = NaN;
    fileV(bads) = NaN;
    
    hit5 = value2Index(fileDepth, 5);
    cU5 = fileU(:, hit5);
    try
        [fcU5_lp, fileTime_lp] = plfilt(cU5, fileTime, 6); 
    catch
        continue
    end
    cV5 = fileV(:, hit5);
    [fcV5_lp, fileTime_lp] = plfilt(cV5, fileTime, 6);
    
    if strcmp(lower(theSite), 'a')
        hit10 = value2Index(fileDepth, 15);
    elseif strcmp(lower(theSite), 'b')
        hit10 = value2Index(fileDepth, 10);
    end
    cU10 = fileU(:, hit10);
    [fcU10_lp, fileTime_lp] = plfilt(cU10, fileTime, 6); 
    cV10 = fileV(:, hit10);
    [fcV10_lp, fileTime_lp] = plfilt(cV10, fileTime, 6);
    
    hit1 = find(fileDepth == max(fileDepth));
    cU1 = fileU(:, hit1);
    [fcU1_lp, fileTime_lp] = plfilt(cU1, fileTime, 6); 
    cV1 = fileV(:, hit1);
    [fcV1_lp, fileTime_lp] = plfilt(cV1, fileTime, 6);
    
    [currTime_lp] = [currTime_lp(:); fileTime_lp(:)];
    [cU5_lp] = [cU5_lp(:); fcU5_lp(:)];
    [cV5_lp] = [cV5_lp(:); fcV5_lp(:)];
    [cU10_lp] = [cU10_lp(:); fcU10_lp(:)];
    [cV10_lp] = [cV10_lp(:); fcV10_lp(:)];
    [cU1_lp] = [cU1_lp(:); fcU1_lp(:)];
    [cV1_lp] = [cV1_lp(:); fcV1_lp(:)];
    
end

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    if (targetYear < 1996 & strcmp(lower(theSite), 'a')) | ...
            (targetYear < 1997 & strcmp(lower(theSite), 'b'))
        continue
    end
    
    figure, orient landscape
    
    xTick = [];
    for indexMonth = 1:12
        xTick = [xTick(:); datenum(targetYear, indexMonth, 1)];
    end
    %Do the plots
    timeVars = {'windTime'; 'currTime_lp'; 'currTime_lp'; 'currTime_lp'};
    uVars = {'tauU'; 'cU5_lp'; 'cU10_lp'; 'cU1_lp'};
    vVars = {'tauV'; 'cV5_lp'; 'cV10_lp'; 'cV1_lp'};
    yLimits = {[-4 4]; [-30 30]; [-30 30]; [-30 30]};
    yTicks = {[-4:2:4]; [-30:10:30]; [-30:10:30]; [-30:10:30]};
    yLabel = {'dynes/cm^2'; 'cm/s'; 'cm/s'; 'cm/s'};
    if strcmp(lower(theSite), 'a')
        plotTitles = {'Wind Stress'; 'Current (5 mbs)'; ...
                'Current (15 mbs)'; 'Current (deepest bin)'};
    else
        plotTitles = {'Wind Stress'; 'Current (5 mbs)'; ...
                'Current (10 mbs)'; 'Current (deepest bin)'};
    end
    
    for indexVar = 1:length(timeVars)
        eval(['hitTime = ' timeVars{indexVar} ';']); 
        eval(['hitU = ' uVars{indexVar} ';']);
        eval(['hitV = ' vVars{indexVar} ';']);
        inYear = find((hitTime >= julian(targetYear, 1, 1, 00)) & ...
            hitTime <= julian(targetYear, 12, 31, 24));
        hitTime = hitTime(inYear); 
        
        [theAngleRad, theLength] = cart2pol(hitU(inYear), hitV(inYear));
        theAngle = theAngleRad * 180/pi;
        
        goods = find(~isnan(theAngle));
        theAngle = theAngle(goods);
        theLength = theLength(goods);
        
        plotX = julian2datenum(hitTime(goods));
        plotY = zeros(size(plotX));
        
        subplot(4,1,indexVar)
        
        if isempty(plotX) & indexVar == 2 & targetYear >= 1998
            t = text(datenum(targetYear, 7, 1), 10, 'No instrumentation.');
            set(t, 'horizontal', 'center')
        end
        
        sticksafe(plotX, plotY, theLength, theAngle)
        ylim(yLimits{indexVar})
        set(gca, 'xtick', xTick, 'ytick', yTicks{indexVar})
        xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
                datenum(targetYear, 12, 31, 23, 59, 59)])
        datetick('x', 3, 'keeplimits', 'keepticks')
        ylabel(yLabel{indexVar})
        title(plotTitles{indexVar})
        l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
                datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
        set(l, 'color', 'k')
        box('on')
        sticksafe
    end
    
    suptitle({'Low-Pass Filtered Wind Stress and Profiling Current Meter Observations'; ...
            ['Site ' upper(theSite) ', ' num2str(targetYear, '%.0f')]})
    
    print(fullfile(plotDirectory, ...
        [lower(theSite) '_whisk_' num2str(targetYear, '%.0f')]), '-dpdf')
    
    close
end