
function [theResult, duplicate_S] = average_and_remove_duplicates_aoml4 (output, outfile)
%
% Find duplicate samples in output array
%
% output = seven or eight column array of... 
%   [output_sta output_niskin output_bottle output_twiceratio output_salinity 
%          1         2            3             4                 5 
%                      output_quality output_cast (OPTIONAL) salinity_correction];
%                           6            7                         8
%
% this routine looks for duplicate output_niskin for each output_sta
% Averages these duplicates together provided they are less than 2*stddev
% of all duplicates away from the mean duplicate difference (i.e. tries to remove outliers
% before averaging).
%
% Output:
% outfile = text file of seven or eight columns as in output array, but with
% duplicates removed.
% Averaged duplicates qC flags are replaced with a "6"
%
% Version 1 called average_and_remove_duplicates for clivar with columns of
%   output different
% Version 2 called average_and_remove_duplicates_aoml for the format as
%   described above.
% Version 3 checks to see if output is 7 or 8 columns wide and deals with
%   the correct print statement for each case.
% Version4 only checks for nine columns that also contain cast number
%   If cast number is present, then only averages duplicates if cast and
%   station number match.
%
Icolumn_S = 5;
Icolumn_ratio = 4;
Icolumn_quality = 6;
Icolumn_salinity_correction = 8;
Icolumn_castnumber = 9;

[pathstr,this_name,ext] = fileparts(outfile);

duplicate_S = [];
index_dups = find_repeats2 (output);
% note that any index may appear more than once in this index_dups if more
% than one duplicate is drawn from a Niskin bottle


if isempty(index_dups)
    fprintf(1,['\nNo Duplicates Found.  Writing ' outfile '\n\n'])
else
    duplicate_S     = zeros (size(index_dups,1), size(index_dups,2));
    duplicate_ratio = zeros (size(index_dups,1), size(index_dups,2));
    for II = 1:size(index_dups,1)
        duplicate_S (II,:) = output(index_dups(II,:),Icolumn_S);
        duplicate_ratio (II,:) = output(index_dups(II,:),Icolumn_ratio);
    end
    if (size(output,2) == 8)
        duplicate_salinity_correction = zeros (size(index_dups,1), size(index_dups,2));
        for II = 1:size(index_dups,1)
            duplicate_salinity_correction(II,:) = output(index_dups(II,:),Icolumn_salinity_correction);
        end
    end
     
    
    %
    % Average duplicates and plot differences
    %
    %    plot(mean(duplicate_S'),diff(duplicate_S'),'+')
    diff_dup_S = diff(duplicate_S');
    Index_good_dups = find(abs(diff_dup_S) <= median(diff_dup_S)                  + 2*std(diff_dup_S));
    Index_good_dups = find(abs(diff_dup_S) <= median(diff_dup_S(Index_good_dups)) + 2*std(diff_dup_S(Index_good_dups)));
    Index_good_dups = find(abs(diff_dup_S) <= median(diff_dup_S(Index_good_dups)) + 2*std(diff_dup_S(Index_good_dups)));

    figure;
        number_dups = [1:length(diff_dup_S)];
        plot(number_dups, diff_dup_S,'+')
        hold on;
        plot(number_dups(Index_good_dups), diff_dup_S(Index_good_dups), 'r+')
        [xlimits]=get(gca,'xlim');
        plot(xlimits, [median(diff_dup_S(Index_good_dups)) median(diff_dup_S(Index_good_dups))],'k--')
        plot(xlimits, [median(diff_dup_S(Index_good_dups))+2*std(diff_dup_S(Index_good_dups)) median(diff_dup_S(Index_good_dups))+2*std(diff_dup_S(Index_good_dups))],'r:')
        plot(xlimits, [median(diff_dup_S(Index_good_dups))-2*std(diff_dup_S(Index_good_dups)) median(diff_dup_S(Index_good_dups))-2*std(diff_dup_S(Index_good_dups))],'r:')
        [ylimits] = get(gca,'ylim');
        text([xlimits(1)+0.1*diff(xlimits)], [ylimits(2)-0.1*diff(ylimits)], ['Median = ' num2str(median(diff_dup_S(Index_good_dups)),5) ' +/- ' num2str(std(diff_dup_S(Index_good_dups)),5)], 'fontsize', 14)
        title (['Duplicate Salinity from ' upper(strrep(this_name,'_',' '))], 'fontsize', 14)
        set(gca,'xlim',[0.5 length(diff_dup_S)+0.5])
    %
    % Average good duplicates and put the average value in the first
    % duplicate location and remove the second one
    %
    avg_S = meanmiss(duplicate_S')';
    avg_ratio = meanmiss(duplicate_ratio')';
    for II = 1: length(Index_good_dups);
        index_dups(Index_good_dups(II),1);
        JJ = find(index_dups(:,1) == index_dups(Index_good_dups(II),1));
        if (length(JJ) == 1)
            output (index_dups(Index_good_dups(II),1),Icolumn_ratio)   = avg_ratio(Index_good_dups(II));
            output (index_dups(Index_good_dups(II),1),Icolumn_S)       = avg_S(Index_good_dups(II));
            output (index_dups(Index_good_dups(II),1),Icolumn_quality) = 6;
            if (size(output,2) == 8)
                output (index_dups(Index_good_dups(II),1),Icolumn_salinity_correction) = meanmiss(duplicate_salinity_correction(Index_good_dups(II),:)')';
            end
        else
            % Then there is more than one row with duplicate values
            temp_index = unique(reshape(index_dups(JJ,:),size(JJ,1)*2,1))
            output (temp_index(1),Icolumn_ratio)   = meanmiss(output(temp_index,Icolumn_ratio));
            output (temp_index(1),Icolumn_S)       = meanmiss(output(temp_index,Icolumn_S));
            output (temp_index(1),Icolumn_quality) = 6;
        end
            
    end;
      
    %
    % Remove the duplicates (by convention remove the second duplicate).
    %
    output(index_dups(:,2),:) = [];
end

fid=fopen(outfile,'wt');
if (size(output,2) == 7)
    fprintf(fid,'%Station\t Cast\t Niskin\t Sample_Bottle\t 2*Cond. Ratio Corrected Salinity\n');
    fprintf(fid,'%d\t%d\t%d\t%d\t%1.5f\t%2f5\t%d\n',output');
else
    fprintf(fid,'%Station\t Cast\t Niskin\t Sample_Bottle\t 2*Cond. Ratio Corrected Salinity\n Salinity Correction \n');
    fprintf(fid,'%d\t%d\t%d\t%d\t%1.5f\t%2f5\t%d\t%d\n',output');
end
      

fclose(fid);

theResult = output;
return;


