function compare_sites(theYears, dataDirectory, plotDirectory, siteA, siteB)
%
%   compare_sites(theYears, dataDirectory, plotDirectory, siteA, siteB)
%
%   Function to create comparison plots of the wind stress, current at
%   sites A and B, and STD of pressure at site A.

%Soupy Alexander, 11/6/03

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

colList = {'r'; 'g'; 'b'; 'c'; 'm'; 'k'; 'r'; 'g'; 'b'; 'c'; 'k'; 'm'};
monthList = {'Jan'; 'Feb'; 'Mar'; 'Apr'; 'May'; 'Jun'; 'Jul'; 'Aug'; 'Sep'; ...
        'Oct'; 'Nov'; 'Dec'};
theSpacing = 15;
theLimit = 15;

theSites = {'A'; 'B'};
for indexSite = 1:2
    eval(['siteArchive = site' upper(theSites{indexSite}) ';']);
    
    timeVar_lp = [];
    currVar_lp = [];
    moorList = cell(0);

    %Load in the data
    for indexFile = 1:length(siteArchive)
        if rem(indexFile,10) == 0
            fprintf(['On file number ' num2str(indexFile, '%.0f') ' of ' ...
                    num2str(length(siteArchive), '%.0f') '.\n']);
        end
        fileName = siteArchive(indexFile).fileList;
        if isempty(strfind(fileName, 'adc'))
            continue
        end
        moorNum = fileName(1:4);
        if strmatch(moorNum, moorList)
            continue
        end
        currFileA = netcdf(fullfile(dataDirectory, siteArchive(indexFile).fileList), 'nowrite');
        fileDepth = currFileA{'depth'}(:);
        fileTime = singlejd(currFileA{'time'}(:), currFileA{'time2'}(:));
        fileU = currFileA{'u_1205'}(:);
        fileV = currFileA{'v_1206'}(:);
        bads = find(fileU > 1e3 | fileV > 1e3);
        fileU(bads) = NaN;
        fileV(bads) = NaN;
        
        theHit = value2Index(fileDepth, 5);
        theU = fileU(:, theHit);
        try
            [theU_lp, fileTime_lp] = plfilt(theU, fileTime, 6); 
        catch
            continue
        end
        theV = fileV(:, theHit);
        [theV_lp, fileTime_lp] = plfilt(theV, fileTime, 6); 
        
        fileCurr = theU + i.*theV;
        fileCurr_lp = theU_lp + i.*theV_lp;
        
        timeVar_lp = [timeVar_lp(:); fileTime_lp(:)];
        currVar_lp = [currVar_lp(:); fileCurr_lp(:)];
        
    end
    eval(['timeVar_lp' upper(theSites{indexSite}) '= timeVar_lp;']);
    eval(['currVar_lp' upper(theSites{indexSite}) '= currVar_lp;']);
end

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

timeSDP = [];
SDP = [];
moorList = cell(0);

%Load in the data
for indexFile = 1:length(siteA)
    if rem(indexFile,10) == 0
        fprintf(['On file number ' num2str(indexFile, '%.0f') ' of ' ...
                num2str(length(siteA), '%.0f') '.\n']);
    end
    fileName = siteA(indexFile).fileList;
    fileVars = siteA(indexFile).theVars;
    if isempty(strmatch('P_4023', fileVars))
        continue
    end
    moorNum = fileName(1:4);
    if strmatch(moorNum, moorList)
        continue
    end
    moorList = [moorList(:); {moorNum}];
    theFileA = netcdf(fullfile(dataDirectory, siteA(indexFile).fileList), 'nowrite');
    fileTime = singlejd(theFileA{'time'}(:), theFileA{'time2'}(:));
    fileData = theFileA{'P_4023'}(:);
    fileData(find(fileData > 1e6 | fileData < 1e-6)) = NaN;
    fileData = fileData - nanmean(fileData);
    try
        [fileData_lp, fileTime_lp] = plfilt(fileData, fileTime, 6);
    catch
        continue
    end
    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(:)];
