function varargout = cond_cal(cruiseid,st1,st2)
% COND_CAL - conductivity calibration calling routine.
%
% CTD Calibration toolbox.
%
% USAGE
% 	[cond_coeff,fh] = cond_cal(cruiseid)
%		or
% 	[offset_out,offset_fit,slope_out,slope_fit,pcor]=cond_cal(cruiseid);
%
% INPUT
%	cruiseid	cruise identification name.
%   
% OUTPUT
%	cond_coeff=[offset_out,offset_fit,slope_out,slope_fit,pcor]
%          
% 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:52:47 $
%	$Id: cond_cal.m,v 1.23 2008/04/14 19:52:47 duncombe Exp duncombe $
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%	* Initial version. 
%
%  2007/12/12 CMDR
%  	* Get cruiseid from a file containing the identifier in the data
%  	directory.
% 	* test for secondary sensors. If they do not exist make a copy of the 
%	primary sensors and use that.
% 	* make input lines easier to read (aesthetic adjustment!)
%	* make functions cond_cal and FS_apply_calibration have the same input
%	and output arguments to facilitate cut and paste
%  2008/01/10 CMDR
% 	* I want the function to output cond_coeff and figure handles.
% 	Incompatible with maintaining backward compatibility. Introduce
% 	varargout to deal with the two situations.
%	* return figure handles
%  2008/01/19 CMDR
%	* put first and last stations in the argument list
%  2008/02/05 CMDR
%	* discourage failure due to invalid input
%  2008/02/18 CMDR
%	* append new conductivity calibration coefficients to the database
%  2008/03/10 CMDR
%	* provide some defaults for st1 and st2; vectorize some loops
%  2008/03/13 CMDR
% 	* chasing bugs and weirdness in pressure correction: complete 
%	reformatting, recommenting, and restyling.  2008/03/19 the weirdness 
% 	in the pressure turned out to be a typo in FS_apply_calibration.m
%  2008/03/16 CMDR
%	* use default for primary sensor
%  2008/04/14 CMDR
%  	* test if quality fields are present. if not remind operator to run
%	compare_cond.m
%

% BUGS:
% 	* Function produces 20 separate figures. We need six per page for
%	the final report.

% function [offset_out,offset_fit,slope_out,slope_fit,pcor]=cond_cal(cruiseid)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%% CMDR get cruiseid from file
if nargin==0, cruiseid=getcruiseid; end;
if isempty(cruiseid), cruiseid=getcruiseid; end;
%%%%%%%%%%%
%%%%%%%%%%% CMDR initialize figure counter
fcnt=0;

% 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('cond_match');

if isempty(n)==1,
	error(['cond_match not found. Run make_cond_sumfile.m ']);
	return;
end;

stations=unique(cond_match.stnnbr);
station_min=min(stations);
station_max=max(stations);
disp('**************************************************************');
disp('The data found in cond_match goes from');
disp(['the station ',num2str(station_min),' to the station ' ...
		num2str(station_max),'.'] );
disp(['The maximum pressure is: ',num2str(max(cond_match.ctdprs))]);
disp('**************************************************************');
%
%type_fit={'calcos0' 'calcos1' 'calcos2' 'calcos3' 'calcos4' 'calcos5' 'calcos6'...
%          'calcop0' 'calcop1' 'calcop2' 'calcop3' 'calcop4' 'calcop5' 'calcop6'...
%          'calcodcdp0' 'calcodcdp1' 'calcodcdp2' 'calcodcdp3' 'calcodcdp4' 'calcodcdp5' 'calcodcdp6'...
%          'calcopdcdp0' 'calcopdcdp1' 'calcopdcdp2' 'calcopdcdp3' 'calcopdcdp4' 'calcopdcdp5' 'calcopdcdp6'};

% 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;
%%%%%%%%%%%%

