function sort_currents(theYears, dataDirectory, plotDirectory, siteA, siteB)
%
%   sort_currents(theYears, dataDirectory, plotDirectory, siteA, siteB)
%
%   Function to sort the currents by speed and direction (see the
%   currents_speed_dir.html web page).  Output files are named
%   SITE_DEPTHa_TIMESTEP_YYYY (for ADCP) and SITE_DEPTH_TIMESTEP_YYYY (for point current
%   meters).

%Modified from Ben Guitterez' "binfix", Soupy Alexander, 11/3/03

%Updated by Soupy, 5/3/04

%Set up the files to analyze - PCM
theFiles = {'a_curr_5mbs.nc'; 'a_curr_10mab.nc'; 'a_curr_1mab.nc'};
timeName = {'time5'; 'time10'; 'time1'};
timeName_lp = {'time5_lp'; 'time10_lp'; 'time1_lp'};
currName = {'curr5'; 'curr10'; 'curr1'};
currName_lp = {'curr5_lp'; 'curr10_lp'; 'curr1_lp'};

for index = 1:3
    fID = netcdf(fullfile(dataDirectory, theFiles{index}), 'nowrite');
    fileTime = singlejd(fID{'time'}(:), fID{'time2'}(:));
    fileU = fID{'u_1205'}(:);
    fileU(find(fileU > 1e3)) = NaN;
    fileV = fID{'v_1206'}(:);
    fileV(find(fileV > 1e3)) = NaN;
    fileCurr = fileU + i.*fileV;
    
    [fileU_lp, fileTime_lp] = plfilt(fileU, fileTime, 6);
    [fileV_lp, fileTime_lp] = plfilt(fileV, fileTime, 6);
    fileCurr_lp = fileU_lp + i.*fileV_lp;
    
    eval([timeName{index} ' = fileTime;']);
    eval([currName{index} ' = fileCurr;']);
    eval([timeName_lp{index} ' = fileTime_lp;']);
    eval([currName_lp{index} ' = fileCurr_lp;']);
end

clear fileU fileV fileCurr fID fileTime fileTime_lp index theFiles
clear timeName timeName_lp currName currName_lp

%Set up the files to analyze - Site A, ADCP
aTime = [];
aTime_lp = [];
a_currS = [];
a_currM = [];
a_currD = [];
a_currS_lp = [];
a_currM_lp = [];
a_currD_lp = [];

goodFiles = cell(0,0);

%Load in the site A data
for indexFile = 1:length(siteA)
    fileVars = siteA(indexFile).theVars;
    fileStart = siteA(indexFile).startTime;
    gregStart = gregorian(fileStart);
    fileEnd = siteA(indexFile).endTime;
    gregEnd = gregorian(fileEnd);
    if ~any(strcmp(fileVars, 'Werr_1201')) | ~any(ismember([gregStart(1) gregEnd(1)], theYears))
        continue
    end
    currFileA = netcdf(fullfile(dataDirectory, siteA(indexFile).fileList), 'nowrite');
    fileDepth = currFileA{'depth'}(:);
    fileTime = singlejd(currFileA{'time'}(:), currFileA{'time2'}(:));
    if (length(fileDepth) == 1) | (median(diff(fileTime)) > 2/24) | (median(diff(fileTime)) < 50/(60*24))
        continue
    end
    goodFiles{length(goodFiles)+1} = siteA(indexFile).fileList;
    
    fileU = currFileA{'u_1205'}(:);
    fileV = currFileA{'v_1206'}(:);
    bads = find(fileU > 1e3 | fileV > 1e3);
    fileU(bads) = NaN;
    fileV(bads) = NaN;
    
    hit5 = value2Index(fileDepth, 5);
    cU5 = fileU(:, hit5);
    try
        [fcU5_lp, fileTime_lp] = plfilt(cU5, fileTime, 6); 
    catch
        continue
    end
    cV5 = fileV(:, hit5);
    [fcV5_lp, fileTime_lp] = plfilt(cV5, fileTime, 6);
    
    hit10 = value2Index(fileDepth, 15);
    cU10 = fileU(:, hit10);
    [fcU10_lp, fileTime_lp] = plfilt(cU10, fileTime, 6); 
    cV10 = fileV(:, hit10);
    [fcV10_lp, fileTime_lp] = plfilt(cV10, fileTime, 6);
    
    hit1 = find(fileDepth == max(fileDepth));
    cU1 = fileU(:, hit1);
    [fcU1_lp, fileTime_lp] = plfilt(cU1, fileTime, 6); 
    cV1 = fileV(:, hit1);
    [fcV1_lp, fileTime_lp] = plfilt(cV1, fileTime, 6);
    
    [aTime] = [aTime(:); fileTime(:)];
    [aTime_lp] = [aTime_lp(:); fileTime_lp(:)];
    [a_currS] = [a_currS(:); cU5(:) + i.*cV5(:)];
    [a_currM] = [a_currM(:); cU10(:) + i.*cV10(:)];
    [a_currD] = [a_currD(:); cU1(:) + i.*cV1(:)];
    [a_currS_lp] = [a_currS_lp(:); fcU5_lp(:) + i.*fcV5_lp(:)];
    [a_currM_lp] = [a_currM_lp(:); fcU10_lp(:) + i.*fcV10_lp(:)];
    [a_currD_lp] = [a_currD_lp(:); fcU1_lp(:) + i.*fcV1_lp(:)];
