function pcm_ell_plots(theYears, dataDirectory, plotDirectory)
%
%   pcm_ell_plots(theYears, dataDirectory, plotDirectory)
%
%   Function to create ellipse plots of east and north components by month 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];

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 = 10;
theLimit = 15;

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

for indexFile = 1:length(theFiles)
    ncID = netcdf(fullfile(dataDirectory, theFiles{indexFile}), 'nowrite');
    fileTime = singlejd(ncID{'time'}(:), ncID{'time2'}(:));
    timeVar{indexFile} = fileTime;
    fileU = ncID{'u_1205'}(:);
    fileV = ncID{'v_1206'}(:);
    currVar{indexFile} = fileU + i.*fileV;
    
    timeNew = min(fileTime):1/24:max(fileTime);
    newU = smart_interp(fileTime, fileU, timeNew, 6);
    newV = smart_interp(fileTime, fileV, timeNew, 6);
    
    [theU, time_lp] = plfilt(newU, timeNew, 6);
    [theV, time_lp] = plfilt(newV, timeNew, 6);
    
    timeVar_lp{indexFile} = time_lp;
    currVar_lp{indexFile} = theU + i.*theV;
    
end

for indexYear = 1:length(theYears)
    
    targetYear = theYears(indexYear);
    
    figure, orient landscape
    
    for indexVar = 1:3
        
        theTime = timeVar{indexVar};
        theU = real(currVar{indexVar});
        theV = imag(currVar{indexVar});
        
        gregTime = gregorian(theTime);
        
        subplot(3,1,indexVar), 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) > 10*24
                [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
        
        if isempty(find(gregTime(:,1) == targetYear))
            t = text(theSpacing*6.5, 0, 'No Instrumentation');
            set(t, 'horiz', 'center')
        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([depthList{indexVar}])
    end
    
    outTitle = {['Site A - ' num2str(targetYear, '%.0f') ': Monthly Mean Current and Variability']; 
        'Point Current Meter Hour-Averaged Data'};
    suptitle(outTitle)
    
    printTitle = ['pcm_ell_h_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    
    close
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    figure, orient landscape
    
    for indexVar = 1:3
        
        theTime = timeVar_lp{indexVar};
        theU = real(currVar_lp{indexVar});
        theV = imag(currVar_lp{indexVar});
        
        gregTime = gregorian(theTime);
        
        subplot(3,1,indexVar), 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) > (10*4)
                [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
        
        if isempty(find(gregTime(:,1) == targetYear))
            t = text(theSpacing*6.5, 0, 'No Instrumentation');
            set(t, 'horiz', 'center')
        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([depthList{indexVar}])
    end
    
    outTitle = {['Site A - ' num2str(targetYear, '%.0f') ': Monthly Mean Current and Variability']; 
        'Point Current Meter Low-Passed Data'};
    suptitle(outTitle)
    
    printTitle = ['pcm_ell_lp_' num2str(targetYear, '%.0f')];
    print('-dpdf', fullfile(plotDirectory, printTitle))
    
    close
end