% Choose the sensor to calibrate pri/sec
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),
while ( 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; % while (isempty
if isempty(pri_sec), pri_sec=0; disp('Using primary sensor'); end;

% loading the data

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


if pri_sec==0,
  ii=find( cond_match.stnnbr>=st1 & cond_match.stnnbr<=st2  & ...
		cond_match.quality1==1);
  ctdcon=cond_match.ctdcon1(ii);
  ctdcon_ori=cond_match.ctdcon1(ii);
  btlcon=sw_c3515*sw_cndr(cond_match.btlsal(ii),cond_match.ctdtmp1(ii),cond_match.ctdprs(ii));
  ctdprs=cond_match.ctdprs(ii);
  ctdsta=cond_match.stnnbr(ii);
  ctdtmp=cond_match.ctdtmp1(ii);
  tmp1=cond_match.ctdtmp1(ii);
  %%%%%%% CMDR test for second sensors
  if ~isnan(cond_match.ctdtmp2),
     tmp2=cond_match.ctdtmp2(ii); 
  else
     disp('There are no secondary sensors. Using primary sensors instead!');
     tmp2=cond_match.ctdtmp1(ii); 
  end;
  %%%%%%% 
  sensor='Primary';
else 
  ii=find(cond_match.stnnbr>=st1 & cond_match.stnnbr<=st2  ...
	& cond_match.quality2==1);
  ctdprs=cond_match.ctdprs(ii);
  ctdsta=cond_match.stnnbr(ii);
  ctdtmp=cond_match.ctdtmp2(ii);
  btlcon=sw_c3515*sw_cndr(cond_match.btlsal(ii),cond_match.ctdtmp12(ii),cond_match.ctdprs(ii));
  ctdcon=cond_match.ctdcon2(ii);
  ctdcon_ori=cond_match.ctdcon2(ii);
  tmp1=cond_match.ctdtmp1(ii);
  tmp2=cond_match.ctdtmp2(ii);
  sensor='Secondary';
end;  % if pri_sec==0,

% preliminary qc of the data using the temperature
disp('*********************************');
disp('Preliminary inspection of the data based on temperature');
disp('*********************************');
    
%% BEGIN FIGURE 1) uncalibrated data
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta,btlcon-ctdcon,'r+');
  grid on;
  hold on;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    plot(i,mean(btlcon(aux)-ctdcon(aux)),'b^','MarkerFaceColor','b');
  end;
  xlabel('Station number');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor -Uncalibrated Data']);
%% END FIGURE
    
%% BEGIN FIGURE 2) vs stn number
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs,btlcon-ctdcon,'ro');
  diff_mean=mean(btlcon-ctdcon);
  diff_stddev=std(btlcon-ctdcon);
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title(['Primary Sensor - Mean Diff: ',num2str(diff_mean), ...
	' Std Dev: ',num2str(diff_stddev)]);
  disp(['Raw data - ',sensor,' Sensor']);
  disp(['Mean residual: ',num2str(diff_mean)]);
  disp(['Standard Deviation: ',num2str(diff_stddev)]);
%% END FIGURE
disp('Hit any key to continue');
pause;

%% Check the temperature data
diff_temp=tmp1-tmp2;

%% BEGIN FIGURE 3) primary-secondary vs station
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta, diff_temp,'*');
  hold on;
  aux=find(ctdprs>=1000);
  plot(ctdsta(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');
%% END FIGURE

%% print 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']);
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);

% Removing bad data according to the temperature
ctdcon=ctdcon(good_data);
ctdcon_aux=ctdcon;
btlcon=btlcon(good_data);
ctdprs=ctdprs(good_data);
ctdsta=ctdsta(good_data);
ctdtmp=ctdtmp(good_data);

fprintf(1,['After the temperature inspection ... \n']);
fprintf(1,['From ',num2str(length(ctdcon_ori)), ' data points  \n']);
fprintf(1,[num2str(length(ctdcon_ori)-length(good_data)), ...
	' data points were removed  \n']);
fprintf(1,[num2str(((length(ctdcon_ori)-length(good_data))./ ...
		length(ctdcon_ori))*100), ...
	' %% of the total points \n']);

%% BEGIN FIGURE 4) after removal of bad points based on temperature data
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta,btlcon-ctdcon,'r+');
  grid on;
  hold on;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    plot(i,mean(btlcon(aux)-ctdcon(aux)),'b^','MarkerFaceColor','b');
  end;
  xlabel('Station number');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - after temperature based removal']);
%% END FIGURE

