function theResult  = correct_autosal_drift2_matlab_structure_no_file (nominal_std, salts, outfile)
% CORRECT_AUTOSAL_DRIFT2_MATLAB_STRUCTURE - correct for autosal drift with pre/post standards 
% 
% CTD Calibration toolbox
% 
% INPUT: 
%   nominal_cond_standard_water: scalar representing the nominal value for
%		standard water conductivity ratio.  This number changes
% 		for each batch of standard seawater. 
%       e.g. nominal_std=2*0.99981;
%       Note:  it is assumed that the "nominal value is the same units as provided in
% 		the input file.
%   salts: structure of input data from the autosal interface
%   outfile: path and name of the output file which will contain the 
%		corrected averaged bottle sample salinity.
%       e.g. outfile = 'a10_014_015_out2.txt';
%
% OUTPUT:
%   theResult: the data structure of corrected values.  This is
%		essentially a duplicate of the information in the infile with an
%		additional varialbe that is the corrected value.
%
% DESCRIPTION:	
% 	This function is specific to the autosal configuration found on the
%	Ron Brown during the 2011 field season.  Modifications may be
%	necessary, but please rename accordingly.
%   This is just like correct_autosal_drift2.m except does not read in txt
%   files seperately.  Assumes txt files have already been read in and
%   placed in structure.
%
% SEE ALSO:
%		read_autosal_dat_raw		
%
% 
% CHANGELOG: 
% Sept/2012 - wrote this
% Feb 2012 - check to make sure plot subdirectory exists.  If not, create
%               it.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


% nominal_std = 2*0.99981
% outfile = 'test.lis'

warning off MATLAB:divideByZero
[pathstr,this_name,ext,versn] = fileparts(outfile);

%
% check that plot subdirectory exists.  If not, create it.
if ~exist('plot','dir')
    mkdir plot
end


fprintf (1, ['\n\n\n Processing ' salts.file '\n\n\n']);

%
% Figure out indeces of the standards
% Note the first standard must be labeled "1000" and the closing standard for
% each run must be labelled "1001".
%
standard1_idx=find(salts.txt_station_id==1000);
standard2_idx=find(salts.txt_station_id==1001);
standard1_idx_all = standard1_idx;
standard2_idx_all = standard2_idx;

%
% Give all stanards a "bad" flag" - reset flag after QC
%
salts.txt_qc(standard1_idx) = 4;
salts.txt_qc(standard2_idx) = 4;

%
% QC standard values
% 
% First get only the best 3 values
%
while length(standard1_idx) > 3
    [value index_max]=max(abs(salts.txt_cond_ratio(standard1_idx') - median(salts.txt_cond_ratio(standard1_idx'))));
    standard1_idx(index_max(1)) = [];
