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

%Soupy Alexander, 11/6/03

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

%PCM
theFiles = {'a_curr_5mbs.nc'; 'a_curr_10mab.nc'; 'a_curr_1mab.nc'};
timeVar = cell(1,3);
currVar = cell(1,3);
depthList = {'5 mbs'; '10 mab'; '1 mab'};

for indexFile = 1:length(theFiles)
    ncID = netcdf(fullfile(dataDirectory, theFiles{indexFile}), 'nowrite');
    timeVar{indexFile} = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    fileU = ncID{'u_1205'}(:);
    fileV = ncID{'v_1206'}(:);
    currVar{indexFile} = fileU + i.*fileV;
end

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    
    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
            t = text(0, 0, 'No instrumentation');
            set(t, 'horiz', 'center')
        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)
            timeNew = min(theTime):1/24:max(theTime);
            newU = smart_interp(theTime, theU, timeNew, 6);
            newV = smart_interp(theTime, theV, timeNew, 6);
            
            [theU] = plfilt(newU, timeNew, 6);
            [theV] = plfilt(newV, timeNew, 6);
            
            bads = find(isnan(theU + theV));
            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
    
    outTitle = {['Site A - ' num2str(targetYear, '%.0f') ': Current and Variability']; 
        'Point Current Meter Data'};
    suptitle(outTitle)
    
    printTitle = ['pcm_scat_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    
    close
end