%% BEGIN FIGURE 5) pressure dependence
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs,btlcon-ctdcon,'ro');
  mdiff=mean(btlcon-ctdcon);
  sdiff=std(btlcon-ctdcon);
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Mean Diff: ',num2str(mdiff),' Std Dev: ', ...
	num2str(sdiff)]);
%% END FIGURE
disp('Hit any key to continue');
pause;

% weird correction % comment in original 
% ctdcon=ctdcon+((-0.015*exp(-ctdprs./600))-0.004);
    
k=1;
for i=st1:1:st2,
  aux=find(ctdsta==i);% & ctdprs>=2000);
  sta(k)=i;
  if isempty(aux)==1,
    %% CMDR test if offset was set, otherwise give it a value
    if exist('offset'), offset(k)=offset(k-1); else offset(k)=NaN; end;
  else    
    aux_off=repmat(1,length(aux),1);
    offset(k)=aux_off\(btlcon(aux)-ctdcon(aux));
    % ctdcon_aux(aux)=ctdcon(aux)+slope(i)*ctdcon(aux);
  end;
  k=k+1;
end; % for i=st1:1:st2

%% BEGIN FIGURE 6) various functions that fit offset to stn no
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(sta,offset,'k*');
  hold on;
  grid on;
  xlabel('Station number');
  ylabel(' Conductivity Offset (mS/cm)');
  title([sensor,' Sensor - Offset Variation']);
% constant
  mean_offset=mean(offset);
  plot(sta,repmat(mean_offset,length(sta),1),'r');
  d0=offset-mean_offset;
  s0=sqrt(sum(d0.^2));
  r0=1-sum(d0.^2)/sum((offset-nanmean(offset)).^2);
  % r0=NaN;
% linear fit
  [p1,s1]=polyfit(sta,offset,1);
  offset_p1=polyval(p1,sta);
  d1=offset-offset_p1;
  r1=1-sum(d1.^2)/sum((offset-mean(offset)).^2);
  plot(sta,offset_p1','g');
% quadratic
  [p2,s2]=polyfit(sta,offset,2);
  offset_p2=polyval(p2,sta);
  d2=offset-offset_p2;
  r2=1-sum(d2.^2)/sum((offset-mean(offset)).^2);
  plot(sta,offset_p2','b');
% cubic
  [p3,s3]=polyfit(sta,offset,3);
  offset_p3=polyval(p3,sta);
  d3=offset-offset_p3;
  r3=1-sum(d3.^2)/sum((offset-mean(offset)).^2);
  plot(sta,offset_p3','y');
%
  legend('offset','constant','linear','quadratic','cubic');
%% END FIGURE

display('Comparisons between fits');
display(    'fit     normr     %variance');
fprintf (1,['constant ', num2str(s0),      ' ', num2str(r0),'\n']);
fprintf (1,['linear ',   num2str(s1.normr),' ', num2str(r1),'\n']);
fprintf (1,['quadratic ',num2str(s2.normr),' ', num2str(r2),'\n']);
fprintf (1,['cubic ',    num2str(s3.normr),' ', num2str(r3),'\n']);
disp('Hit any key to continue');
pause;

%% CMDR dont fail if invalid input
offset_fit=[];
while isempty(offset_fit) | offset_fit>3 | offset_fit<0 | ischar(offset_fit),
  disp(' 0 = constant' );
  disp(' 1 = linear' );
  disp(' 2 = quadratic' );
  disp(' 3 = cubic' );
  offset_fit=input('What offset fit do you prefer? ');
end; % while
if offset_fit==0,
  ctdcon=ctdcon+nanmean(offset);
  offset_fit=nanmean(offset);
  offset_out=nanmean(offset);
else
  if offset_fit==1,
    offset_fit=p1;    
    offset_aux=offset_p1;
  elseif offset_fit==2,
    offset_fit=p2;    
    offset_aux=offset_p2;    
  elseif offset_fit==3,
    offset_fit=p3;    
    offset_aux=offset_p3;    
  end; % if offset_fit==1,
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    idx=find(sta==i);
    ctdcon(aux)=ctdcon(aux)+offset_aux(idx);
  end; % for i=st1:1:st2,
  offset_out=offset_aux;
end; % if offset_fit==0

%% BEGIN FIGURE 7) corrected offset 
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta,ctdcon-btlcon,'b+');
  grid on;
  xlabel('Station number');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Final offset Correction']);
  hold on;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    plot(i,mean(ctdcon(aux)-btlcon(aux)),'r^','MarkerFaceColor','r');
  end; % for i=st1:1:st2,
%% END FIGURE

%% BEGIN FIGURE 8) pressure dependence
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs,ctdcon-btlcon,'bo');
  mdiff=mean(ctdcon-btlcon);
  sdiff=std(ctdcon-btlcon);
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Mean Diff: ',num2str(mdiff),' Std Dev: ', ...
		num2str(sdiff)]);
  disp('Hit any key to continue');
  pause;
  k=1;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    slope(k)=ctdcon(aux)\(btlcon(aux)-ctdcon(aux));
    ctdcon_aux(aux)=ctdcon(aux)+slope(k)*ctdcon(aux);
    k=k+1;
  end; % for i=st1:1:st2,

