function mixed_plot_b(theYears, dataDirectory, plotDirectory, concatDirectory)
%
%   mixed_plot_b(theYears, dataDirectory, plotDirectory, concatDirectory)
%
%   Creates the plots of significant wave height,
%   sediment trapping rate, and current speed (1 mab) at Site B.  Plots are
%   named b_mix_YYYY.

%This m-file written by Soupy Alexander (11/2003) for U.S. Geological 
%Survey Data Series 74, Version 2.0. Updated by Soupy Dalyander (8/2008) 
%for U.S. Geological Survey Data Series 74, Version 3.0 (citation below).
%Soupy Dalyander previously Soupy Alexander.  

%This Matlab m-file was used to create portions of U.S. Geological Survey 
%Data Series 74, Version 3.0. Although this program has been used by the 
%USGS, no warranty, expressed or implied, is made by the USGS or the 
%United States Government as to the accuracy and functioning of the 
%program and related program material nor shall the fact of distribution 
%constitute any such warranty, and no responsibility is assumed by the 
%USGS in connection therewith.

%Citation:  Butman, Bradford, Dalyander, P.S., Bothner, M.H., Borden,
%Jonathan, Casso, M.A., Gutierrez, B.T., Hastings, M.E., Lightsom, F.L., 
%Martini, M.A., Montgomery, E.T, Rendigs, R.R., and Strahle, W.S., 2008, 
%Long-Term Oceanographic Observations in Massachusetts Bay, 1989—2006: 
%U.S. Geological Survey Data Series 74, Version 3.0, DVD-ROM. 
%Also online at http://pubs.usgs.gov/ds/74/.

%Wave data
waveFile = netcdf(fullfile(dataDirectory, '44013-84t06-cal.nc'), 'nowrite');
waveTime = singlejd(waveFile{'time'}(:), waveFile{'time2'}(:));
waveData = waveFile{'wh_4061'}(:);
waveData(waveData > 1e3) = NaN;
[waveTime, waveData] = nanfill_2(waveTime, waveData);
%Current Speed
currTime = [];
currData = [];

[theNums, theText] = xlsread...
    (fullfile(concatDirectory, 'adcpB_c'));

toUse = theNums(:,1);
fileList = theText(2:end,1);

%Load in the site A data
for indexFile = 1:length(toUse)
    if isnan(toUse(indexFile))
        continue
    end

    currFileA = netcdf(fullfile([dataDirectory '\MBLT_B'], fileList{indexFile}), 'nowrite');
    fileDepth = currFileA{'depth'}(:);
    fileTime = singlejd(currFileA{'time'}(:), currFileA{'time2'}(:));
    fileU = currFileA{'u_1205'}(:)/100;
    fileV = currFileA{'v_1206'}(:)/100;
    bads = find(fileU > 1e3 | fileV > 1e3);
    fileU(bads) = NaN;
    fileV(bads) = NaN;

    %theHit = find(fileDepth == max(fileDepth));
    theHit = value2index(fileDepth, 17);
    fileU = fileU(:, theHit);
    fileV = fileV(:, theHit);
    fileData = sqrt(fileU.^2 + fileV.^2);

    currTime = [currTime(:); fileTime(:)];
    currData = [currData(:); fileData(:)];

    close(currFileA)
    clear currFileA
end

newTime = min(currTime):1/24:max(currTime);
[newData] = smart_interp(currTime, currData, newTime, 5);
currTime = newTime;
currData = newData;


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

honjoNums = theNums(theHonjos,:);
honjoText = theText(theHonjos,:);

honjoStartText = honjoText(:,3);
for index = 1:length(honjoStartText)
    try
        honjoStart(index) = datenum(honjoStartText{index});
    catch
        honjoStart(index) = NaN;
    end
end

honjoEndText = honjoText(:,4);
for index = 1:length(honjoEndText)
    try
        honjoEnd(index) = datenum(honjoEndText{index});
    catch
        honjoEnd(index) = NaN;
    end
end
honjoFlux = honjoNums(:, 11);

bads = find(isnan(honjoFlux(:)) | isnan(honjoStart(:)) | isnan(honjoEnd(:)));
honjoStart(bads) = [];
honjoEnd(bads) = [];
honjoFlux(bads) = [];

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

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'; 'curr'};
plotTitles = {'\it\bfA    \rm SIGNIFICANT WAVE HEIGHT'; ...
    '\it\bfC    \rm CURRENT SPEED AT ~17 M'};
plotNums = [1 3];
yRanges = {[0 10]; [0 0.40]};
yLabels = {'HEIGHT, IN M'; ...
    {'SPEED, IN M/S'}};
yTicks = {0:5:10; 0:0.10:0.40};

for indexYear = 1:length(theYears)
    targetYear = theYears(indexYear);
    
    if (targetYear < 1997) || (targetYear > 2004)
        continue
    end
    
    tickLoc = [];
    for indexMonth = 1:12
        tickLoc = [tickLoc; datenum(targetYear, indexMonth, 1)];
    end
    
    figure, orient landscape

    indexVar = 1;
    eval(['theTime = julian2datenum(' plotVars{indexVar} 'Time);']);
    eval(['theData = ' plotVars{indexVar} 'Data;']);
    usgs_plot(targetYear, 3, plotNums(indexVar), plotTitles{indexVar}, ...
        yLabels{indexVar}, yRanges{indexVar}, ...
        yTicks{indexVar});
    [plotX, plotY] = nanfill_2(theTime, theData);
    plot(plotX, plotY)

    plotLoc2 = repmat([10 13 16], 1, length(theXs));
    usgs_plot(targetYear, 3, 2, '\it\bfB    \rm SEDIMENT-TRAPPING RATE, TIME-SERIES SEDIMENT TRAP, ~4 M ABOVE BOTTOM', ...
        {'TRAPPING RATE,'; 'IN G/M^ ^2/DAY'}, [0 20], ...
        0:5:20)
    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)), 18, ...
                num2str(max(theYs(:,index)), '%.0f'));
            set(t, 'fontname', 'arial', 'fontsize', 6, 'color', 'k', ...
                'clipping', 'off', 'horiz', 'center', 'vert', 'middle')
            tPos = get(t, 'position');
            if tPos(1) < min(xlim)
                tPos(1) = min(xlim);
                set(t, 'horiz', 'left');
            elseif tPos(1) > max(xlim)
                tPos(1) = max(xlim);
                set(t, 'horiz', 'right');
            end
        end
    end

    indexVar = 2;
    eval(['theTime = julian2datenum(' plotVars{indexVar} 'Time);']);
    eval(['theData = ' plotVars{indexVar} 'Data;']);
    usgs_plot(targetYear, 3, plotNums(indexVar), plotTitles{indexVar}, ...
        yLabels{indexVar}, yRanges{indexVar}, ...
        yTicks{indexVar});
    [plotX, plotY] = nanfill_2(theTime, theData);
    plot(plotX, plotY)
    if targetYear == 2004
        t = text(datenum(targetYear, 7, 1), mean(ylim), 'No instrumentation.');
        set(t, 'horizontal', 'center', 'fontname', 'arial', 'fontsize', 8)
    end


    plotTitle = {'Significant Wave Height at NOAA Buoy 44013, and Pressure Standard Deviation, Beam Attenuation,'; ...
            'Sediment-Trapping Rate, and Current Speed at LT-B'; num2str(targetYear, '%.0f')};
    s = suptitle(plotTitle);
    set(s, 'fontname', 'arial', 'fontsize', 12, 'position', [0.5 -0.08 0]);
    
    printTitle = ['b_mix_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    close
end
