function [oxy_coeff,fh]=oxy_new_equation(cruiseid)
% OXY_CAL - oxygen calibration calling routine. 
% 
% CTD Calibration toolbox.
%
% USAGE
% 	[oxy_coeff,fh]=oxy_new_equation(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
%	Uses a non linear fit to calculate new oxygen coefficients of sensor
%	calibration.  Must change the oxygen coefficients inside the program to
%	account for the sensors used on the cruise.
%
% FUNCTIONS USED
%   prelim_oxy_plots
%   ox_units
%   oxy_fit_new_equation
%   oxy_fit_new_equation_linear

% AUTHOR:
%         Jay Hooper
%         UM/CIMAS
%         NOAA/AOML/PhOD
%         Thursday May 17 2012

%close all;
%Load database 
db_file=([cruiseid,'_db.mat']);
load([cruiseid,'_db.mat']);


%Check to see proper directory is present
n=who('oxy_match');
if isempty(n)==1
	error(['oxy_match not found. Run make_oxy_sumfile_jay2.m ']);
	return;
end;


disp('**************************************************************************************')
disp('This is the program that applies the calibration coefficients to the oxy_match data.')
disp('**************************************************************************************')

%%%%%%%%%%%%%%%%%%% Preliminary look at the Oxygen data %%%%%%%%%%%%%%%%%%%
[fh] = prelim_oxy_plots(oxy_match);
disp('Hit Enter to continue to station selection.')
pause
%%%%%%%%%%%%%%%%%%% End Preliminary look at the Oxygen data %%%%%%%%%%%%%%%


disp('*********************************************************************')
disp('Choose desired Station interval to calculate new Oxygen coefficients.')
disp('*********************************************************************')
staa=input('Input first station: ');
stab=input('Input last station: ');

%%%%%%%%%%%%%%%%%%%%%% Good data statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
disp('*******************************************')
disp('Initial statistical evaluation of good data')
disp('*******************************************')
good_data=find(oxy_match.quality1 == 2 & oxy_match.quality2 == 2 ...
                & oxy_match.stnnbr>= staa & oxy_match.stnnbr  <= stab);
disp(['Number of good data points found: ',num2str(length(good_data))])
aux_gd=find( oxy_match.stnnbr>= staa & oxy_match.stnnbr  <= stab);
per_gd=(length(good_data)/length(aux_gd))*100;
disp(['This represents ', num2str(per_gd),'% of the data']); 
%%%%%%%%%%%%%%%%%%%%%%% End good statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%% Begin New equation %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


disp('**************************************************************************************')
disp(' Do you want apply the calibration to the primary or to the secondary sensor?' )
disp(' 1 = use the primary sensor (DEFAULT)' )
disp(' 2 = use the secondary sensor' )
oxy_sensor=input('Oxygen sensor: ');
while (isempty(oxy_sensor)==1 | oxy_sensor<=0 | oxy_sensor>2)
	disp('Please answer 1 for primary or  2 for secondary ')
	oxy_sensor=input('Do you want to use the primary or the secondary sensor: ');
end
    %staa=2;
    %stab=28;
    %oxy_sensor=2;
    
    %Define quality control filter
    ii=find(oxy_match.quality1 == 2 & oxy_match.stnnbr >= staa & ...
            oxy_match.stnnbr <= stab);
        
    disp('*********************************************************') 
    disp('*********************************************************')
    disp('** Did you change the oxygen coefficients in the file? **')
    disp('*********************************************************')
    disp('*********************************************************')
    pause
    
    disp('****************************************************************') 
    disp('****************************************************************')
    disp('** Swapped O2 sensors 2->1 and 1->2 to calculate new equation **')
    disp('****************************************************************')
    disp('****************************************************************')
    pause
    
    if oxy_sensor == 1,
        %creat variables for primary oxy_match
        sensor='Primary';
        stn=oxy_match.stnnbr(ii);
        oxy_volt=oxy_match.ctdoxv2(ii);
        oxy_mll=oxy_match.ctdoxy2(ii);
        temp=oxy_match.ctdtmp1(ii);
        sal=oxy_match.ctdsal1(ii);
        prs=oxy_match.ctdprs(ii);
        boxy=oxy_match.btloxy(ii);
        cond=oxy_match.ctdcon1(ii);
        tot=length(boxy);        
        
        %********************************%
        % Use cnv.ox to get coefficients %
        %********************************%
        
        %Primary oxygen sensor 1348 ab1202 post cal
        fprintf('Primary sensor cal coeff\t')   
        if 0,
        %s/n 1348
        Soc = 0.5514
        Voffset = -0.5114
        Tau20 = 1.34
        A = -3.6914e-3
        B = 1.7009e-5
        C = -3.0484e-6
        E = 0.036
        end;
        if 1,
        %s/n 2082
        Soc = 0.4223
        Voffset = -0.5233
        Tau20 = 1.47
        A = -2.8059e-3
        B = 1.4396e-5
        C = -2.0227e-7
        E = 0.036
        end;
        
    else
        %creat variables for secondary oxy_match
        sensor='Secondary';
        stn=oxy_match.stnnbr(ii);
        oxy_volt=oxy_match.ctdoxv1(ii);
        oxy_mll=oxy_match.ctdoxy1(ii);
        temp=oxy_match.ctdtmp2(ii);
        sal=oxy_match.ctdsal2(ii);
        prs=oxy_match.ctdprs(ii);
        boxy=oxy_match.btloxy(ii);
        cond=oxy_match.ctdcon2(ii);
        tot=length(boxy);
        
        %Secondary oxygen sensor 2082 ab1202 post cal
        fprintf('Secondary sensor cal coeff\t')
        %s/n 2082
        Soc = 0.4223
        Voffset = -0.5233
        Tau20 = 1.47
        A = -2.8059e-3
        B = 1.4396e-5
        C = -2.0227e-7
        E = 0.036
        if 0,
        %s/n 1266
        Soc = 0.5343
        Voffset = -0.5402
        tau20 = 1.23
        A = -3.3869e-3
        B = 8.5060e-5
        C = -1.5406e-6
        E = 0.036
        end;
    end;
    
    %calculate the saturated oxygen
    oxsat=sw_satO2(temp,sal);
    %Convert bottle oxygen to mL/L
    boxyL=ox_units(boxy,sal,sw_ptmp(sal,temp,prs,0),'mmk2mll');
    
    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                    Calculate new Oxygen ml/l                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 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(oxy_volt,oxsat,temp,prs,stn,boxyL,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(oxy_volt,oxsat,temp,prs,stn,boxyL,Soc,Voffset,A,B,C,E,p1);
    end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                End Calculate new Oxygen ml/l                            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    %create oxy coefficient variable
    oxy_coeff=zeros(7,1);
    
    %%%%%%%%%% calculating final statistics
    sub=length(new_oxy);
    d=boxy-new_oxy;
    m=nanmean(d);
    m_abs=nanmean(abs(d)./boxy)*100;
    s=nanstd(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(nanmean(abs(d)),4)]);
    disp(['fit mean absolute residual (%)  ',num2str((nanmean(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
    
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               End of Calibration Calculations                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

   
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                   Begin Plotting Calibrated data                        %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    %Convert back to micromoles
    new_oxmm = ox_units(new_oxy,sal,sw_ptmp(sal,temp,prs,0),'mll2mmk');
    
    

    %Oxygen vs. Pressure
    fcnt=0;
    fcnt=fcnt+1;
    fh(fcnt)=figure;
    subplot(2,1,1)
    %%%% Plot Oxygen Difference vs pressure (full profile)
    plot(prs,boxy-new_oxmm,'go');
    
       diff_o2 = boxy-new_oxmm;
       diff_o2_pct=diff_o2./boxy;
       mumol=nanmean(diff_o2);
       sumol=nanstd(diff_o2);
       grid on;
       ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
       title(['Calibrated',sensor,'Sensor',' Mean Diff = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
    
     
     subplot(2,1,2)
     %Find variable values below 500 db
     deep=find(prs>500);
     %%%% Plot Oxygen Difference vs pressure (<500 db)
     plot(prs(deep),boxy(deep)-new_oxmm(deep),'bo');        
        
        diff_o2 = boxy(deep)-new_oxmm(deep);
        diff_o2_pct=diff_o2./boxy(deep);
        mumol=nanmean(diff_o2);
        sumol=nanstd(diff_o2);
        grid on;
        ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
        xlabel('Pressure');
        title(['Calibrated',sensor,'Sensor',' Mean Diff (Below 500 db) = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
      
    %Oxygen vs. Station  
    fcnt=fcnt+1;
    fh(fcnt)=figure;
    subplot(2,1,1)
    %%%% Plot Oxygen Difference vs pressure (full profile)
    plot(stn,boxy-new_oxmm,'go');
    
       diff_o2 = boxy-new_oxmm;
       diff_o2_pct=diff_o2./boxy;
       mumol=nanmean(diff_o2);
       sumol=nanstd(diff_o2);
       grid on;
       ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
       %xlabel('Station');
       title(['Calibrated',sensor,'Sensor',' Mean Diff = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
    
     
     subplot(2,1,2)
     %Find variable values below 500 db
     deep=find(prs>500);
     %%%% Plot Oxygen Difference vs pressure (<500 db)
     plot(oxy_match.stnnbr(deep),boxy(deep)-new_oxmm(deep),'bo');        
        
        diff_o2 = boxy(deep)-new_oxmm(deep);
        diff_o2_pct=diff_o2./boxy(deep);
        mumol=nanmean(diff_o2);
        sumol=nanstd(diff_o2);
        grid on;
        ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
        xlabel('Station');
        title(['Calibrated',sensor,'Sensor',' Mean Diff (Below 500 db) = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
     
return;