end

%Set up the files to analyze - Site B, ADCP
bTime = [];
bTime_lp = [];
b_currS = [];
b_currM = [];
b_currD = [];
b_currS_lp = [];
b_currM_lp = [];
b_currD_lp = [];

goodFiles = cell(0,0);
%Load in the site B data
for indexFile = 1:length(siteB)
    fileVars = siteB(indexFile).theVars;
    fileStart = siteB(indexFile).startTime;
    gregStart = gregorian(fileStart);
    fileEnd = siteB(indexFile).endTime;
    gregEnd = gregorian(fileEnd);
    if ~any(strcmp(fileVars, 'Werr_1201')) | ~any(ismember([gregStart(1) gregEnd(1)], theYears))
        continue
    end
    currFileA = netcdf(fullfile(dataDirectory, siteB(indexFile).fileList), 'nowrite');
    fileDepth = currFileA{'depth'}(:);
    fileTime = singlejd(currFileA{'time'}(:), currFileA{'time2'}(:));
    thisFile = siteB(indexFile).fileList;
    %Removing non-ADCP files and the NOPP/redundent ADCP's
    if (length(fileDepth) == 1) | ~isempty(findstr('6151adc-a.nc', thisFile)) | ...
            ~isempty(findstr('6281', thisFile)) | ~isempty(findstr('6361', thisFile)) | ...
            ~isempty(findstr('6871', thisFile)) | ~isempty(findstr('6131', thisFile))
        continue
    end
    
    goodFiles{length(goodFiles)+1} = siteB(indexFile).fileList;

    fileU = currFileA{'u_1205'}(:);
    fileV = currFileA{'v_1206'}(:);
    bads = find(fileU > 1e3 | fileV > 1e3);
    fileU(bads) = NaN;
    fileV(bads) = NaN;
    
    hit5 = value2Index(fileDepth, 5);
    cU5 = fileU(:, hit5);
    try
        [fcU5_lp, fileTime_lp] = plfilt(cU5, fileTime, 6); 
    catch
        continue
    end
    cV5 = fileV(:, hit5);
    [fcV5_lp, fileTime_lp] = plfilt(cV5, fileTime, 6);
    
    hit10 = value2Index(fileDepth, 10);
    cU10 = fileU(:, hit10);
    [fcU10_lp, fileTime_lp] = plfilt(cU10, fileTime, 6); 
    cV10 = fileV(:, hit10);
    [fcV10_lp, fileTime_lp] = plfilt(cV10, fileTime, 6);
    
    hit1 = find(fileDepth == max(fileDepth));
    cU1 = fileU(:, hit1);
    [fcU1_lp, fileTime_lp] = plfilt(cU1, fileTime, 6); 
    cV1 = fileV(:, hit1);
    [fcV1_lp, fileTime_lp] = plfilt(cV1, fileTime, 6);
    
    [bTime] = [bTime(:); fileTime(:)];
    [bTime_lp] = [bTime_lp(:); fileTime_lp(:)];
    [b_currS] = [b_currS(:); cU5(:) + i.*cV5(:)];
    [b_currM] = [b_currM(:); cU10(:) + i.*cV10(:)];
    [b_currD] = [b_currD(:); cU1(:) + i.*cV1(:)];
    [b_currS_lp] = [b_currS_lp(:); fcU5_lp(:) + i.*fcV5_lp(:)];
    [b_currM_lp] = [b_currM_lp(:); fcU10_lp(:) + i.*fcV10_lp(:)];
    [b_currD_lp] = [b_currD_lp(:); fcU1_lp(:) + i.*fcV1_lp(:)];
