function pcm_whisk(theYears, dataDirectory, plotDirectory)
%
%   pcm_whisk(theYears, dataDirectory, plotDirectory)
%
%   Function to create a plot of the low passed wind stress and current
%   observations from the Site A point current meters for theYears.
%   Files are output as pcm_whisk_YYYY.

%Soupy Alexander, 10/29/03

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

%Load in the concatenated data from 5 mbs, 10 mab, and 1 mab
cFile5 = netcdf(fullfile(dataDirectory, 'a_curr_5mbs.nc'), 'nowrite');
cTime5 = singlejd(cFile5{'time'}(:), cFile5{'time2'}(:));
cU5 = cFile5{'u_1205'}(:);
cV5 = cFile5{'v_1206'}(:);

cFile10 = netcdf(fullfile(dataDirectory, 'a_curr_10mab.nc'), 'nowrite');
cTime10 = singlejd(cFile10{'time'}(:), cFile10{'time2'}(:));
cU10 = cFile10{'u_1205'}(:);
cV10 = cFile10{'v_1206'}(:);

cFile1 = netcdf(fullfile(dataDirectory, 'a_curr_1mab.nc'), 'nowrite');
cTime1 = singlejd(cFile1{'time'}(:), cFile1{'time2'}(:));
cU1 = cFile1{'u_1205'}(:);
cV1 = cFile1{'v_1206'}(:);

currTime = min([cTime5(:); cTime10(:); cTime1(:)]):1/24:...
    max([cTime5(:); cTime10(:); cTime1(:)]);

[cU5_i] = smart_interp(cTime5, cU5, currTime, 5);
[cU5_lp, currTime_lp] = plfilt(cU5_i, currTime, 6);

[cV5_i] = smart_interp(cTime5, cV5, currTime, 5);
[cV5_lp, currTime_lp] = plfilt(cV5_i, currTime, 6);

[cU10_i] = smart_interp(cTime10, cU10, currTime, 5);
[cU10_lp, currTime_lp] = plfilt(cU10_i, currTime, 6);

[cV10_i] = smart_interp(cTime10, cV10, currTime, 5);
[cV10_lp, currTime_lp] = plfilt(cV10_i, currTime, 6);

[cU1_i] = smart_interp(cTime1, cU1, currTime, 5);
[cU1_lp, currTime_lp] = plfilt(cU1_i, currTime, 6);

[cV1_i] = smart_interp(cTime1, cV1, currTime, 5);
[cV1_lp, currTime_lp] = plfilt(cV1_i, currTime, 6);

%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'};
plotTitles = {'Wind Stress'; 'Current (5 mbs)'; ...
        'Current (10 mab)'; 'Current (1 mab)'};

for indexYear = 1:length(theYears)
    targetYear = theYears(indexYear);
    figure, orient landscape
    
    xTick = [];
    for indexMonth = 1:12
        xTick = [xTick(:); datenum(targetYear, indexMonth, 1)];
    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 Point Current Meters Observations'; ...
        num2str(targetYear, '%.0f')})
    
    print(fullfile(plotDirectory, ['pcm_whisk_' num2str(targetYear, '%.0f')]), ...
        '-dpdf')

    close
end