%% BEGIN FIGURE 9) slope variation vs stn no
fcnt=fcnt+1;
fh(fcnt)=figure;
  sta=st1:1:st2;
  plot(sta,slope,'k*');
  hold on;
  grid on;
  xlabel('Station number');
  ylabel(' Conductivity Slope (mS/cm)');
  title([sensor,' Sensor - Slope Variation']);
% constant
  mean_slope=mean(slope);
  d0=slope-mean_slope;
  s0=sqrt(sum(d0.^2));
  r0=1-sum(d0.^2)/sum((slope-mean(slope)).^2);
  plot(sta,repmat(mean_slope,length(sta),1),'r');
% linear fit
  [p1,s1]=polyfit(sta,slope,1);
  slope_p1=polyval(p1,sta);
  d1=slope-slope_p1;
  r1=1-sum(d1.^2)/sum((slope-mean(slope)).^2);
  plot(sta,slope_p1','g');
% quadratic
  [p2,s2]=polyfit(sta,slope,2);
  slope_p2=polyval(p2,sta);
  d2=slope-slope_p2;
  r2=1-sum(d2.^2)/sum((slope-mean(slope)).^2);
  plot(sta,slope_p2','b');
% cubic
  [p3,s3]=polyfit(sta,slope,3);
  slope_p3=polyval(p3,sta);
  d3=slope-slope_p3;
  r3=1-sum(d3.^2)/sum((slope-mean(slope)).^2);
  plot(sta,slope_p3','y');
  legend('slope','constant','linear','quadratic','cubic');
%% END FIGURE

display('Comparisons between fits');
display('fit     normr     %variance');
fprintf (1,['constant ',num2str(s0),' NaN \n']);
fprintf (1,['linear ',num2str(s1.normr),' ', num2str(r1),'\n']);
fprintf (1,['quadratic ',num2str(s2.normr),' ', num2str(r2),'\n']);
fprintf (1,['cubic ',num2str(s3.normr),' ', num2str(r3),'\n']);
disp('Hit any key to continue');
pause;

%% CMDR dont fail if invalid input
slope_fit=[];
while isempty(slope_fit) | slope_fit>3 | slope_fit<0 | ischar(slope_fit),
  disp(' 0 = constant' );
  disp(' 1 = linear' );
  disp(' 2 = quadratic' );
  disp(' 3 = cubic' );
  slope_fit=input('What slope fit do you prefer? ');
end; % while
if slope_fit==0,
  ctdcon=ctdcon+nanmean(slope)*ctdcon;
  slope_fit=nanmean(slope);
  slope_out=nanmean(slope);
else
  if slope_fit==1,
    slope_fit=p1;    
    slope_aux=slope_p1;
  elseif slope_fit==2,
    slope_fit=p2;
    slope_aux=slope_p2;    
  elseif slope_fit==3,
    slope_fit=p3;    
    slope_aux=slope_p3;    
  end;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    idx=find(sta==i);
    ctdcon(aux)=ctdcon(aux)+slope_aux(idx)*ctdcon(aux);
  end; %  for i=st1:1:st2,
  slope_out=slope_aux;
end;  % if slope_fit==0,

%% BEGIN FIGURE 10) slope correction
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta,ctdcon-btlcon,'b+');
  grid on;
  xlabel('Station number');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Final Slope Correction']);
  hold on;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    plot(i,mean(ctdcon(aux)-btlcon(aux)),'r^','MarkerFaceColor','r');
  end; % for i=st1:1:st2
%% END FIGURE

%% BEGIN FIGURE 11) after slope correction
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs,ctdcon-btlcon,'bo');
  mdiff=mean(ctdcon-btlcon);
  sdiff=std(ctdcon-btlcon);
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Mean Diff: ',num2str(mdiff),' Std Dev: ', ...
	num2str(sdiff)]);
