function [beta0]=cond_calfit(cruiseid)
% cond_cal - Conductivity calibration calling routine.
%
% CTD Calibration toolbox.
%
% USAGE:
% [cond_coeffs]=cond_cal(cruiseid);
%
% INPUT: cruiseid: Cruise identification name.
%   
% OUTPUT: ...
%          
% DESCRIPTION: ...  
% 
% AUTHORS:
%			Carlos Fonseca 
%         UM/CIMAS
%         Derrick Snowden
%         NOAA/AOML/PhOD
%         Tue May 11 11:23:37 EDT 2004
%
% CHANGELOG: 
%   23-Jul-2004, Version 1.0
%		* Initial version. 
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% 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
st1=input('Enter first station number: ');
st2=input('Enter last station number: ');
 % Choose the sensor to calibrate pri/sec
disp(' 1 = use the primary sensor (DEFAULT)' )
disp(' 2 = 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>2)
        disp('Please answer 1 for primary or  2 for secondary ')
        pri_sec=input('Do you want to use the primary or the secondary sensor: ');
    end 

%choose a depth for quality control
%disp('You have an option to color code the output selecting a depth')
disp('This depth will be also used in the criteria to quality control the data' )
col_code=input('Do you want select a depth? (0=no 1=yes): ');
    while (isempty(col_code)==1 | col_code<0 | col_code>1)
        disp('Please answer 0 (zero) for no or  1(one) for yes')
        col_code=input('Do you want select a depth? (0=no 1=yes): ');
    end;
    if col_code==0
        q_depth=0;
    end;
    if col_code==1
        q_depth=input('Enter the depth for color coding the output: ');
    end;

% loading the data
if pri_sec==1
 ii=find(cond_match.stnnbr>=st1 & cond_match.stnnbr<=st2  & cond_match.quality1==2 & cond_match.quality2==2);
 length(ii)   
 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);
    tmp2=cond_match.ctdtmp2(ii);
    sensor='Primary';
else
    ii=find(cond_match.stnnbr>=st1 & cond_match.stnnbr<=st2  & cond_match.quality1==2 & cond_match.quality2==2);
    length(ii) 
    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.ctdtmp2(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    
    
    % preliminary qc of the data using the temperature
disp('*********************************')
disp('Preliminary inspection of the data based on temperature')
disp('*********************************')    
    
    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']) 
    
    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([sensor,' 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)]);
    disp('Hit ant key to continue')
    pause
    %% Check the temperature data
    diff_temp=tmp1-tmp2;
    figure
    plot(ctdsta, diff_temp,'*')
    hold on
    aux=find(ctdprs>=q_depth);
    plot(ctdsta(aux),diff_temp(aux),'r*')
    legend(['all','below',num2str(q_depth)])
    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 setion')
    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',num2str(q_depth),'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 accordingly with 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'])
    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'])  
    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)])  

    disp('Hit ant key to continue')
    pause
 
%%%% Start aplying corrections to the data    
aux=find(ctdprs>=q_depth);
%keyboard
w(1:length(ctdcon),1)=1;
w(aux)=7;
w=w/mean(w);
btlcon_w=sqrt(w).*btlcon;
x=[ctdcon ctdprs ctdsta w];
beta=[1 0 -0.02 1e-5]
beta0=non_lin_fit(x,btlcon_w,'cpcor_fit_w',beta);
ctdcon=(beta0(1)*ctdcon)+(beta0(2).*ctdsta)+beta0(3)+beta0(4)*ctdprs;
%
 
 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
 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)])  
c3515 = sw_c3515;
ctdsal=sw_salt(ctdcon./c3515,ctdtmp,ctdprs);
btlsal=sw_salt(btlcon./c3515,ctdtmp,ctdprs);
figure
plot(ctdsta,ctdsal-btlsal,'b+')
grid on;
hold on
for i=st1:1:st2
 aux=find(ctdsta==i);
 plot(i,mean(ctdsal(aux)-btlsal(aux)),'r^','MarkerFaceColor','r')
end
xlabel('Station number')
ylabel('Bottle Salinity- CTD Salinity (psu)')
title([sensor,' Sensor - Pressure Correction']) 
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)])  

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']) 

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)])  

figure
aux2=find(ctdprs>=3000);
plot(ctdsta(aux2),ctdsal(aux2)-btlsal(aux2),'b+')
grid on;
xlabel('Station number')
ylabel('Bottle Salinity- CTD Salinity (psu)')
title([sensor,' Sensor - Pressure Correction']) 

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)])  