end

%Do the actual analysis and plotting
timeVars = {'time5'; 'time10'; 'time1'; 'aTime'; 'aTime'; 'aTime'; 'bTime'; ...
        'bTime'; 'bTime'; 'time5_lp'; 'time10_lp'; 'time1_lp'; 'aTime_lp'; ...
        'aTime_lp'; 'aTime_lp'; 'bTime_lp'; 'bTime_lp'; 'bTime_lp'};
currVars = {'curr5'; 'curr10'; 'curr1'; 'a_currS'; 'a_currM'; 'a_currD'; ...
        'b_currS'; 'b_currM'; 'b_currD'; 'curr5_lp'; 'curr10_lp'; 'curr1_lp'; ...
        'a_currS_lp'; 'a_currM_lp'; 'a_currD_lp'; 'b_currS_lp'; 'b_currM_lp'; ...
        'b_currD_lp'};
site = {'A'; 'A'; 'A'; 'A'; 'A'; 'A'; 'B'; 'B'; 'B'; 'A'; 'A'; 'A'; 'A'; 'A'; 'A'; 'B'; 'B'; 'B'};
depth = {'5 mbs'; '10 mab'; '1 mab'; '5 mbs'; '15 mbs'; 'Deepest Bin'; '5 mbs'; ...
        '10 mbs'; 'Deepest Bin'; '5 mbs'; '10 mab'; '1 mab'; '5 mbs'; '15 mbs'; ...
        'Deepest Bin'; '5 mbs'; '10 mbs'; 'Deepest Bin'};
timeStep = {'Hourly Averaged'; 'Hourly Averaged'; 'Hourly Averaged'; ...
        'Hourly Averaged'; 'Hourly Averaged'; 'Hourly Averaged'; ...
        'Hourly Averaged'; 'Hourly Averaged'; 'Hourly Averaged'; ...
        'Low-passed'; 'Low-passed'; 'Low-passed'; ...
        'Low-passed'; 'Low-passed'; 'Low-passed'; ...
        'Low-passed'; 'Low-passed'; 'Low-passed'};
printTitle = {'a_5_h_'; 'a_10_h_'; 'a_1_h_'; 'a_5a_h_'; 'a_15a_h_'; ...
        'a_1a_h_'; 'b_5a_h_'; 'b_10a_h_'; 'b_1a_h_'; 'a_5_lp_'; 'a_10_lp_'; ...
        'a_1_lp_'; 'a_5a_lp_'; 'a_15a_lp_'; 'a_1a_lp_'; 'b_5a_lp_'; ...
        'b_10a_lp_'; 'b_1a_lp_'}; 