%% END FIGURE
disp('Hit any key to continue');
pause;

%Pressure Correction
k=1;
for i=st1:1:st2,
%% this was restricted in depth. Try without the restriction.
%% aux=find(ctdsta==i & ctdprs>=1000);
  aux=find(ctdsta==i);
  if (isempty(aux)==1 ),
    press_slope(k)=NaN;
  else
    press_slope(k)=ctdprs(aux)\(btlcon(aux)-ctdcon(aux));
  end;
  k=k+1;
end; % for i=st1:1:st2,

%% BEGIN FIGURE 12) pressure correction (deep stations)
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(sta,press_slope,'k*');
  hold on;
  grid on;
  xlabel('Station number');
  ylabel(' Pressure Slope (S/m)');
  title([sensor,' Sensor - Pressure Correction']);
%% END FIGURE
disp('Hit any key to continue');
pause;

%% This aux_p assignment was repeated redundantly below. Why?
aux_p=find(~isnan(press_slope));

%% CMDR dont fail if invalid input
pcor_fit=[];
while isempty(pcor_fit) | pcor_fit>1 | pcor_fit<0 | ischar(pcor_fit),
  disp(' 0 = no (DEFAULT)' );
  disp(' 1 = yes' );
  pcor_fit=input('Do you want to apply a linear fit for pcor? ');
end; % while isempty
if pcor_fit,
  % aux_p=find(~isnan(press_slope));
  [p1,s1]=polyfit(sta(aux_p),press_slope(aux_p),1);
  pcor_p1=polyval(p1,sta);
  plot(sta,pcor_p1,'r');
  disp('Hit any key to continue');
  pause;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    idx=find(sta==i);
    ctdcon(aux)=ctdcon(aux)+pcor_p1(idx)*ctdprs(aux);
  end; % for i=
  pcor=p1;
else
  pcor=0;  
%   % aux_p=find(~isnan(press_slope));
%   if isempty(aux_p),
%     pcor=0;
%   else
%     pcor=trimmean(press_slope(aux_p),10);
%     if isnan(pcor)==1,
%       pcor=-9.57e-8;
%     end; % if isnan
  end; % if isempty(aux_p)
  ctdcon=ctdcon+pcor*ctdprs;
end; % if pcor_fit,

%% BEGIN FIGURE 13) after deep pressure correction?
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta,ctdcon-btlcon,'b+');
  grid on;
  xlabel('Station number');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Pressure Correction']);
  hold on;
  for i=st1:1:st2,
    aux=find(ctdsta==i);
    plot(i,mean(ctdcon(aux)-btlcon(aux)),'r^','MarkerFaceColor','r');
  end;
%% END FIGURE

%% BEGIN FIGURE 14) conductivity differences
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs,ctdcon-btlcon,'bo');
  mdiff=mean(ctdcon-btlcon);
  sdiff=std(ctdcon-btlcon);
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Conductivity  - CTD Conductivity (mS/cm)');
  title([sensor,' Sensor - Mean Diff: ',num2str(mdiff),' Std Dev: ', ...
		num2str(sdiff)]);
%% END FIGURE

ctdsal=sw_salt(ctdcon./sw_c3515,ctdtmp,ctdprs);
btlsal=sw_salt(btlcon./sw_c3515,ctdtmp,ctdprs);

%% BEGIN FIGURE 15) salinity differences
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdsta,ctdsal-btlsal,'b+');
  grid on;
  xlabel('Station number');
  ylabel('Bottle Salinity- CTD Salinity (psu)');
  title([sensor,' Sensor - Pressure Correction']);
