function adcp_scat_plots(theYears, dataDirectory, plotDirectory, siteArchive, theSite)
%
%   adcp_scat_plots(theYears, dataDirectory, plotDirectory, siteArchive, theSite)
%
%   Function to create yearly scatter plots of east vs. north component for
%   hourly averaged and low-passed ADCP measurements for theYears.
%   Output files are named SITE_scat_YYYY.

%Soupy Alexander, 11/6/03

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

%ADCP
timeVar = cell(1,3);
timeVar_lp = cell(1,3);
currVar = cell(1,3);
currVar_lp = cell(1,3);

if strcmp(lower(theSite), 'a')
    depthList = {'5 mbs'; '15 mbs'; 'Deepest Bin'};
else
    depthList = {'5 mbs'; '10 mbs'; 'Deepest Bin'};
end

%Load in the site A 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
    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;
    
    oldTime = timeVar{1};
    timeVar{1} = [oldTime(:); fileTime(:)]; 
    oldTime_lp = timeVar_lp{1};
    timeVar_lp{1} = [oldTime_lp(:); fileTime_lp(:)];
    
    oldCurr = currVar{1};
    currVar{1} = [oldCurr(:); fileCurr(:)]; 
    oldCurr_lp = currVar_lp{1};
    currVar_lp{1} = [oldCurr_lp(:); fileCurr_lp(:)];
    
    
    if strcmp(lower(theSite), 'a')
        theHit = value2Index(fileDepth, 15);
    elseif strcmp(lower(theSite), 'b')
        theHit = value2Index(fileDepth, 10);
    end
    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;
    
    oldTime = timeVar{2};
    timeVar{2} = [oldTime(:); fileTime(:)]; 
    oldTime_lp = timeVar_lp{2};
    timeVar_lp{2} = [oldTime_lp(:); fileTime_lp(:)];
    
    oldCurr = currVar{2};
    currVar{2} = [oldCurr(:); fileCurr(:)]; 
    oldCurr_lp = currVar_lp{2};
    currVar_lp{2} = [oldCurr_lp(:); fileCurr_lp(:)];
    
    
    theHit = find(fileDepth == max(fileDepth));
        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;
    
    oldTime = timeVar{3};
    timeVar{3} = [oldTime(:); fileTime(:)]; 
    oldTime_lp = timeVar_lp{3};
    timeVar_lp{3} = [oldTime_lp(:); fileTime_lp(:)];
    
    oldCurr = currVar{3};
    currVar{3} = [oldCurr(:); fileCurr(:)]; 
    oldCurr_lp = currVar_lp{3};
    currVar_lp{3} = [oldCurr_lp(:); fileCurr_lp(:)];
end

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    
    f = figure; orient tall
    
    for indexVar = 1:3
        
        theTime = timeVar{indexVar};
        theU = real(currVar{indexVar});
        theV = imag(currVar{indexVar});
        
        gregTime = gregorian(theTime);
        inYear = find(gregTime(:,1) == targetYear);
        
        theTime = theTime(inYear);
        theU = theU(inYear);
        theV = theV(inYear);
        
        bads = find(isnan(theU + theV));
        theTime(bads) = [];
        theU(bads) = [];
        theV(bads) = [];
        
        subplot(3,2,2*indexVar-1)
        
        if ~isempty(theU)
            [majAx, majAz, minAx, minAz, theEllip] = pcaben(theU, theV);
            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(theU);
            meanV = mean(theV);
            
            plot(theU, theV, '.')
            hold on
            plot(xx2 + meanU, yy2 + meanV, 'r', 'linewidth', 2);
            l = line([-50 50], [0 0]);
            l2 = line([0 0], [-50 50]);
            set(l, 'color', 'k', 'linestyle','-.');
            set(l2, 'color', 'k', 'linestyle','-.');
        else
            close
            clear f
            continue
        end
        
        axis equal
        axis tight
        axis([-50 50 -50 50])
        set(gca, 'xtick', [-50:10:50]);
        set(gca, 'ytick', [-50:10:50]);
        xlabel('Northward Component (cm/s)')
        ylabel('Eastward Component (cm/s)')
        box('on')
        
        title(['Hourly Averaged - ' depthList{indexVar}])
        
        subplot(3,2,2*indexVar)
        
        if ~isempty(theU)
            theTime = timeVar_lp{indexVar};
            theU = real(currVar_lp{indexVar});
            theV = imag(currVar_lp{indexVar});
            
            gregTime = gregorian(theTime);
            inYear = find(gregTime(:,1) == targetYear);
            
            theTime = theTime(inYear);
            theU = theU(inYear);
            theV = theV(inYear);
            
            bads = find(isnan(theU + theV));
            theTime(bads) = [];
            theU(bads) = [];
            theV(bads) = [];
            
            [majAx, majAz, minAx, minAz, theEllip] = pcaben(theU, theV);
            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(theU);
            meanV = mean(theV);
            
            plot(theU, theV, '.')
            hold on
            plot(xx2 + meanU, yy2 + meanV, 'r', 'linewidth', 2);
            
            l = line([-50 50], [0 0]);
            l2 = line([0 0], [-50 50]);
            set(l, 'color', 'k', 'linestyle','-.');
            set(l2, 'color', 'k', 'linestyle','-.');
        else
            t = text(0, 0, 'No instrumentation');
            set(t, 'horiz', 'center')
        end
        
        axis equal
        axis tight
        axis([-50 50 -50 50])
        labels = [-50 -40 -30 -20 -10 0 10 20 30 40 50];
        set(gca, 'xtick', [-50:10:50]);
        set(gca, 'ytick', [-50:10:50]);
        xlabel('Northward Component (cm/s)')
        ylabel('Eastward Component (cm/s)')
        box('on')
        
        title(['Low-Passed - ' depthList{indexVar}])
  
    end
    
    if ~exist('f')
        continue
    end
    
    outTitle = {['Site ' upper(theSite) ' - ' num2str(targetYear, '%.0f') ...
                ': Current and Variability']; 
        'ADCP'};
    suptitle(outTitle)
    
    printTitle = [lower(theSite) '_scat_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    
    close
    clear f
end