function [sta,slope,bias,newbotco,newctdco]=calcos1(ctdst,ctdco,ctdpr,botco,stddev,stamin,stamax);

% function [sta,slope,bias,newbotco,newctdco]=calcos1(ctdst,ctdco,ctdpr,botco,stddev,stamin,stamax);
% fed vectors of station number [ctdst], ctd conductivity [ctdco],
% ctd pressure [ctdpr],
% bottle conductivity [botco], the number of standard deviations 
% beyond which to reject outliers [stddev], as well as a minimum (stamin) 
% and maximum [stamax] station number, the routine spits back the a vector 
% of stations [sta], a station-dependent conductivity slope [slope], a 
% single conductivity bias value [bias], the bottles with outliers 
% rejected [newbotco], and the calibrated ctd values [newctdco].  
% uses a linear (1st order) station-dependent fit for co slope.

% set up various indices

ctdst=ctdst(:); ctdco=ctdco(:); botco=botco(:);
ctdones=1+0*ctdco; rejects=[]; count=0; 
tot=length(find(finite(botco)==1));

jj=1;

% repeat fit until standard deviation criteria fits
% store rejects indices in rejects


while length(jj)>0;

% fit slope adjustment as a linear function of station number
% with a uniform BIAS

count=count+1;
ii=find(finite(botco)==1);
c=[ctdco(ii) ctdco(ii).*ctdst(ii) ctdones(ii)]\botco(ii);
c=c';

% calculate standard deviation and remove outliers from 
% the bottle file

d=c(1)*ctdco+c(2)*ctdco.*ctdst+c(3)*ctdones-botco;
s=std(d(ii));
jj=find(abs(d)>stddev*s);
botco(jj)=NaN*jj;

end % while

% find reject indices

% get station numbers, but only one for each station

sta=[round(stamin):round(stamax)];

% calculate conductivity slope and bias for each station

slope=[c(1)+c(2)*sta];
bias=c(3);

% display some statistics

sub=length(find(finite(botco)==1));
disp(['number of points used   ',int2str(sub)]);
disp(['total number of points  ',int2str(tot)]);
disp(['% of points used in fit ',num2str(100*sub/tot,4)]);
disp(['fit standard deviation  ',num2str(s,4)]);
disp(['fit bias                ',num2str(bias,8)]);
disp(['min fit slope           ',num2str(min(slope),8)]);
disp(['max fit slope           ',num2str(max(slope),8)]);

% get new bottle and ctd conds

newbotco=botco;
newctdco=c(1)*ctdco+c(2)*ctdco.*ctdst+c(3);

% plot the fit

%close all;                         % for single plots
%figure; startup;                   % for single plots 
figure; hold on;                    % use when over plotting 
orient portrait;
plot(ctdst,newctdco-botco,'ko');
%a = axis;                          % for single plots 
%axis([a(1) a(2) -.03 .03]);        % for single plots 
grid on;
xlabel('Station Number');
ylabel('Calibrated CTD-Bottle conductivity');

%figure; startup;                   % for single plots 
figure; hold on;                    % use when over plotting
orient portrait;
plot(ctdpr,newctdco-botco,'ko');
%b = axis;                          % for single plots 
%axis([b(1) b(2) -.03 .03]);        % for single plots 
grid on;
xlabel('Pressure [dbar]');
ylabel('Calibrated CTD-Bottle conductivity');