end


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

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    
    f = figure; orient landscape
    
    theTime = windTime;
    theU = tauU;
    theV = tauV;
    
    bads = find(isnan(theU + theV));
    theTime(bads) = [];
    theU(bads) = [];
    theV(bads) = [];
    
    gregTime = gregorian(theTime);
    
    if isempty(find(gregTime(:,1) == targetYear))
        close
        continue
    end
    
    subplot(4,1,1), hold on
    
    for indexMonth = 1:12
        inYear = find(gregTime(:,1) == targetYear & ...
            gregTime(:,2) == indexMonth);
        
        theInTime = theTime(inYear);
        theInU = theU(inYear);
        theInV = theV(inYear);
        
        bads = find(isnan(theInU + theInV));
        theInTime(bads) = [];
        theInU(bads) = [];
        theInV(bads) = [];
        
        
        if length(theInU) > 0
            [majAx, majAz, minAx, minAz, theEllip] = 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);
            
            q = quiver(theSpacing/5*indexMonth, 0, meanU, meanV, colList{indexMonth});
            hold on
            plot(theSpacing/5*indexMonth+xx2 + meanU, yy2 + meanV, colList{indexMonth});
            
            l = line([-10 theSpacing/5*15], [0 0]);
            set(l, 'color', 'k')
        end  
    end
    
    axis equal
    axis([0 theSpacing/5*13 -3 3])
    set(gca, 'xtick', [theSpacing/5:theSpacing/5:12*theSpacing/5]);
    set(gca, 'xticklabel', monthList)
    set(gca, 'ytick', [-3:1:3]);
    ylabel('dynes/cm^2')
    box('on')
    
    title('Wind Stress')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    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);
    
    if isempty(find(gregTime(:,1) == targetYear))
        close
        continue
    end
    
    subplot(4,1,2), hold on
    
    for indexMonth = 1:12
        inYear = find(gregTime(:,1) == targetYear & ...
            gregTime(:,2) == indexMonth);
        
        theInTime = theTime(inYear);
        theInU = theU(inYear);
        theInV = theV(inYear);
        
        bads = find(isnan(theInU + theInV));
        theInTime(bads) = [];
        theInU(bads) = [];
        theInV(bads) = [];
        
        
        if length(theInU) > 0
            [majAx, majAz, minAx, minAz, theEllip] = 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);
            
            q = quiver(theSpacing*indexMonth, 0, meanU, meanV, colList{indexMonth});
            hold on
            plot(theSpacing*indexMonth+xx2 + meanU, yy2 + meanV, colList{indexMonth});
            
            l = line([-10 theSpacing*15], [0 0]);
            set(l, 'color', 'k')
        end  
    end
    
    axis equal
    axis([0 theSpacing*13 -15 15])
    set(gca, 'xtick', [theSpacing:theSpacing:12*theSpacing]);
    set(gca, 'xticklabel', monthList)
    set(gca, 'ytick', [-15:5:15]);
    ylabel('cm/s')
    box('on')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    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);
    
    if isempty(find(gregTime(:,1) == targetYear))
        close
        continue
    end
    
    subplot(4,1,3), hold on
    
    for indexMonth = 1:12
        inYear = find(gregTime(:,1) == targetYear & ...
            gregTime(:,2) == indexMonth);
        
        theInTime = theTime(inYear);
        theInU = theU(inYear);
        theInV = theV(inYear);
        
        bads = find(isnan(theInU + theInV));
        theInTime(bads) = [];
        theInU(bads) = [];
        theInV(bads) = [];
        
        
        if length(theInU) > 0
            [majAx, majAz, minAx, minAz, theEllip] = 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);
            
            q = quiver(theSpacing*indexMonth, 0, meanU, meanV, colList{indexMonth});
            hold on
            plot(theSpacing*indexMonth+xx2 + meanU, yy2 + meanV, colList{indexMonth});
            
            l = line([-10 theSpacing*15], [0 0]);
            set(l, 'color', 'k')
        end  
    end
    
    axis equal
    axis([0 theSpacing*13 -15 15])
    set(gca, 'xtick', [theSpacing:theSpacing:12*theSpacing]);
    set(gca, 'xticklabel', monthList)
    set(gca, 'ytick', [-15:5:15]);
    ylabel('cm/s')
    box('on')
    
    
    title('Current, Site B, 5 mbs')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%
    %Standard deviation of pressure
    theTime = timeSDP;
    theData = SDP;
    
    gregTime = gregorian(theTime);
    
    subplot(4,1,4), hold on
    
    for indexMonth = 1:12
        inYear = find(gregTime(:,1) == targetYear & ...
            gregTime(:,2) == indexMonth);
        
        theInTime = theTime(inYear);
        theInData = theData(inYear);
        
        hitData(indexMonth) = nanstd(theInData);
    end
    
    theXs = [theSpacing*2/3:theSpacing*2/3:12*theSpacing*2/3];
    b = 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', monthList)
    set(gca, 'ytick', [0:5:20]);
    ylabel('mb')
    box('on')
    title('Standard Deviation of Low-Passed Pressure')
    
    theTitle = {'Monthly Mean and Variability of Low-Passed Wind Stress, ADCP Current (5 mbs) from Sites A & B'; ...
        ['and St. Deviation of Pressure from Site A, ' num2str(targetYear, '%.0f')]};
    suptitle(theTitle)
    printTitle = ['compare_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    close
end