for index = 16:18 %1:length(timeVars)
    for indexYear = 1:length(theYears)
        targetYear = theYears(indexYear);
        
        eval(['theTime = ' timeVars{index} ';']);
        eval(['theU = real(' currVars{index} ');']);
        eval(['theV = imag(' currVars{index} ');']);
        
        gregTime = gregorian(theTime);
        inYear = find(gregTime(:,1) == targetYear);
        if isempty(inYear)
            continue
        end
        figure, orient landscape
        
        [cdir_x, cspd] = cart2pol(theU(inYear), theV(inYear));
        cdir = angleXtoCompass(cdir_x*180/pi);
        lessZero = find(cdir <= 0);
        cdir(lessZero) = 360 + cdir(lessZero);
        cspd = varify(cspd(:));
        cdir = varify(cdir(:));
        
        xbin = [0:20:360];
        ybin = [0:2:50];
        
        %hspotx = zeros(length(ybin)-1,length(xbin)-1);
        %hspoty = zeros(length(ybin)-1,length(xbin)-1);
        
        %for j = 1:length(ybin),
        %    for i = 1:length(xbin),
        %        
        %        hspotx(j,i) = xbin(i);
        %        hspoty(j,i) = ybin(j);
        %        
        %    end
        % end
        
        hspotx = repmat([0:20:360], length(ybin), 1);
        hspoty = repmat([0:2:50]', 1, length(xbin));
        
        
        [a,h] = binsort(cdir(:),cspd(:),cspd(:),xbin(:),ybin(:));
        
        d_xbin = (mean(diff(xbin))/2);
        d_ybin = (mean(diff(ybin))/2);
        
        [m,n] = size(hspotx);
        
        a2 = a(2:end,2:end);
        h2 = h(2:end,2:end);
        
        hmax = max(max(h));
        hmin = min(min(h));
        
        cm = colormap;
        cm(1,:) = [1 1 1];
        colormap(cm)
        
        axis([0 380 0 52])
        
        h2(end+1,:) = 0;
        h2(:,end+1) = 0;
        
        
        theRange = [hmin hmax];
        theSpacing = diff(theRange)/64;   %Find the colorbar divisions
        
        toPlot = h2;
        toBump = find(toPlot > 0 & toPlot <= (hmin+theSpacing));
        toPlot(toBump) = hmin+theSpacing;
        
        H = surface(hspotx,hspoty,toPlot)
        set(H, 'edgealpha', 0.9)
        shading faceted
        ax1 = gca;
        set(ax1,'ydir','reverse')
        set(ax1,'Box','on')
        set(ax1,'XTick',xbin,'FontSize',[8])
        set(ax1,'YTick',ybin,'FontSize',[8])
        
        
        %hpos = find(h2 > (hmin+theSpacing));
        hpos = find(h2 > 0);
        hfin = h2(hpos);
        hxpos = hspotx(hpos) + (d_xbin);
        hypos = hspoty(hpos) + (d_ybin);
        
        
        
        for i = 1:length(hpos),
            hname = hfin(i);
            if hname < 1
                hval = sprintf('%4.2f',hname);
            else
                hval = sprintf('%4.1f',hname);
            end
            
            if hname > (hmin+10*theSpacing) & hname < (hmin+60*theSpacing)
                t(i) = text(hxpos(i),hypos(i),hfin(i)+1,hval,'color',[0 0 0],'HorizontalAlignment','center',...
                    'VerticalAlignment','middle','FontSize',[7],'FontWeight','bold');
            else
                t(i) = text(hxpos(i),hypos(i),hfin(i)+1,hval,'color',[1 1 1],'HorizontalAlignment','center',...
                    'VerticalAlignment','middle','FontSize',[7],'FontWeight','bold');
            end
            set(t(i), 'erase', 'none')
        end
        
        for i = 1:size(h2,2)-1
            text(i*20-10, 51, sprintf('%4.2f', sum(h2(:,i))), 'color', 'k', ...
                'HorizontalAlignment','center', 'VerticalAlignment','middle', ...
                'FontSize',[7],'FontWeight','bold')
        end
        
        for i = 1:size(h2,1)-1
            text(370, i*2-1, sprintf('%4.2f', sum(h2(i,:))), 'color', 'k', ...
                'HorizontalAlignment','center', 'VerticalAlignment','middle', ...
                'FontSize',[7],'FontWeight','bold')
        end
        
        text(370, 51, sprintf('%4.2f', sum(sum(h2))), 'color', 'k', ...
            'HorizontalAlignment','center', 'VerticalAlignment','middle', ...
                'FontSize',[7],'FontWeight','bold')
                
        
        title({['Site ' upper(site{index}) ', Sorted Current Speed and Direction']; ...
                [timeStep{index} ' Current, ' depth{index} ': ' num2str(targetYear, '%.0f')]})
        
        cb = colorbar;
        
        if rem(targetYear,4) == 0
            yearDays = 366;
        else
            yearDays = 365;
        end
        
        if strcmp(timeStep{index}, 'Hourly Averaged')
            totalNumber = 24*yearDays;
        else
            totalNumber = 24*yearDays/6;
        end
        
        percentGood = 100*length(find(~isnan(cspd)))/totalNumber;
        
        xlabel({'Current Direction (compass degrees)'; ...
                ['Total Points Sorted = ' num2str(length(find(~isnan(cspd))), '%.0f') ...
                    ', ' num2str(percentGood, '%.1f') '% of Possible Observations']})
        ylabel('Current Speed (cm/s)')
        axes(cb)
        
        
        ylabel('Percent of data')
        print('-dpdf', fullfile(plotDirectory, [printTitle{index} num2str(targetYear, '%.0f')]));
        close
    end
end
