function ctd_cal=apply_calibration(cruiseid,coeff)
% APPLY_CALIBRATION - apply oxygen and conductivity calibration
%
% CTD Calibration toolbox.
%
% USAGE:
% cond_match=apply_calibration(cruiseid,coeff);
%
% INPUT: cruiseid: Cruise identification name.
%        coeff: matrix with coefficients from conductivity calibration
%   
% OUTPUT: .... 
%          
% DESCRIPTION: ...  
% 

% AUTHORS:
%			Carlos Fonseca 
%         UM/CIMAS
%         Derrick Snowden
%         NOAA/AOML/PhOD
%         Tue May 11 11:23:37 EDT 2004
%
% 
% VERSION: 
%	$Revision: 1.7 $
%	$Date: 2008/02/18 19:30:53 $
%
% CHANGELOG: 
%   07-Oct-2004, Version 1.0
%		* Initial version. 
%  2007/12/17 - CMDR
%	* read cruise id from the file in the current directory
%	* use a function for oxygen saturation calculation
%  2008/01/11 CMDR
%	* use the oxygen saturation function from the seawater toolbox
% 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

tcor = 0.0053903; 
pcor =  0.00013621;
Boc=0.2106;
Soc=0.31655; 
Voffset=-1.1596;
cal=0;
gama=0;
po1=0.00049081;
po2=0;
po3=0;
po4=0;
po5=0;
po6=0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
load_model_types;
%
flag_message1=0;
flag_message2=0;
flag_conv1=0;
flag_conv2=0;
%
%%%%%%%%%%% CMDR read cruiseid if we can
if nargin==0, cruiseid=getcruiseid; end;
if isempty(cruiseid), cruiseid=getcruiseid; end;
%%%%%%%%%%%
% 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('cnv_dlp');

if isempty(n)==1
	error(['cnv_dlp not found. Run deloop.m ']);
	return;
end
disp(['****************************************************************'...
	'**********************'])
disp(['This is the program that applies the calibration coefficients to ' ...
	'the .cnv delooped data.'])
disp(['Edit the values for the oxygen coefficients in the begining of ' ...
	'this m-file.'])
disp('Rename and save in your cruise folder.')
disp(['****************************************************************'...
	'**********************'])

disp(['****************************************************************'...
	'**********************'])
