function attn(theYears, dataDirectory, plotDirectory)
%
%   attn(theYears, dataDirectory, plotDirectory)
%
%   Function to create a plot of the low passed wind stress and beam
%   attenuation from Site A  for theYears.
%   Files are output as attn_YYYY.

%Soupy Alexander, 10/29/03

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

waveFile = netcdf(fullfile(dataDirectory, 'lt_wave.nc'), 'nowrite');
waveTime = singlejd(waveFile{'time'}(:), waveFile{'time2'}(:));
waveData = waveFile{'wh_4061'}(:);

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

cFile10 = netcdf(fullfile(dataDirectory, 'a_attn_10mab.nc'), 'nowrite');
cTime10 = singlejd(cFile10{'time'}(:), cFile10{'time2'}(:));
cU10 = cFile10{'ATTN_55'}(:);

cFile1 = netcdf(fullfile(dataDirectory, 'a_attn_1mab.nc'), 'nowrite');
cTime1 = singlejd(cFile1{'time'}(:), cFile1{'time2'}(:));
cU1 = cFile1{'ATTN_55'}(:);


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

[cU5_i] = smart_interp(cTime5, cU5, attnTime, 5);
[cU10_i] = smart_interp(cTime10, cU10, attnTime, 5);
[cU1_i] = smart_interp(cTime1, cU1, attnTime, 5);


for indexYear = 1:length(theYears)
    targetYear = theYears(indexYear);
    figure, orient landscape
    
    xTick = [];
    for indexMonth = 1:12
        xTick = [xTick(:); datenum(targetYear, indexMonth, 1)];
    end
    
    hitTime = windTime;
    hitU = tauU;
    hitV = tauV;
    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(5,1,1)
    sticksafe(plotX, plotY, theLength, theAngle)
    ylim([-4 4])
    set(gca, 'xtick', xTick, 'ytick', [-4:2:4])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('dynes/cm^2')
    t = title('Low-passed Wind Stress');   
    set(t, 'fontsize', 8)
    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
    
    subplot(5,1,2)
    plot(julian2datenum(waveTime), waveData)
    ylim([0 10])
    set(gca, 'xtick', xTick, 'ytick', [0:2:10])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('m')
    t = title('Hourly Averaged Significant Wave Height');   
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    subplot(5,1,3)
    if targetYear >= 1998
        t = text(datenum(targetYear, 7, 1), 2.5, 'No instrumentation.');
        set(t, 'horizontal', 'center')
    else    
        plot(julian2datenum(attnTime), cU5_i)
    end
    ylim([0 5])
    set(gca, 'xtick', xTick, 'ytick', [0:2:10])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('m^-^1')
    t = title('Beam Attenuation:  5 mbs');
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    subplot(5,1,4)
    plot(julian2datenum(attnTime), cU10_i)
    ylim([0 5])
    set(gca, 'xtick', xTick, 'ytick', [0:2:10])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('m^-^1')
    t = title('Beam Attenuation:  10 mab');    
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    subplot(5,1,5)
    plot(julian2datenum(attnTime), cU1_i)
    ylim([0 5])
    set(gca, 'xtick', xTick, 'ytick', [0:2:10])
    xlim([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)])
    datetick('x', 3, 'keeplimits', 'keepticks')
    ylabel('m^-^1')
    t = title('Beam Attenuation:  1 mab'); 
    set(t, 'fontsize', 8)
    l = line([datenum(targetYear, 1, 1, 00, 00, 00) ...
            datenum(targetYear, 12, 31, 23, 59, 59)], [0 0]);
    set(l, 'color', 'k')
    box('on')
    
    
    suptitle({'Low-Pass Filtered Wind Stress, Hourly Averaged Significant Wave Height, and Beam Attenuation'; ...
        ['Site A, ' num2str(targetYear, '%.0f')]})
    
    print(fullfile(plotDirectory, ['attn_' num2str(targetYear, '%.0f')]), ...
        '-dpdf')

    close
end