end
%
% Next check that there is not still one value that is an outlier
%
if max(salts.txt_cond_ratio(standard1_idx))-min(salts.txt_cond_ratio(standard1_idx)) > 0.00002
    fprintf (1,['WARNING (' salts.file ') :  Range of Initial Standard water values is too large \n   = ' num2str(std(salts.txt_cond_ratio(standard1_idx')),13) ...
            '\n   Removing outlier and recompute.\n \n'])
    % Then check the last two values
    if (length(standard1_idx) >= 3)
        if abs(salts.txt_cond_ratio(standard1_idx(3)) - salts.txt_cond_ratio(standard1_idx(2))) < 0.00002
            % then use only the last two values
            standard1_idx(1)=[];
        else
            % get rid of the outlier whereever it is
            [value index_max] = max(abs(salts.txt_cond_ratio(standard1_idx') - median(salts.txt_cond_ratio(standard1_idx'))));
            standard1_idx(index_max(1)) = [];
        end
    end
end

%
% Now do the same thing for the final standards
%
while length(standard2_idx) > 3
    [value index_max] = max(abs(salts.txt_cond_ratio(standard2_idx') - median(salts.txt_cond_ratio(standard2_idx'))));
    standard2_idx(index_max(1)) = [];
end
%
% Next check that there is not still one value that is an outlier
%
if max(salts.txt_cond_ratio(standard2_idx))-min(salts.txt_cond_ratio(standard2_idx)) > 0.00002
    fprintf (1,['WARNING (' salts.file ') :  Range of Final Standard water values is too large \n   = ' num2str(std(salts.txt_cond_ratio(standard1_idx')),13) ...
            '\n   Removing outlier and recompute.\n \n'])
    % Then check the last two values
    if (length(standard2_idx) >= 3)
        if abs(salts.txt_cond_ratio(standard2_idx(3))- salts.txt_cond_ratio(standard2_idx(2))) < 0.00002
            % then use only the last two values
            standard2_idx(1)=[];
        else
            % get rid of the outlier whereever it is
            [value index_max]=max(abs(salts.txt_cond_ratio(standard2_idx') - median(salts.txt_cond_ratio(standard2_idx'))))
            standard2_idx(index_max(1)) = [];
        end
    end
end

% 
% Set standards that were actually used to "2" = good
%
salts.txt_qc(standard1_idx) = 2;
salts.txt_qc(standard2_idx) = 2;

%
% Make a variable containing only the standards
% x = the date/time of the standards
% y = the uncorrected standard ratios
%
x = salts.txt_date ([standard1_idx' standard2_idx']');
y = salts.txt_cond_ratio ([standard1_idx' standard2_idx']');
x_all = salts.txt_date ([standard1_idx_all' standard2_idx_all']');
y_all = salts.txt_cond_ratio ([standard1_idx_all' standard2_idx_all']');

%
% Make sure they are column vectors
%
x=x(:);
y = y(:);
x_all = x_all(:);
y_all = y_all(:);

%
% make sure there are no NaNs
%
II = find(isnan(y));
x(II)=[]; y(II) = [];

%
% Determine the standard water correction that needs to be applied.
%
if isempty(standard1_idx)
    fpprintf (1, ['WARNING (' salts.file ') :  No starting standard found\n'])
    % Look for ending standard
    if isempty(standard2_idx)
        fpprintf (1, ['WARNING (' salts.file ') : No standards at all are found.  Cannot adjust values.\n'])
        btl_correction = zeros(length(salts.txt_cond_ratio),1);
    else
        fpprintf (1, ['WARNING (' salts.file ') : Only ending Standard is found.  Can only adjust mean offset for all bottles.\n'])
        btl_correction  = (mean  (salts.txt_cond_ratio(standard2_idx')) - nominal_std) * ones(length(salts.txt_cond_ratio),1);
        btl_correction2 = (median(salts.txt_cond_ratio(standard2_idx')) - nominal_std) * ones(length(salts.txt_cond_ratio),1);
        %
        % Do some very simple error checking:
        %
        if (std(salts.txt_cond_ratio(standard2_idx')) > 0.00001*2)
             fprintf (1,['WARNING (' salts.file ') :  Ending Standard water values variance too large \n   = ' num2str(std(salts.txt_cond_ratio(standard2_idx')),13) '\n'])
        end
    end
else
    if isempty(standard2_idx)
       fpprintf (1, ['WARNING (' salts.file ') : Only beginning Standard is found.  Can only adjust mean offset for all bottles.\n'])
       btl_correction  = (mean  (salts.txt_cond_ratio(standard1_idx')) - nominal_std) * ones(length(salts.txt_cond_ratio),1);
       btl_correction2 = (median(salts.txt_cond_ratio(standard1_idx')) - nominal_std) * ones(length(salts.txt_cond_ratio),1);
        %
        % Do some very simple error checking:
        %
        if (std(salts.txt_cond_ratio(standard2_idx')) > 0.00001*2)
             fprintf (1,['WARNING (' salts.file ') :  Ending Standard water values variance too large \n   = ' num2str(std(salts.txt_cond_ratio(standard2_idx')),13) '\n'])
        end
    else
        % Can do a linear offset correction as a funciton of time
        %
        % fit a linear poynomial
        %
        % Many ways to do this e.g.:
        %   [coeff, S, tmp_best_fit_Rsq, tmp_best_fit_f, t] = molly_polyfit (x, y, 1, 1, infile);
        %   calculate_and_plot_linear_fit2 (x, y, ['standards for ' infile], 'date', ' autosal ratio', [infile '.ps'], 1)
        %
        %   Most robust method computes all possible varieties: do_mbari_fits2 (below)
        %       Note that the third and final argument is the unit number to write
        %       all the results to.  Unit=1 is the screen.  Use this unless you
        %       have already opened a logfile (otherwise there will be an error message).
        %
        fprintf(1,'\nUsing ALL beginning and ending standards with least squares fit:\n');
        [m, b, r, sm, sb, rms_err, mean_per_err, mean_per_err_2] = do_mbari_fits2 (x, y, 1);
        % 
        % Also could make one median start and stop standard:
        fprintf(1,'\nUsing Median of all beginning and ending standards only:\n');
        x2 = [ mean(salts.txt_date(standard1_idx')) mean(salts.txt_date(standard2_idx'))];
        y2 = [ median(salts.txt_cond_ratio(standard1_idx')) median(salts.txt_cond_ratio(standard2_idx'))];
        [m2, b2, r2, sm2, sb2, rms_err2, mean_per_err2, mean_per_err_22] = do_mbari_fits2 (x2, y2, 1);
        
        % Also could make one mean start and stop standard:
        fprintf(1,'\nUsing Mean of all beginning and ending standards only:\n');
        x3 = [ mean(salts.txt_date(standard1_idx')) mean(salts.txt_date(standard2_idx'))];
        y3 = [ mean(salts.txt_cond_ratio(standard1_idx')) mean(salts.txt_cond_ratio(standard2_idx'))];
        [m3, b3, r3, sm3, sb3, rms_err3, mean_per_err3, mean_per_err_23] = do_mbari_fits2 (x3, y3, 1);
        
   
        m =fillmiss (m, NaN, 10, 0.);
        m2 =fillmiss (m2, NaN, 10, 0.);
        m3 =fillmiss (m3, NaN, 10, 0.);
        b  = fillmiss  (b, NaN, 10, mean(y));
        b2 = fillmiss  (b2, NaN, 10, mean(y2));
        b3 = fillmiss  (b3, NaN, 10, mean(y3));       
         
        %
        %  for the standards, the correction should give you
        %  y_orig - yfit = nominal_stnd
        %  where yfit = m*x + b
        %  hence y_new = y_orig - correction
        %  where correction = (yfit - nominal_stnd)
        %  y_new = y_orig - (yfit-nominal_stnd) = y_orig - yfit + nominal_stnd
        %
        % btl_correction = m(5)* date_cond + b(5) - nominal_std;
        % After looking at the first set of standards through station 36,
        % it looks best to use a median fit for all.
        
        btl_correction = m(5)* salts.txt_date + b(5) - nominal_std;
        
        %
        % Do some very simple error checking:
        %
        if (std(salts.txt_cond_ratio(standard1_idx')) > 0.00001*2)
             fprintf (1,['WARNING (' salts.file ') :  Beginning Standard water values variance too large \n   = ' num2str(std(salts.txt_cond_ratio(standard1_idx')),13) ...
                 '\n   Consider removing outlier and recompute.\n For now use slope from median standards.\n'])
             btl_correction = m2(5)* salts.txt_date + b2(5) - nominal_std;
        end
        if (std(salts.txt_cond_ratio(standard2_idx')) > 0.00001*2)
             fprintf (1,['WARNING (' salts.file ') :  Ending Standard water values variance too large \n   = ' num2str(std(salts.txt_cond_ratio(standard1_idx')),13) ...
                 '\n   Consider removing outlier and recompute.\n For now use slope from median standards.\n'])
             btl_correction = m2(5)* salts.txt_date + b2(5) - nominal_std;
        end
        if ((max(btl_correction) - min(btl_correction)) > 0.00006)
            fprintf (1,['WARNING (' salts.file ') :  Standard water values changed too much during run.  \n   '  ... 
                'Check autosal peformance, room temperature, calibration values and scatter before and after run.\n'])
        end
        if (abs(m2(5)) > (abs(m(5))+sm(5)))
            fprintf (1,['WARNING (' salts.file ') :  Differences between slopes using Median and Mean reads is out of bounds.  \n   ' ...
                'Check calibration values scatter.  \n   Consider removing outlier and recompute.\n For now use slope from median standards.\n'])
            btl_correction = m2(5)* salts.txt_date + b2(5) - nominal_std;
        end

        
        figure;
            subplot(2,1,1);
                plot(x_all,100*y_all,'ok','Markerfacecolor','k','MarkerSize',10)
                hold on;
                plot(x, 100*y,'ok','Markerfacecolor','b','MarkerSize',10)
                plot(x, 100*( m(5)*x+b(5)), 'b', 'linewidth', 2)
                plot(x, 100*(m2(5)*x+b2(5)), 'k--', 'linewidth', 1.5)
                plot(x, 100*(m3(5)*x+b3(5)), 'r:', 'linewidth', 1.5)
                legend ('All Standards', 'QC''d Standards', 'Fit from LSQ Reads', 'Fit from Median Reads', 'Fit from Mean Reads','Location','Best') 
                datetick
                title (['Twice Standard Ratio times 100: ' upper(strrep(salts.file,'_',' '))])
  %              set(gca,'Box','off');
  %              axesposition=get(gca,'Position');
  %              ylimits = get(gca,'ylim');
  %              hNewAxes = axes('Position', axesposition, 'color','none','Ylim',[sw_sals(ylimits(1)/100/2,24) sw_sals(ylimits(2)/100/2,24)],'Yaxislocation','right','xtick',[],'box','off');

            subplot(2,1,2);
                plot(x_all, 100000*(y_all-(m(5)*x_all+b(5))), 'ok','Markerfacecolor','k','MarkerSize',10)
                hold on;
                plot(x,     100000*(y              -(m(5)*x+b(5))), 'ok','Markerfacecolor','b','MarkerSize',10)
                plot(x,     100000*((m(5)*x+b(5))  -(m(5)*x+b(5))), 'b', 'linewidth', 2)
                plot(x,     100000*((m2(5)*x+b2(5))-(m(5)*x+b(5))), 'k--', 'linewidth', 1.5)
                plot(x,     100000*((m3(5)*x+b3(5))-(m(5)*x+b(5))), 'r--', 'linewidth', 1.5)
                datetick
                title ('Residual times 100000')
  %              set(gca,'Box','off');
  %              axesposition=get(gca,'Position');
  %              ylimits = get(gca,'ylim');
  %              hNewAxes = axes('Position', axesposition, 'color','none','Ylim',[sw_sals((nominal_std+ylimits(1)/100000)/2,24) sw_sals((nominal_std+ylimits(2)/100000)/2,24)],'Yaxislocation','right','xtick',[],'box','off');
                eval(['print -depsc plot/cal_' salts.file '.ps'])
                pause
 %      figure;
 %           plot(salts.txt_date, 100*salts.txt_cond_ratio, 'ok','Markerfacecolor','b','MarkerSize',10)
 %           hold on;
 %           plot(salts.txt_date, 100*(salts.txt_cond_ratio-btl_correction),' ok','Markerfacecolor','r','MarkerSize',5)
 %           plot(x, 100*(y - (m(5)*x + b(5) - nominal_std)), ' *k','Markerfacecolor','g','MarkerSize',5) 
 %           datetick
 %           title (['Twice Standard Ratio times 100' upper(strrep(salts.file,'_',' '))])
 %           set(gca,'Box','off');
 %           axesposition=get(gca,'Position');
 %           ylimits = get(gca,'ylim');
 %           hNewAxes = axes('Position', axesposition, 'color','none','Ylim',[sw_sals((ylimits(1)/100)/2,24) sw_sals((ylimits(2)/100)/2,24)],'Yaxislocation','right','xtick',[],'box','off');
    end
end

btl_correction = btl_correction(:);
cond_ratio2_correct = salts.txt_cond_ratio-btl_correction;

output_correct=[salts.txt_station_id, salts.txt_niskin_bottle , ...
                salts.txt_samp_nbr  , cond_ratio2_correct     ,...
                sw_sals(cond_ratio2_correct/2, salts.txt_tank_temp)];

%fid=fopen(outfile,'wt');
%fprintf(fid,'%Station\tNiskin\tSample_Bottle\t2*Cond. Ratio\tCorrected Salinity\n');
%fprintf(fid,'%d\t1\t%d\t%d\t%1.5f\t%2f5\t2\n', output_correct');
%fclose(fid);

theResult = salts;
theResult.txt_cond_ratio_correct = cond_ratio2_correct;
theResult.txt_salinity_correct = sw_sals(cond_ratio2_correct/2, salts.txt_tank_temp);

return