%% END FIGURE

%% BEGIN FIGURE 16) showing nominal sensor precision
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs,ctdsal-btlsal,'bo');
  mdiff=mean(ctdsal-btlsal);
  sdiff=std(ctdsal-btlsal);
  hold on;
  plot(ctdprs,repmat(0.002,1,length(ctdprs)),'r-');
  plot(ctdprs,repmat(-0.002,1,length(ctdprs)),'r-');
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Salinity- CTD Salinity (psu)');
  title([sensor,' Sensor - Mean Diff: ',num2str(mdiff),' Std Dev: ', ...
	num2str(sdiff)]);

%% BEGIN FIGURE 17) deep salinities vs station
fcnt=fcnt+1;
fh(fcnt)=figure;
  aux1=find(ctdprs>=500);
  plot(ctdsta(aux1),ctdsal(aux1)-btlsal(aux1),'b+');
  grid on;
  xlabel('Station number');
  ylabel('Bottle Salinity- CTD Salinity (psu)');
  title([sensor,' Sensor - Pressure Correction']);
%% END FIGURE

%% BEGIN FIGURE 18) deep salinities vs pressure
fcnt=fcnt+1;
fh(fcnt)=figure;
  plot(ctdprs(aux1),ctdsal(aux1)-btlsal(aux1),'bo');
  mdiff=mean(ctdsal(aux1)-btlsal(aux1));
  sdiff=std(ctdsal(aux1)-btlsal(aux1));
  hold on;
  plot(ctdprs,repmat(0.002,1,length(ctdprs)),'r-');
  plot(ctdprs,repmat(-0.002,1,length(ctdprs)),'r-');
  grid on;
  xlabel('Pressure [dbar]');
  ylabel('Bottle Salinity- CTD Salinity (psu)');
  title([sensor,' Sensor - Mean Diff: ',num2str(mdiff), ...
	' Std Dev: ',num2str(sdiff)]);
%% END FIGURE

%% test for deeper data, and skip the last two figures if there are none
aux2=find(ctdprs>=3000);
if ~isempty(aux2), 

%% BEGIN FIGURE 19) deeper salinities 
  fcnt=fcnt+1;
  fh(fcnt)=figure;
    plot(ctdsta(aux2),ctdsal(aux2)-btlsal(aux2),'b+');
    grid on;
    xlabel('Station number');
    ylabel('Bottle Salinity- CTD Salinity (psu)');
    title([sensor,' Sensor - Pressure Correction']);
%% END FIGURE

%% BEGIN FIGURE 20) deeper salinities (skip this if there are none)
  fcnt=fcnt+1;
  fh(fcnt)=figure;
    plot(ctdprs(aux2),ctdsal(aux2)-btlsal(aux2),'bo');
    mdiff=mean(ctdsal(aux2)-btlsal(aux2));
    sdiff=std(ctdsal(aux2)-btlsal(aux2));
    hold on;
    plot(ctdprs,repmat(0.002,1,length(ctdprs)),'r-');
    plot(ctdprs,repmat(-0.002,1,length(ctdprs)),'r-');
    grid on;
    xlabel('Pressure [dbar]');
    ylabel('Bottle Salinity- CTD Salinity (psu)');
    title([sensor,' Sensor - Mean Diff: ',num2str(mdiff), ...
	' Std Dev: ',num2str(sdiff)]);
%% END FIGURE

end; % if ~isempty(aux2),

%% CMDR - I want to save the new coefficients to the data base
%% 	We don't have to do anything more than append it, because
%%	we loaded data from it earlier.

cond_coeff={offset_out, offset_fit, slope_out, slope_fit, pcor};
save(db_file,'cond_coeff','-append');
%%

if nargout==0, assignin('caller','ans',cond_coeff); end;
if nargout>=1 & nargout<=2, varargout{1}=cond_coeff; end;
if nargout==2, varargout{2}=fh; end;
if nargout>=3,
	varargout{1}=offset_out;
	varargout{2}=offset_fit;
	varargout{3}=slope_out;
end;
if nargout>=4, varargout{4}=slope_fit; end;
if nargout==5, varargout{5}=pcor; end;

save('calibdata','ctdcon','btlcon');

return;

