function mixed_plot(theYears, dataDirectory, plotDirectory)
%
%   mixed_plot(theYears, dataDirectory, plotDirectory)
%
%   Creates the plots of significant wave height, psdev, beam attenuation,
%   sediment trapping rate, and current speed (1 mab) at Site A.  Plots are
%   named a_mix_YYYY.

%Soupy Alexander, 11/12/03

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

%Standard deviation of pressure
sdpFile = netcdf(fullfile(dataDirectory, 'a_sdp_1mab.nc'), 'nowrite');
sdpTime = singlejd(sdpFile{'time'}(:), sdpFile{'time2'}(:));
sdpData = sdpFile{'SDP_850'}(:);

%Beam Attenuation
attnFile = netcdf(fullfile(dataDirectory, 'a_attn_1mab.nc'), 'nowrite');
attnTime = singlejd(attnFile{'time'}(:), attnFile{'time2'}(:));
attnData = attnFile{'ATTN_55'}(:);

%Current Speed
currFile = netcdf(fullfile(dataDirectory, 'a_curr_1mab.nc'), 'nowrite');
currTime = singlejd(currFile{'time'}(:), currFile{'time2'}(:));
currU = currFile{'u_1205'}(:);
currV = currFile{'v_1206'}(:);
currData = sqrt(currU.^2 + currV.^2);

%Sediment trap data
[theNums, theText] = xlsread(fullfile(dataDirectory, 'all_traps.xls'), ...
    'LNB Traps');
theNums = [repmat(NaN, 1, size(theText, 2)); theNums];
trapType = theText(:,19);
for index = 1:length(trapType)
    if isempty(trapType{index})
        trapType(index) = {' '};
    end
end
theMAB = theNums(:,23);
theHonjos = find(strcmp('honjo', lower(trapType)) & (theMAB <= 4.5));

honjoStart = theNums(theHonjos,4) + datenum('30-Dec-1899');
honjoEnd = theNums(theHonjos, 6) + datenum('30-Dec-1899');
honjoFlux = theNums(theHonjos, 46);

theXs = [];
theYs = [];
allCol = {};

for indexTrap = 1:length(honjoFlux)
    if ~isnan(honjoFlux(indexTrap))
        thisX = [honjoStart(indexTrap) honjoStart(indexTrap) ...
                honjoEnd(indexTrap) honjoEnd(indexTrap)];
        thisY = [0 honjoFlux(indexTrap) honjoFlux(indexTrap) 0];
        if honjoFlux(indexTrap) >= 20
            thisCol = 'r';
        else
            thisCol = 'b';
        end
        theXs = [theXs thisX(:)];
        theYs = [theYs thisY(:)];
        allCol = [allCol {thisCol}];
    end
end

plotVars = {'wave'; 'sdp'; 'attn'; 'curr'};
plotTitles = {'Significant Wave Height'; 'Pressure Standard Deviation (1 mab)'; ...
        'Beam Attenuation (1 mab)'; 'Current Speed (1 mab)'};
plotNums = [1 2 3 5];
yRanges = {[0 10]; [0 50]; [0 5]; [0 30]};
yLabels = {'m'; 'mb'; 'm^-^1'; 'cm/s'};

theLabels = {'Jan'; 'Feb'; 'Mar'; 'Apr'; 'May'; 'Jun'; 'Jul'; 'Aug'; 'Sep'; 'Oct'; 'Nov'; ...
        'Dec'};

for indexYear = 1:length(theYears)
    targetYear = theYears(indexYear);
    
    tickLoc = [];
    for indexMonth = 1:12
        tickLoc = [tickLoc; datenum(targetYear, indexMonth, 1)];
    end
    
    figure, orient landscape
    
    for indexVar = 1:length(plotVars);
        eval(['theTime = julian2datenum(' plotVars{indexVar} 'Time);']);
        eval(['theData = ' plotVars{indexVar} 'Data;']);
        
        [plotX, plotY] = nanfill_2(theTime, theData);
        
        subplot(5,1,plotNums(indexVar))
        set(gca, 'fontsize', 8)
        plot(plotX, plotY)
        axis([datenum(targetYear, 1, 1) datenum(targetYear, 12, 31) ...
                yRanges{indexVar}])
        set(gca, 'xtick', tickLoc, 'xticklabel', theLabels)
        ylabel(yLabels{indexVar})
        title(plotTitles{indexVar})
        box('on')
    end
    
    plotLoc2 = repmat([10 13 16], 1, length(theXs));
    
    subplot(5,1,4), hold on
    set(gca, 'fontsize', 8)
    for index = 1:length(allCol)
        p = patch(theXs(:,index), theYs(:,index), allCol{index});
    end
    for index = 1:length(allCol)
        if strcmp(allCol{index}, 'r') & (max(theXs(:,index)) > datenum(targetYear, 1, 1)) & ...
                (min(theXs(:,index)) < datenum(targetYear, 12, 31))
            t = text(mean(theXs(:,index))+10, plotLoc2(index), ...
                num2str(max(theYs(:,index)), '%.2f'));
            q = quiver(mean(theXs(:,index))+10, plotLoc2(index), -10, 0);
            set(q, 'color', 'k', 'clipping', 'off')
            set(t, 'color', 'k', 'fontsize', 8, 'clipping', 'off');
        end
    end
    axis([datenum(targetYear, 1, 1) datenum(targetYear, 12, 31) 0 20])
    set(gca, 'xtick', tickLoc, 'xticklabel', theLabels)
    ylabel('g/m^2/day')
    title('Sediment Trapping Rate (Honjo Trap, 4 mab)')
    box('on')
    
    plotTitle = {'Significant Wave Height, Pressure Standard Deviation, Beam Attenuation,'; ...
            ['Sediment Trapping Rate, and Current Speed:  Site A, ' num2str(targetYear, '%.0f')]};
    suptitle(plotTitle)
    
    printTitle = ['a_mix_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    close
end
