function [oxy_coeff,fh]=oxy_cal(cruiseid,st1,st2)
% OXY_CAL - oxygen calibration calling routine.
%
% CTD Calibration toolbox.
%
% USAGE
% 	[oxy_coeff,fh]=oxy_cal(cruiseid,st1,st2);
%
% INPUT
%	cruiseid	cruise identification name
% 	st1		first station
%	st2		last station
%
% OUTPUT
%	oxy_coeff	oxygen calibration coefficients
%	fh		figure handles (in case we need them again)
%
% DESCRIPTION
%	still to write
%

% AUTHORS:
%			Carlos Fonseca
%         UM/CIMAS
%         Derrick Snowden
%         NOAA/AOML/PhOD
%         Tue May 11 11:23:37 EDT 2004
%
% VERSION:
%	$Revision: 1.23 $
%	$Date: 2008/04/14 19:45:11 $
%	$Id: oxy_cal.m,v 1.23 2008/04/14 19:45:11 duncombe Exp duncombe $
%
% CHANGELOG:
%   23-Jul-2004, Version 1.0
%	* Initial version.
%  2007/12/17 - CMDR
%	* read cruise id from the file in the current directory
% 	* test for secondary sensors. If they do not exist make a copy of the
%	primary sensors and use that.
%	* some statements look like they could be made more readable by
%	utilizing Matlab absolute value function (abs)
%	* use a function for oxygen saturation calculation
% 	* make input lines easier to read (aesthetic adjustment!)
%	* make functions oxy_cal and FS_apply_calibration have the same input
%	and output arguments to facilitate cut and paste
%	* return figure handles
%  2008/01/11 CMDR
%	* use the oxygen saturation function from the seawater toolbox
%  2008/01/19 CMDR
%	* put first and last stations in the argument list
%	* use nanmean for plotting residual means
%	* NIST recommended symbol for the non-SI unit of volume liter is
%	capital L (http://physics.nist.gov/cuu/Units/outside.html)
%  2008/01/20 CMDR
%	* flag missing secondary sensors and do not plot them.
%	* plot normalized differences against bottle number.
%  2008/02/18 CMDR
%	* append new oxygen calibration coefficients to the database
%  2008/03/10 CMDR
%	* provide some defaults for st1 and st2; vectorize some loops
%  2008/03/11 CMDR
%	* there is some weirdness going on with the figures. There seems
%	to be a figure statement left out (boxy-coxy was plotted on top
%	of coxy-boxy; and there was no hold on statement). Put in a
%	figure statement. Mean to look into the rest later.
%  2008/04/14 CMDR
%  	* test if quality fields are present. if not remind operator to run
%	compare_oxy.m
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%Parameters**************************
% Applied in lines 514 and 519
%Primary oxygen sensor 0730 ab1202
if 0, 
Soc = 0.4946;		
Voffset = -0.5102;	
tau20 = 1.16;
A = -3.5069e-3;
B = 1.3800e-4;
C = -2.6719e-6;
E = 0.036;
end;

%Secondary oxygen sensor 1266 ab1202
if 0,
Soc = 0.5110
Voffset = -0.5418
tau20 = 1.13
A = -3.0146e-3
B = 1.3378e-4
C = -2.6688e-6
E = 0.036
end;

% sensor 1329 marh 2009
if 1,
Soc=0.3694;		
Voffset=-0.4965;	
tau20=1.14;
A=-5.4290e-4;
B=1.4123e-4;
C=-2.5131e-6;
E=0.036;
end;
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%% CMDR read cruiseid if we can
if nargin==0, cruiseid=getcruiseid; end;
if isempty(cruiseid), cruiseid=getcruiseid; end;
%%%%%%%%%%%
fcnt=0; 	% count for the figures
t2flag=1;	% flags for 2nd sensors
s2flag=1;	%
o2vflag=1;	%
o2cflag=1;	%
c2flag=1;	%

% Retrieve the necessary preferences for this cruise
group = group_name(cruiseid);

% cruise_dir is the main calibration data directory

if ispref(group,'cal_data_dir')
	cruise_dir=getpref(group,'cal_data_dir');
else
	error(['The cal_data_dir preference was not set for cruise: ' ...
		cruiseid '.  Run register_cruise.m']);
end;

% Load data from the cruise database
db_file = fullfile(cruise_dir,[cruiseid '_db.mat']);
if exist(db_file)==2
	load(db_file);
else
	error(['The cruise data base file was not found on the path.  ' ...
		db_file]);
end;
n=who('oxy_match');

if isempty(n)==1
	error(['oxy_match not found. Run make_oxy_sumfile.m ']);
	return;
end;
warning off MATLAB:divideByZero;
% set up output variable
oxy_coeff=zeros(7,1);
stations=unique(oxy_match.stnnbr);
station_min=min(stations);
station_max=max(stations);
disp('**************************************************************')
disp('The data found in oxy_match goes from')
disp(['the station ',num2str(station_min),' to the station ', ...
	num2str(station_max),'.'] )
disp(['The maximum pressure is: ',num2str(max(oxy_match.ctdprs))])
disp('**************************************************************')
%
% Station interval
%%%%%%%%%%%% CMDR
% st1=input('Enter first station number: ');
% st2=input('Enter last station number: ');
if ~exist('st1'), st1=input('Enter first station number: '); end;
if ~exist('st2'), st2=input('Enter last station number: '); end;
if isempty(st1), st1=station_min; end;
if isempty(st2), st2=station_max; end;
%%%%%%%%%%%%
% preliminary qc of the data using the temperature
disp('*********************************')
disp('Preliminary inspection of the data based on temperature')
disp('*********************************')

% %% CMDR test if the quality fields are present
% if ~any(isfield(oxy_match,{'quality1','quality2'})),
% 	error('oxy_match.quality fields are missing. Run compare_oxy.m');
% 	return;
% end;
%%

% loading the data
    ii=find(oxy_match.stnnbr>=st1 & oxy_match.stnnbr<=st2  & ...
		oxy_match.quality1==2 & oxy_match.quality2==2);
  % common variables
    stn=oxy_match.stnnbr(ii);
    nis=oxy_match.niskinbtlpos(ii);
    o2_bot_num=oxy_match.sampnbr(ii);
    cprs=oxy_match.ctdprs(ii);
    boxy=oxy_match.btloxy(ii);
    tot=length(boxy);
  % primary sensors
     ctem1=oxy_match.ctdtmp1(ii);
     csal1=oxy_match.ctdsal1(ii);
     coxv1=oxy_match.ctdoxv1(ii);
     coxy1=oxy_match.ctdoxy1(ii);
     ccon1=oxy_match.ctdcon1(ii);
  % secondary sensors
%%%%%%%%%%% CMDR test for second sensors. Set flags about the 2nd sensors
% ctem2
     if ~isnan(oxy_match.ctdtmp2), ctem2=oxy_match.ctdtmp2(ii);
		else ctem2=NaN; t2flag=0; end;
% csal2
     if ~isnan('oxy_match.ctdsal2'), csal2=oxy_match.ctdsal2(ii);
		else csal2=NaN; s2flag=0; end;
% coxv2
     if  ~isnan('oxy_match.ctdoxv2'), coxv2=oxy_match.ctdoxv2(ii);
		else  coxv2=NaN; o2vflag=0; end;
% coxy2
     if ~isnan(oxy_match.ctdoxy2), coxy2=oxy_match.ctdoxy2(ii);
		else  coxy2=NaN; o2cflag=0; end;
% ccon2
     if ~isnan(oxy_match.ctdcon2), ccon2=oxy_match.ctdcon2(ii);
		else  ccon2=NaN; c2flag=0; end;
%%%%%%%%%
    coxy_ori=coxy1;
  % First look in the data
     % Primary Sensor
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
    plot(stn,coxy1-boxy,'r+');  % plot diff betw 1st sensor and bottles
    grid on;
    hold on;
    for i=st1:1:st2,
    	aux=find(stn==i);
%%%%%% CMDR use nanmean to prevent any NaNs spoiling the result
    	plot(i,nanmean(coxy1(aux)-boxy(aux)),'b^','MarkerFaceColor','b');
    end;
    xlabel('Station number');
%% CMDR NIST recommended symbol for liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title('Primary Sensor - Uncalibrated Data');
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
    plot(cprs,coxy1-boxy,'ro');  % plotting diff betw 1st sensor and bottles
%% CMDR NaNs make a mess of the result
    diff_mean=nanmean(coxy1-boxy);
    diff_stddev=nanstd(coxy1-boxy);
%%
    grid on;
    xlabel('Pressure [dbar]');
%% CMDR NIST recommended symbol for liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title(['Primary Sensor - Mean Diff: ',num2str(diff_mean), ...
		' Std Dev: ',num2str(diff_stddev)]);
    disp('Raw data - Primary Sensor');
    disp(['Mean residual: ',num2str(diff_mean)]);
    disp(['Standard Deviation: ',num2str(diff_stddev)]);
    disp('Hit any key to continue');
    pause;
    % Secondary Sensor
%%%%%%%%%%% CMDR test for second sensors
if o2cflag,
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
    plot(stn,coxy2-boxy,'r+');
    grid on;
    hold on;
    for i=st1:1:st2,
    aux=find(stn==i);
%% CMDR ignore nans they make a mess
    plot(i,nanmean(coxy2(aux)-boxy(aux)),'b^','MarkerFaceColor','b');
    end;
    xlabel('Station number');
%% CMDR NIST recommended symbol for liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title('Secondary Sensor -Uncalibrated Data');
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
    plot(cprs,coxy2-boxy,'ro');
%% CMDR ignore nans they make a mess
    diff_mean=nanmean(coxy2-boxy);
    diff_stddev=nanstd(coxy2-boxy);
    grid on;
    xlabel('Pressure [dbar]');
%% CMDR NIST recommended symbol for liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title(['Secondary Sensor - Mean Diff: ',num2str(diff_mean), ...
		' Std Dev: ',num2str(diff_stddev)]);
    disp('Raw data - Primary Sensor');
    disp(['Mean residual: ',num2str(diff_mean)]);
    disp(['Standard Deviation: ',num2str(diff_stddev)]);
else
    disp(['Secondary oxygen concentration not available: oxy_match.ctdoxy2']);
end; % if o2cflag
    disp('Hit any key to continue');
    pause;

%%%%%%%%%%% CMDR test for second sensors
if t2flag,
    %% Check the temperature data
    %%
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
    diff_temp=ctem1-ctem2;
    plot(stn, diff_temp,'*');
    hold on;
    aux=find(cprs>=1000);
    plot(stn(aux),diff_temp(aux),'r*');
    legend('all','below 1000');
    plot([st1 st2],[mean(diff_temp) mean(diff_temp)],'r');
    plot([st1 st2],[    mean(diff_temp)+3*std(diff_temp) ...
			mean(diff_temp)+3*std(diff_temp)],'r--');
    plot([st1 st2],[	mean(diff_temp)-3*std(diff_temp) ...
			mean(diff_temp)-3*std(diff_temp)],'r--');
    grid on;
    xlabel('Station');
    ylabel('Temperature (^oC)');
    title('Differences between temperature sensors');
    %% some statistics
    fprintf(1,['********************************** \n']);
    disp('Temperature section');
    fprintf(1,['********************************** \n']);
    disp('All data');
    fprintf(1,['Mean difference: ', num2str(mean(diff_temp)),'\n']);
    fprintf(1,['Std. Dev : ', num2str(std(diff_temp)),'\n']);
    disp('Data below 1000 dbar');
    fprintf(1,['Mean difference: ', num2str(mean(diff_temp(aux))),'\n']);
    fprintf(1,['Std. Dev : ', num2str(std(diff_temp(aux))),'\n']);
else
    disp(['Secondary temperature not available: oxy_match.ctdtmp2']);
%% CMDR set temperature diffs to zero
    diff_temp=ctem1.*0;
end; 	% if t2flag
    disp('Hit any key to continue...');
    pause;

    dt_max=mean(diff_temp)+3*std(diff_temp);
    dt_min=mean(diff_temp)-3*std(diff_temp);
    good_data=find(diff_temp>=dt_min & diff_temp<=dt_max);
%% CMDR if some of our sensors are missing we can skip questions about
%% whether to use primary or secondary sensors
  % Choose the sensor to calibrate pri/sec
if t2flag & o2vflag & o2cflag,
    disp(' 0 = use the primary sensor (DEFAULT)' )
    disp(' 1 = use the secondary sensor' )
    pri_sec=input('Do you want to use the primary or the secondary sensor: ');
    while (isempty(pri_sec)==1 | pri_sec<0 | pri_sec>1),
	 disp('Please answer 0 (zero) for no or 1 (one) for yes ');
	 pri_sec=input(['Do you want to use the primary or ' ...
		'the secondary sensor: ']);
    end;
else
	pri_sec=0;
end;
    stn=stn(good_data);
    nis=nis(good_data);
    o2_bot_num=o2_bot_num(good_data);
    cprs=cprs(good_data);
    boxy=boxy(good_data);

    if pri_sec==0
     ctem=ctem1(good_data);
     csal=csal1(good_data);
     coxv=coxv1(good_data);
     coxy=coxy1(good_data);
     ccon=ccon1(good_data);
    sensor='Primary';
   else
     ctem=ctem2(good_data);
     csal=csal2(good_data);
     coxv=coxv2(good_data);
     coxy=coxy2(good_data);
     ccon=ccon2(good_data);
    sensor='Secondary';
end;
    fprintf(1,['After the temperature inspection ... \n'])
%% CMDR give some more information
    fprintf(1,[num2str(length(coxy_ori)-length(good_data)), ' data points ' ...
		'out of ' num2str(length(coxy_ori)) ' were removed  \n']);
    fprintf(1, ...
	[num2str(((length(coxy_ori)-length(good_data))/length(coxy_ori))*100),...
		'%% of the total points \n']);
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
    plot(stn,coxy-boxy,'r+');
    grid on;
    hold on;
    for i=st1:1:st2,
    aux=find(stn==i);
%%%%%% CMDR use nanmean to prevent any NaNs spoiling the result
    plot(i,nanmean(coxy(aux)-boxy(aux)),'b^','MarkerFaceColor','b');
    end;
    xlabel('Station number');
%% CMDR NIST recommended symbol liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title([sensor,' Sensor - after temperature based removal']);
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
%% here we do coxy-boxy and ....
    plot(cprs,coxy-boxy,'ro');
    mdiff=mean(coxy-boxy);
    sdiff=std(coxy-boxy);
    grid on;
	hold on;
    xlabel('Pressure [dbar]');
%% CMDR NIST recommended symbol liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title([sensor,' Sensor - Mean Diff: ',num2str(mdiff), ...
	' Std Dev: ',num2str(sdiff)]);
    disp('Hit any key to continue');
    pause;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%% Initial Quality Control of the data %%%%%%%%%%%%%%%%
%% ...here we do boxy-coxy!!! Is this a new figure??
%% I guess it is a new figure! Make a new figure.
fcnt=fcnt+1;
	fh(fcnt)=figure;
plot(boxy-coxy,'k*');
d=boxy-coxy;
desv=nan_std(d);
%%%%%%%% CMDR use absolute value function
% aux_bad=find(d > nan_mean(d)+3*desv | d < nan_mean(d)-3*desv);
aux_bad=find(abs(d) > nan_mean(d)+3*desv);
%%%%%%%%
disp(['There are ',num2str(length(aux_bad)),' bad points.']);
boxy_aux=boxy;
boxy_aux(aux_bad)=NaN;
hold on;
plot(boxy_aux-coxy,'r*');
a=unique(o2_bot_num);
for i=1:1:length(a),
	aux=find(o2_bot_num==a(i));
	out(i,1)=nanmean(d(aux));
	out(i,2)=nanstd(d(aux));
	out(i,3)=length(aux);
end;
out(:,4)=a;
aux_bad=1;
while length(aux_bad)>0,
   d_mean=nanmean(out(:,1));
   d_std=nanstd(out(:,1));
%%%%%%%%% CMDR use absolute value function
%    aux_bad=find(out(:,1)>=d_mean+3*d_std);
%    out(aux_bad,1)=NaN;
%    out(aux_bad,2)=NaN;
%    aux_bad=find(out(:,1)<=d_mean-3*d_std);
%    out(aux_bad,1)=NaN;
%    out(aux_bad,2)=NaN;
    aux_bad=find(abs(out(:,1))>=d_mean+3*d_std);
    out(aux_bad,1:2)=NaN.*out(aux_bad,1:2);
%%%%%%%%%
end;
aux=find(isnan(out(:,1))==1);
aux1=find(out(aux,3)>=floor(mean(out(:,3))));
bad_bot=out(aux(aux1),4);

%%%%%% CMDR plot the bad bottles
%%%%%% I'm not sure how to find out which are the bad bottles. wait till we
%%%%%% get some, then test some more.
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%% plot normalized differences in standard deviation units
	normd=(d-nanmean(d))./desv;
	plot(o2_bot_num,normd,'.');
	hold on;
 	plot(o2_bot_num(aux_bad),normd(aux_bad),'r+');
%	plot(o2_bot_num(aux_bad),d(aux_bad),'r+');
%	plot([min(o2_bot_num) max(o2_bot_num)],[nan_mean(d) nan_mean(d)],'r--');
%	plot([min(o2_bot_num);max(o2_bot_num)]*ones(1,2), ...
%		nan_mean(d)+desv.*ones(2,1)*[-1 1],'r:');
	title(['Bad bottles']);
	xlabel(['Sample bottle number']);
	ylabel(['Normalized Bottle-CTD Oxygen (standard deviation)']);
%%%%%%
%%%%%%
%%%%%%
%%%%%%

disp('***************************************************** ')
disp('This program suggests to remove the sample bottles:')
disp(num2str(bad_bot))
disp('***************************************************** ')
%
%Removing the bad sample bottles
counter=1;
flag_bottle=0;
bad_btl(counter)=input(['Enter oxygen bottle sample number to remove or ', ...
	'0 to skip: ']);
 while bad_btl(counter)~=0,
 counter=counter+1;
 bad_btl(counter)=input(['Enter oxygen bottle sample number to remove or ', ...
	'0 to skip: ']);
 end;
 if (counter > 1),
 flag_bottle=1;
 bad_bottle=bad_btl(1:length(bad_btl)-1);
 disp('Sample Bottles marked as bad');
 bad_bottle
 end;
 if (flag_bottle==1),
    for b=1:length(bad_bottle),
    bad_index=find(o2_bot_num==bad_bottle(b));
    boxy(bad_index)=NaN;
    end;
end;
jj_nan=find(isnan(boxy)==0);
stn=stn(jj_nan);
nis=nis(jj_nan);
o2_bot_num=o2_bot_num(jj_nan);
cprs=cprs(jj_nan);
ctem=ctem(jj_nan);
ccon=ccon(jj_nan);
csal=csal(jj_nan);
coxv=coxv(jj_nan);
coxy=coxy(jj_nan);
boxy=boxy(jj_nan);
%%%%%%
%Calculating the oxygen saturation

%%%%%%%%%%% CMDR use seawater O2 saturation function
[msize,nsize] = size(boxy);
oxy=zeros(msize,nsize);
oxsat=zeros(msize,nsize);
%
oxsat=sw_satO2(ctem,csal);
%%%%%%%%%%%%%%
keyboard
% choosing the fit
disp('Possible equations to fit: ');
disp('  0 = None - no fit for the station number (DEFAULT)');
disp('  1 = First Order:  linear fit [p1*sta] ');
method=input('What kind of fit for the station number: ');
if (isempty(method)==1),
    method=0;
end;
if (method==0),
[new_oxy,Soc,Voffset,A,B,C,E]=...
    	oxy_fit_new_equation(coxv,oxsat,ctem,cprs,stn,boxy,Soc,Voffset,A,B,C,E);
end
if (method==1)
p1=0;
    [new_oxy,Soc,Voffset,A,B,C,E,p1]=...
    	oxy_fit_new_equation_linear(coxv,oxsat,ctem,cprs,stn,boxy,Soc,Voffset,A,B,C,E,p1);
end
%%%% plotting the result
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
    plot(stn,new_oxy-boxy,'r+');
    grid on;
    hold on;
    for i=st1:1:st2
    aux=find(stn==i);
    plot(i,mean(new_oxy(aux)-boxy(aux)),'b^','MarkerFaceColor','b');
    end;
    xlabel('Station number');
%% CMDR NIST recommended symbol liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title([sensor,' Sensor - Final Result']);
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
    plot(cprs,new_oxy-boxy,'ro');
    mdiff=mean(new_oxy-boxy);
    sdiff=std(new_oxy-boxy);
    grid on;
    xlabel('Pressure [dbar]');
%% CMDR NIST recommended symbol liter is capital L
    ylabel('CTD Oxygen  - Bottle Oxygen (mL/L)');
    title([sensor,' Sensor - Mean Diff: ',num2str(mdiff), ...
		'Std Dev: ',num2str(sdiff)]);
    disp('Hit any key to continue');
pause;
% calculating final statistics
sub=length(new_oxy);
d=boxy-new_oxy;
m=mean(d);
m_abs=mean(abs(d)./boxy)*100;
s=std(d);

% Displaying the new values for the coefficients
disp([sensor,' Sensor']);
disp(['Number of points used   ',int2str(sub)]);
disp(['Total number of points  ',int2str(tot)]);
disp(['Percent of points used  ',num2str(100*sub/tot,4)]);
disp(['fit standard deviation  ',num2str(s,4)]);
disp(['fit mean residual  ',num2str(m,4)]);
disp(['fit mean absolute residual  ',num2str(mean(abs(d)),4)]);
disp(['fit mean absolute residual (%)  ',num2str((mean(abs(d)./boxy)*100),5)]);
disp(['New Soc = ', num2str(Soc)]);
disp(['New Voffset = ', num2str(Voffset)]);
disp(['New A = ', num2str(A)]);
disp(['New B = ', num2str(B)]);
disp(['New C = ', num2str(C)]);
disp(['New E = ', num2str(E)]);
oxy_coeff(1)=Soc;
oxy_coeff(2)=Voffset;
oxy_coeff(3)=A;
oxy_coeff(4)=B;
oxy_coeff(5)=C;
oxy_coeff(6)=E;
if (method==1),
	disp(['New p1 = ', num2str(p1)]);
	oxy_coeff(7)=p1;
end
% Another relevant plot
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
plot(stn,ox_units(boxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk')- ...
	ox_units(new_oxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk'),'bo');
grid on;
title(['Dissolved Oxygen Calibration - Final Result']);
ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
xlabel('Station number');
diff_o2=ox_units(boxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk')- ...
	ox_units(new_oxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk');
diff_o2_pct=diff_o2./ox_units(boxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk');
mumol=nan_mean(diff_o2);
sumol=nan_std(diff_o2);
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
plot(cprs,ox_units(boxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk')- ...
	ox_units(new_oxy,csal,sw_ptmp(ctem,csal,cprs,0),'mll2mmk'),'bo');
grid on;
title([' mean residual = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
ylabel('Bottle - Primary CTD Oxygen (\mu mol/kg)');
xlabel('Pressure');
hold on;
disp('Hit any key to continue');
pause;

aux_contour=[stn  cprs diff_o2 diff_o2_pct*100];
pos_idx_shallow=find(aux_contour(:,4)>=0 & aux_contour(:,2)<=1000);
neg_idx_shallow=find(aux_contour(:,4)<0 & aux_contour(:,2)<=1000);
pos_idx_deep=find(aux_contour(:,4)>=0 & aux_contour(:,2)>1000);
neg_idx_deep=find(aux_contour(:,4)<0 & aux_contour(:,2)>1000);
%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
plot(aux_contour(pos_idx_shallow,1),aux_contour(pos_idx_shallow,2), ...
	'r+','MarkerSize',1)
ax1=gca;
set(ax1,'ydir','reverse')
hold on
plot(aux_contour(neg_idx_shallow,1),aux_contour(neg_idx_shallow,2), ...
	'k+','MarkerSize',1)
text(aux_contour(pos_idx_shallow,1),aux_contour(pos_idx_shallow,2), ...
	num2str(aux_contour(pos_idx_shallow,4),1),'Color','r')
text(aux_contour(neg_idx_shallow,1),aux_contour(neg_idx_shallow,2), ...
	num2str(aux_contour(neg_idx_shallow,4),1),'Color','k')
title('Final Result Upper 1000 dbar : Percentage Error (%)')
xlabel('Station')
ylabel('Pressure (dbar)')

%%%%%% CMDR collect figure handles
    fcnt=fcnt+1;
	fh(fcnt)=figure;
%%%%%%
plot(aux_contour(pos_idx_deep,1),aux_contour(pos_idx_deep,2), ...
	'r+','MarkerSize',1)
ax1=gca;
set(ax1,'ydir','reverse')
hold on
plot(aux_contour(neg_idx_deep,1),aux_contour(neg_idx_deep,2), ...
	'k+','MarkerSize',1)
text(aux_contour(pos_idx_deep,1),aux_contour(pos_idx_deep,2), ...
	num2str(aux_contour(pos_idx_deep,4),1),'Color','r')
text(aux_contour(neg_idx_deep,1),aux_contour(neg_idx_deep,2), ...
	num2str(aux_contour(neg_idx_deep,4),1),'Color','k')
title('Final Result Below 1000 dbar : Percentage Error (%)')
xlabel('Station')
ylabel('Pressure (dbar)')

%% CMDR save the coefficients to the database file before we return
%% 	we don't have to do anything more than save it, because we loaded
%%	data from it earlier.
% save(db_file,'oxy_coeff','-append');
%%

return;