disp(' 0 = use the primary sensor (DEFAULT)' )
disp(' 1 = use the secondary sensor' )
pri_sec=input(['Do you want apply the calibration to the primary or to '...
	'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
disp(['****************************************************************' ...
	'***********'])
fit_choice=input(['Enter the number of the conductivity fit that you want '...
	'to use: ']);
disp(' 0 = none (Seabird Equation)' )
disp(' 1 = exponential' )
disp(' 2 = polynomial')
oxy_choice=input('What method did you choose for the oxygen calibration? ');
while (isempty(oxy_choice)==1 | oxy_choice<0 | oxy_choice>3)
	disp('Invalid oxygen method')
	oxy_choice=input(['What method did you choose for the oxygen ' ...
		'calibration? ']);
end

staa=input('Input first station: ');
stab=input('Input last station: ');
%%%%%%%%%%%%%%%%%%%%
%max_cnv=max(size(cnv));
%for i=1:max_cnv;
pscount=-1; 
  for j=1:1:length(cnv_dlp)
       if staa==cnv_dlp(j).station
       sta_min_idx=j;
       end
       if stab==cnv_dlp(j).station
       sta_max_idx=j;
       end
  end
for i=sta_min_idx:1:sta_max_idx;
    pscount=pscount+1;
    press_aux=cnv_dlp(i).pressure;
    if pri_sec==1
     cond_aux=cnv_dlp(i).condpri*10;
     aux_dcdp=gradient(cond_aux);
     temp_aux=cnv_dlp(i).temppri;
     oxyvolt_aux=cnv_dlp(i).oxvpri;
     oxyvolt_aux=smooth(press_aux,oxyvolt_aux, 0.02,'loess');
    else
     cond_aux=cnv_dlp(i).condsec*10;
     aux_dcdp=gradient(cond_aux);
     temp_aux=cnv_dlp(i).tempsec;
     oxyvolt_aux=cnv_dlp(i).oxvsec;
     oxyvolt_aux=smooth(press_aux,oxyvolt_aux, 0.02,'loess');
    end
  % correcting the conductivity
        if (fit_choice<8)
         ldx=['X=s',num2str(fit_choice-1), ...
		'(cond_aux,repmat(sta_idx,length(cond_aux),1),', ...
		'press_aux,aux_dcdp);'];
        elseif(fit_choice>7 &fit_choice<15)
         ldx=['X=p',num2str(fit_choice-8), ...
		'(cond_aux,repmat(sta_idx,length(cond_aux),1),', ...
		'press_aux,aux_dcdp);'];  
        elseif(fit_choice>14 & fit_choice<22)
         ldx=['X=dcdp',num2str(fit_choice-15), ...
		'(cond_aux,repmat(sta_idx,length(cond_aux),1),', ...
		'press_aux,aux_dcdp);'];
        else
         ldx=['X=pdcdp',num2str(fit_choice-22), ...
		'(cond_aux,repmat(sta_idx,length(cond_aux),1),', ...
		'press_aux,aux_dcdp);'];  
        end
        eval(ldx);
        c=coeff(fit_choice,2:1+length(X(1,:)));
        condc=X*c';
      %condc=weim(11,'hann',cond_cor)';
        sal_aux=sw_salt(condc/sw_c3515,temp_aux,press_aux); 
   
   % correcting the oxygen (oxy fit 3.01 std dev)
   % calculating oxygen saturation
%%%%%%%%%%% CMDR use seawater O2 saturation function
	oxsat=sw_satO2(temp_aux,sal_aux);
%%%%%%%%%%% 
   % calibration coefficients for the oxygen
  sta=repmat(sta_idx,length(oxyvolt_aux),1); 
  if (oxy_choice==0)
     new_oxy=((Soc.*(oxyvolt_aux+Voffset))+Boc.*exp(-0.03.*temp_aux)+...
           (cal.*sta_idx)).*exp(tcor.*temp_aux).*oxsat.*exp(pcor.*press_aux);
  end
  if (oxy_choice==1)
    new_oxy=((Soc.*(oxyvolt_aux+Voffset))+Boc.*exp(-0.03.*temp_aux)+...
		cal.*exp(gama*sta))...
             .*exp(tcor.*temp_aux).*oxsat.*exp(pcor.*press_aux);  
      
  end
  if (oxy_choice==2)
  new_oxy=((Soc*(oxyvolt_aux+Voffset))+Boc*exp(-0.03.*temp_aux)+ ...
		(po1*sta+(po2*sta.^2)+...
           po3*sta.^3)+(po4*sta.^4)+(po5*sta.^5)+(po6*sta.^6))...
           .*exp(tcor.*temp_aux).*oxsat.*exp(pcor.*press_aux);
  end     
  % saving the variables
  ctd_cal(i).pressure(:,1)=press_aux;
  ctd_cal(i).depth(:,1)=sw_dpth(press_aux,cnv(i).lat);
  ctd_cal(i).temp(:,1)=temp_aux;
  ctd_cal(i).pottemp(:,1)=sw_ptmp(sal_aux,temp_aux,press_aux,0);
  ctd_cal(i).salt(:,1)=sal_aux;
  ctd_cal(i).cond(:,1)=condc;
  % converting the oxygen to micromol/kg
  oxy_aux=ox_units(new_oxy,sal_aux, ...
	sw_ptmp(sal_aux,temp_aux,press_aux,0),'mll2mmk');
  ctd_cal(i).oxyumol(:,1)=oxy_aux;
  ctd_cal(i).oxymll(:,1)=new_oxy;
  % calculating the dynamic height
  svan = sw_svan(sal_aux,temp_aux,press_aux);
  clear dh
  for k=1:length(svan); 
    dh(k)=sum(svan(1:k))*1000;
  end
  ctd_cal(i).dynheight(:,1)=dh;     
  % calculating sigma-theta
  pden = sw_pden(sal_aux,temp_aux,press_aux,0)-1000;
  ctd_cal(i).sigma0(:,1)=pden; 
 %%%%%%%%% 
  %  selecting the bottles of the station

 blist=find(cond_match.stnnbr==sta);
 olist=find(oxy_match.stnnbr==sta);
 

  %plotting the data
  
  laty=cnv_dlp(i).latitude;
  lonx=cnv_dlp(i).longitude;
  botstr=num2str(floor(cond_match.ctdprs(blist)));
  hr=floor(cnv_dlp(i).gtime(4));
  mn=floor(cnv_dlp(i).gtime(5));
  sc=floor(cnv_dlp(i).gtime(6));
  date0=datenum(cnv_dlp(i).gtime(1),cnv_dlp(i).gtime(2), ...
		cnv_dlp(i).gtime(3),hr,mn,sc);
  dstr=datestr(date0,0);
  dstr0=dstr(1:17);

  latstr0=sprintf('%7.3f',laty);
  lonstr0=sprintf('%7.3f',abs(lonx));


  ctd0=sprintf('%s','AMMA01 NOAA Ship Ronald H. Brown');
  sname=strcat('CTD',sprintf('%03d',i));
  ctd4=['Latitude',latstr0, ' N'];
  ctd1=[sprintf('%s','CTD Station ',int2str(i)),'  (' , sname, ')'];
  ctd5=['Longitude',lonstr0, ' W'];
  ctd3=sprintf('%s',dstr0,' Z');


figure
hl1=plot(ctd_cal(i).pottemp,ctd_cal(i).pressure,'k');
min_y=0;
max_y=max(ctd_cal(i).pressure+100);
axis([0 30 min_y max_y]);
set(gca,'position',[0.1 0.1 0.45 0.8]);
set(hl1,'linewidth',2);
ax1=gca;
grid on;

set(ax1,'XMinorTick','on','ydir','reverse', 'box','on', ...
	'xcolor',get(hl1,'color'))
xlabel('Potential Temperature (Deg. C)');
ylabel('Pressure (dbar)');
text(5,-0.120*max_y,[ctd0],'FontWeight','bold');  %Cruise Number
text(9,-0.09*max_y,[ctd1],'FontWeight','bold');  %CTD Station Number
text(9,-0.03*max_y,[ctd3],'FontWeight','bold');   %Time
text(5,-0.060*max_y,[ctd4],'FontWeight','bold');  %Latitude
text(16,-0.060*max_y,[ctd5],'FontWeight','bold'); %Longitude
hold on


for oo=1:length(blist);
   bhan(oo)=text(30.1,cond_match.ctdprs(blist(oo)),botstr(oo,:), ...
	'FontSize',6); %bottle pressures
end

[hl2,ax2] = floatAxisX(ctd_cal(i).salt,ctd_cal(i).pressure,'k', ...
	'Salinity (PSS-78)',[33.7 37.3 min_y max_y]);
%set(ax2,'ydir','reverse');
set(ax2,'ydir','reverse','ylim',[min_y max_y],'xlim',[33.7 37.3]);
set(hl2,'linewidth',0.5);

[hl3,ax3] = floatAxisX(ctd_cal(i).oxyumol,ctd_cal(i).pressure,'k-.', ...
	'Oxygen (umol/kg)',[75 280 min_y max_y]);
%set(ax3,'ydir','reverse');
set(ax3,'ydir','reverse','ylim',[min_y max_y],'xlim',[75 280]);
set(hl3,'linewidth',2,'color',[0.4 0.4 0.4]);

[h14,ax4] = floatAxisX(ctd_cal(i).sigma0,ctd_cal(i).pressure,'k--',...
	'Sigma-Theta (kg/m3)',[22 28 min_y max_y]);
%set(ax4,'ydir','reverse');
set(ax4,'ydir','reverse','ylim',[min_y max_y],'xlim',[22 28]);
set(h14,'linewidth',0.5);

tax=subplot('position',[.45 .22 .08 .01]);
tpx=[0:.01:1];
tpy=ones(size(tpx)).*.5;
ttp=plot(tpx,tpy,'k');
%set(ttp,'markersize',2);
set(ttp,'linewidth',2);
set(gca,'xtick',[],'xcolor','w');
set(gca,'ytick',[],'ycolor','w');
set(tax,'xcolor','w','box','off','visible','off');
set(tax,'ycolor','w','box','off');

sax=subplot('position',[.45 .17 .08 .01]);
tts=plot([0 1],[.5 .5],'k');
set(tts,'linewidth',0.5);
set(gca,'xtick',[],'xcolor','w');
set(gca,'ytick',[],'ycolor','w');
set(sax,'xcolor','w');
set(sax,'ycolor','w','box','off','visible','off');

oax=subplot('position',[.45 .12 .08 .01]);
tto=plot([0 1],[.5 .5],'k-.');
%set(tto,'markersize',2);
set(tto,'linewidth',2,'color',[0.4 0.4 0.4]);
set(gca,'xtick',[],'xcolor','w');
set(gca,'ytick',[],'ycolor','w');
set(oax,'xcolor','w');
set(oax,'ycolor','w','box','off','visible','off');

sgax=subplot('position',[.45 .07 .08 .01]);
plot([0 1],[.5 .5],'k--');
set(gca,'xtick',[],'xcolor','w');
set(gca,'ytick',[],'ycolor','w');
set(sgax,'xcolor','w','box','off');
set(sgax,'ycolor','w','visible','off');


ts=subplot('position',[.67 .70 .22 .22]);
plot(ctd_cal(i).salt,ctd_cal(i).pottemp,'k');
xlabel('Salinity (PSS-78)');ylabel('Potential Temperature (deg C)');
grid on
hold on
plot(cond_match.btlsal(blist),sw_ptmp( ...
	cond_match.btlsal(blist), ...
	cond_match.ctdtmp1(blist), ...
	cond_match.ctdprs(blist),0),'k.');
set(ts,'xlim',[33.7 37.3],'ylim',[0 30]);

to=subplot('position',[.67 .41 .22 .22]);
plot(ctd_cal(i).oxyumol,ctd_cal(i).pottemp,'k');
xlabel('Oxygen (umol/kg)');ylabel('Potential Temperature (deg C)');
grid on
hold on
plot(ox_units(oxy_match.btloxy(olist),oxy_match.ctdsal1(olist), ...
	sw_ptmp( ...
	oxy_match.ctdsal1(olist), ...
	oxy_match.ctdtmp1(olist), ...
	oxy_match.ctdprs(olist),0),'mll2mmk'), ...
	sw_ptmp(oxy_match.ctdsal1(olist),oxy_match.ctdtmp1(olist), ...
	oxy_match.ctdprs(olist),0),'k.');
set(to,'xlim',[75 280],'ylim',[0 30]);


so=subplot('position',[.67 .12 .22 .22]);
plot(ctd_cal(i).oxyumol,ctd_cal(i).salt,'k');
xlabel('Oxygen (umol/kg)');ylabel('Salinity (PSS-78)');
grid on
hold on
for t=1:length(blist)    
aux_samples=find(oxy_match.btlnbr(olist)==cond_match.btlnbr(blist(t)));
if isempty(aux_samples)==1
aux_sample(t)=NaN;   
else
aux_sample(t)=aux_samples(1);
end
end
plot(oxy_match.btloxy(aux_sample),cond_match.btlsal(blist),'k.');
set(so,'xlim',[75 280],'ylim',[33.7 37.3]);
orient tall

%creating table for specific pressures
ctda=['Latitude',latstr0, 'N'];
ctdb=[sprintf('%s','CTD Station ',int2str(sta_idx)),' (' , sname, ')'];
ctdc=['Longitude',lonstr0, 'W'];
ctdd=sprintf('%s',dstr0,'Z');

p=[sprintf('%s\t','Pressure','Temp90','PoTemp90','Salinity','Oxygen', ...
	'DynHt','SigT')];
pp=[sprintf('%s\t','dbar','Deg C','Deg C','PSS-78','umol/kg','m3/kg-1', ...
	'kg/m3')];


table=[1 10:10:30 50 75 100 125 150 200 250 300:100:1500 1750 2000:500:6000];

fid=fopen(['abaco03_',sprintf('%03d',sta_idx),'.tab'],'w');
fprintf(fid,'ABACO-03 NOAA Ship Ronald H. Brown');
fprintf(fid,'%s',' ',ctdb,' ',ctda,' ',ctdc,' ',ctdd);
fprintf(fid,'\n');
fprintf(fid,p);
fprintf(fid,'\n');
fprintf(fid,pp);
fprintf(fid,'\n');

for j=1:length(table);
   for k=1:length(ctd_cal(sta_idx).pressure);
      if(table(j)==ctd_cal(sta_idx).pressure(k));
         fprintf(fid,'%5.0f\t%8.3f\t%8.3f\t%8.3f\t%6.1f\t%8.3f\t%8.3f\n', ...
	ctd_cal(sta_idx).pressure(k),...
         ctd_cal(sta_idx).temp(k),ctd_cal(sta_idx).pottemp(k), ...
	ctd_cal(sta_idx).salt(k),ctd_cal(sta_idx).oxyumol(k),...
         ctd_cal(sta_idx).dynheight(k),ctd_cal(sta_idx).sigma0(k));
      end
   end
end
fprintf(fid,'\n');
dd=[cond_match.ctdprs(blist),cond_match.ctdprs(blist), ...
	cond_match.ctdtmp1(blist),sw_ptmp(cond_match.btlsal(blist), ...
	cond_match.ctdtmp1(blist),cond_match.ctdprs(blist),0),...
        cond_match.btlsal(blist),oxy_match.btloxy(aux_sample)];
bp=[sprintf('%s\t','Pressure','Niskin','Temp90','PoTemp90','Salinity', ...
	'Oxygen')];
bpp=[sprintf('%s\t','dbar','  ','Deg C','Deg C','PSS-78','umol/kg')];
fprintf(fid,bp);
fprintf(fid,'\n');
fprintf(fid,bpp);
fprintf(fid,'\n');

fprintf(fid,'%5.0f\t%5.0f\t%8.3f\t%8.3f\t%8.3f\t%8.1f\n',dd');
fclose(fid);
orient tall
if pscount==0
    print -dps cal_ctd.ps
 else
    print -dps -append cal_ctd.ps
 end
close all 
  
  
  
  
  
  
  
  
  
  
  
%
clear aux_sample
end
    

