function compare_sites(theYears, dataDirectory, plotDirectory, concatFile, ...
    concatDirectory)
%
%   compare_sites(theYears, dataDirectory, plotDirectory, concatFile
%       concatDirectory)
%
%   Function to create comparison plots of the wind stress, current at
%   sites A and B, and STD of pressure at LT-A.

%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/.

theta=2*pi*(1:64)/64;
theta=[0 theta];

xLabels = {'JAN.'; 'FEB.'; 'MAR.'; 'APR.'; 'MAY'; 'JUNE'; 'JULY'; 'AUG.'; ...
    'SEPT.'; 'OCT.'; 'NOV.'; 'DEC.'};
axShift = [0 0.026 0.056 0.086];
titLocs = [0.785 0.606 0.432 0.257];

colList = {'r'; 'g'; 'b'; 'c'; 'm'; 'k'; 'r'; 'g'; 'b'; 'c'; 'k'; 'm'};
theSpacing = 15;

n = 0.18;

theSites = {'A'; 'B'};

for indexSite = 1:2
    theSite = theSites{indexSite};
    timeVar_lp = [];
    currVar_lp = [];

    [theNums, theText] = xlsread...
        (fullfile(concatDirectory, ['adcp' upper(theSite) '_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 '\LT_' upper(theSite)], 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 > 1e2 | fileV > 1e2);
        fileU(bads) = NaN;
        fileV(bads) = NaN;

        hit5 = value2Index(fileDepth, 5);
        cU5 = fileU(:, hit5);
        try
            [fcU5_lp] = plfilt(cU5, fileTime, 6);
        catch
            %File too short to LP
            continue
        end
        cV5 = fileV(:, hit5);
        [fcV5_lp, fileTime_lp] = plfilt(cV5, fileTime, 6);
        fileCurr_lp = fcU5_lp + i.*fcV5_lp;

        timeVar_lp = [timeVar_lp(:); fileTime_lp(:)];
        currVar_lp = [currVar_lp(:); fileCurr_lp(:)];

        close(currFileA)
        clear currFileA
    end
    eval(['timeVar_lp' upper(theSites{indexSite}) '= timeVar_lp;']);
    eval(['currVar_lp' upper(theSites{indexSite}) '= currVar_lp;']);
end

%Load in the PCM observations from LT-A
fileData = load(concatFile, ...
    'A_current_S_time_LP', 'A_currentU_S_data_LP', 'A_currentV_S_data_LP');
cTime1 = fileData.A_current_S_time_LP;
cU = fileData.A_currentU_S_data_LP/100;
cV = fileData.A_currentV_S_data_LP/100;
pcm = cU + i*cV;

crossPoint = value2index(cTime1, julian(1997, 01, 19, 00));
timeVar_lpA = [cTime1(1:crossPoint); timeVar_lpA];
currVar_lpA = [pcm(1:crossPoint); currVar_lpA];


%%%%%%%%%%%%%%%%%%%%%%%%%

%Standard deviation of pressure
[nums, theText] = xlsread...
    (fullfile(concatDirectory, 'pressureA_c.xls', '1-2 mab'));

text2 = theText(2:end,:);
fileList = text2(nums(:,1)==1,1);
timeSDP = [];
SDP = [];
for index4 = 1:length(fileList)
    thisFile = fileList{index4};
    fileVar = 'P_4023';
    
    ncID = netcdf([dataDirectory '\LT_A\' thisFile], 'nowrite');
    fileTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    eval(['fileData = ncID{ ''' fileVar '''}(:);']);
    if isempty(fileData)
        fileData = ncID{'P_4022'}(:);
        fileVar = 'P_4022';
    end
    if isempty(fileData)
        fileData = ncID{'P_1'}(:);
        fileVar = 'P_1';
    end

    if strcmp(thisFile, '7173advs-calq.cdf');
        stopTime = value2index(fileTime, julian(2004, 02, 05, 00));
        fileTime(stopTime:end) = [];
        fileData(stopTime:end) = [];
    end
    
    fileData(fileData > 1e9) = NaN;
    fileData(end) = NaN;
    fileData = fileData - nanmean(fileData);
    if strcmp(fileVar, 'P_1')
        fileData = fileData*100;
    end
    [fileData_lp, fileTime_lp] = plfilt(fileData, fileTime, 6);


    bads = find(abs(diff(fileData_lp)) > 25);
    if length(bads) > 2
        fileData_lp = [(fileData_lp(1:(min(bads)-1)) - ...
            nanmean(fileData_lp(1:(min(bads)-1)))); ...
            (fileData_lp(max(bads)+1:end) - ...
            nanmean(fileData_lp(max(bads)+1:end)))];
        fileTime_lp = [fileTime_lp(1:(min(bads)-1)); ...
            fileTime_lp(max(bads)+1:end)];
    end

    timeSDP = [timeSDP(:); fileTime_lp(:)];
    SDP = [SDP(:); fileData_lp(:)];
    close(ncID)
    clear ncID
end

%Load in the wind stress
[windTime, tauU, tauV] = mblt_windstress(dataDirectory);
tauU = tauU/10; %Convert to Pa
tauV = tauV/10;

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    figure; orient landscape
    
    theTime = windTime;
    theU = tauU;
    theV = tauV;
    
    bads = find(isnan(theU + theV));
    theTime(bads) = [];
    theU(bads) = [];
    theV(bads) = [];
    
    gregTime = gregorian(theTime);

    subplot(4,1,1), hold on
    set(gca,'fontname', 'arial', 'fontsize', 8);
    thePos = get(gca, 'position');
    
    if ~isempty(find(gregTime(:,1) == targetYear, 1))

        for indexMonth = 1:12
            inYear = find(gregTime(:,1) == targetYear & ...
                gregTime(:,2) == indexMonth);

            theInU = theU(inYear);
            theInV = theV(inYear);

            bads = find(isnan(theInU + theInV));
            theInU(bads) = [];
            theInV(bads) = [];


            if ~isempty(theInU)
                [majAx, majAz, minAx] = pcaben(theInU, theInV);
                x2=(majAx*cos(theta))';
                y2=(minAx*sin(theta))';
                theAngle=-1*majAz*pi/180+pi/2;
                xx2 = x2(:).*cos(theAngle) - y2.*sin(theAngle);
                yy2 = x2(:).*sin(theAngle) + y2.*cos(theAngle);

                meanU = mean(theInU);
                meanV = mean(theInV);

                quiver(theSpacing/50*indexMonth, 0, meanU, meanV, colList{indexMonth});
                hold on
                plot(theSpacing/50*indexMonth+xx2 + meanU, yy2 + meanV, colList{indexMonth});
            end
        end
        l = line([-10 theSpacing/50*15], [0 0]);
        set(l, 'color', 'k')
    else
        t = text(theSpacing/100*13, 0, 'No instrumentation.');
        set(t, 'fontname', 'arial', 'fontsize', 8, 'horiz', 'center')
    end

    axis equal
    axis([0 theSpacing/50*13 -0.3 0.3])
    set(gca, 'xtick', theSpacing/50:theSpacing/50:12*theSpacing/50);
    set(gca, 'xticklabel', [])
    set(gca, 'ytick', -0.3:0.1:0.3);
    ylabel('STRESS, IN PA')
    annotation('textbox', [n titLocs(1) 0.1 0.1], 'string', '\it\bfA    \rm WIND STRESS AT NOAA BUOY 44013', ...
        'linestyle', 'none', 'fontname', 'arial', 'fontsize', 8);
    thePos(2) = thePos(2) + axShift(1);
    set(gca, 'position', thePos)
    add_right_ticks(gca)


    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    theTime = timeVar_lpA;
    theU = real(currVar_lpA);
    theV = imag(currVar_lpA);
    
    bads = find(isnan(theU + theV));
    theTime(bads) = [];
    theU(bads) = [];
    theV(bads) = [];
    
    gregTime = gregorian(theTime);
    
    subplot(4,1,2), hold on
    set(gca, 'fontname', 'arial', 'fontsize', 8)
    thePos = get(gca, 'position');

    if ~isempty(find(gregTime(:,1) == targetYear, 1))

        for indexMonth = 1:12
            inYear = find(gregTime(:,1) == targetYear & ...
                gregTime(:,2) == indexMonth);

            theInU = theU(inYear);
            theInV = theV(inYear);

            bads = find(isnan(theInU + theInV));
            theInU(bads) = [];
            theInV(bads) = [];

            
            if ~isempty(theInU)
                [majAx, majAz, minAx] = pcaben(theInU, theInV);
                x2=(majAx*cos(theta))';
                y2=(minAx*sin(theta))';
                theAngle=-1*majAz*pi/180+pi/2;
                xx2 = x2(:).*cos(theAngle) - y2.*sin(theAngle);
                yy2 = x2(:).*sin(theAngle) + y2.*cos(theAngle);

                meanU = mean(theInU);
                meanV = mean(theInV);

                quiver(theSpacing/100*indexMonth, 0, meanU, meanV, colList{indexMonth});
                hold on
                plot(theSpacing/100*indexMonth+xx2 + meanU, yy2 + meanV, colList{indexMonth});

                l = line([-10 theSpacing*15], [0 0]);
                set(l, 'color', 'k')
            end
        end
    else
        t = text(15/200*13, 0, 'No instrumentation.');
        set(t, 'fontname', 'arial', 'fontsize', 8, 'horiz', 'center')
    end
    axis equal
    axis([0 theSpacing/100*13 -0.15 0.15])
    set(gca, 'xtick', theSpacing/100:theSpacing/100:12*theSpacing/100);
    set(gca, 'xticklabel', [])
    set(gca, 'ytick', -0.15:0.05:0.15);
    ylabel('CURRENT, IN M/S')
    annotation('textbox', [n titLocs(2) 0.1 0.1], 'string', '\it\bfB    \rm NEAR-SURFACE CURRENT (~5 M) AT LT-A', ...
        'linestyle', 'none', 'fontname', 'arial', 'fontsize', 8);
    thePos(2) = thePos(2) + axShift(2);
    set(gca, 'position', thePos)
    add_right_ticks(gca)
    
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    theTime = timeVar_lpB;
    theU = real(currVar_lpB);
    theV = imag(currVar_lpB);
    
    bads = find(isnan(theU + theV));
    theTime(bads) = [];
    theU(bads) = [];
    theV(bads) = [];

    gregTime = gregorian(theTime);

    subplot(4,1,3), hold on
    set(gca, 'fontsize', 8, 'fontname', 'arial')
    thePos = get(gca, 'position');

    if ~isempty(find(gregTime(:,1) == targetYear, 1))

        for indexMonth = 1:12
            inYear = find(gregTime(:,1) == targetYear & ...
                gregTime(:,2) == indexMonth);

            theInU = theU(inYear);
            theInV = theV(inYear);

        bads = find(isnan(theInU + theInV));
        theInU(bads) = [];
        theInV(bads) = [];
        
        
        if ~isempty(theInU)
            [majAx, majAz, minAx] = pcaben(theInU, theInV);
            x2=(majAx*cos(theta))';
            y2=(minAx*sin(theta))';
            theAngle=-1*majAz*pi/180+pi/2;
            xx2 = x2(:).*cos(theAngle) - y2.*sin(theAngle);
            yy2 = x2(:).*sin(theAngle) + y2.*cos(theAngle);
            
            meanU = mean(theInU);
            meanV = mean(theInV);

            quiver(theSpacing/100*indexMonth, 0, meanU, meanV, colList{indexMonth});
            hold on
            plot(theSpacing/100*indexMonth+xx2 + meanU, yy2 + meanV, colList{indexMonth});

            l = line([-10 theSpacing*15], [0 0]);
            set(l, 'color', 'k')
        end
        end
    else
        t = text(15/200*13, 0, 'No instrumentation.');
        set(t, 'fontname', 'arial', 'fontsize', 8, 'horiz', 'center')
    end


    axis equal
    axis([0 theSpacing/100*13 -0.15 0.15])
    set(gca, 'xtick', theSpacing/100:theSpacing/100:12*theSpacing/100);
    set(gca, 'xticklabel', [])
    set(gca, 'ytick', -0.15:0.05:0.15);
    ylabel('CURRENT, IN M/S')
    annotation('textbox', [n titLocs(3) 0.1 0.1], 'string', '\it\bfC    \rm NEAR-SURFACE CURRENT (~5 M) AT LT-B', ...
        'linestyle', 'none', 'fontname', 'arial', 'fontsize', 8);
    thePos(2) = thePos(2) + axShift(3);
    set(gca, 'position', thePos)
    add_right_ticks(gca)

    
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    %Standard deviation of pressure
    theTime = timeSDP;
    theData = SDP;
    
    gregTime = gregorian(theTime);
    
    subplot(4,1,4), hold on
    set(gca, 'fontname', 'arial', 'fontsize', 8)
    thePos = get(gca, 'position');
    
    hitData = repmat(NaN, 1, 12);
    for indexMonth = 1:12
        
        theInData = theData(gregTime(:,1) == targetYear & ...
            gregTime(:,2) == indexMonth);
        
        hitData(indexMonth) = nanstd(theInData);
    end
    
    theXs = theSpacing*2/3:theSpacing*2/3:12*theSpacing*2/3;
    bar(theXs, hitData);
    
    axis equal
    axis([0 theSpacing*2/3*13 0 20])
    set(gca, 'xtick', theSpacing*2/3:theSpacing*2/3:12*theSpacing*2/3);
    set(gca, 'xticklabel', xLabels)
    set(gca, 'ytick', 0:5:20);
    ylabel({'PRESSURE STANDARD'; 'DEVIATION, IN MB'})
    annotation('textbox', [n titLocs(4) 0.1 0.1], 'string', '\it\bfD    \rm STANDARD DEVIATION OF LOW-PASS FILTERED PRESSURE AT LT-A', ...
        'linestyle', 'none', 'fontname', 'arial', 'fontsize', 8);
    thePos(2) = thePos(2) + axShift(4);
    set(gca, 'position', thePos)
    add_right_ticks(gca)

    annotation('textbox', [0.427 0.17 1 0.1], 'string', num2str(targetYear), ...
        'linestyle', 'none', 'fontname', 'arial', 'fontsize', 8, ...
        'horiz', 'center', 'vert', 'bottom');

    theTitle = {'Monthly Mean and Variability of Low-Pass-Filtered Wind Stress at NOAA Buoy 44013 and'; ...
        ['Current (~5 m) at LT-A and LT-B, and Standard Deviation of Low-Pass-Filtered Pressure at LT-A, ' num2str(targetYear, '%.0f')]};
    s = suptitle(theTitle);
    set(s, 'fontsize', 12, 'fontname', 'arial', 'position', [0.5 -0.08 0]);
    printTitle = ['compare_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    close
end
