function [fh]=eval_calibration_jay(cruiseid,oxy_coeff)

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

disp('**************************************************************************************')
disp(' 1 = use the primary sensor (DEFAULT)' )
disp(' 2 = use the secondary sensor' )
disp('For the oxygen sensor:')
oxy_sensor=input('Do you want apply the calibration to the primary or to the secondary sensor: ');
while (isempty(oxy_sensor)==1 | oxy_sensor<=0 | oxy_sensor>2)
	disp('Please answer 0 (zero) for no or  1(one) for yes ')
	oxy_sensor=input('Do you want to use the primary or the secondary sensor: ');
end
disp('***************************************************************************')
staa=input('Input first station: ');
stab=input('Input last station: ');

%%%%%%%%%%%%%%%%%%%%%% Good data statistics %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%for stt=staa:stab,
    %stn=stt;
    %Define quality control filter
    %ii=find(oxy_match.quality1 == 2 & oxy_match.stnnbr == stt);
    ii=find(oxy_match.quality1 == 2 & oxy_match.stnnbr >= staa & ...
            oxy_match.stnnbr <= stab);
        
    if oxy_sensor == 1,
        %creat variables for primary oxy_match
        %disp('Primary Sensor')
        stn=oxy_match.stnnbr(ii);
        oxy_volt=oxy_match.ctdoxv1(ii);
        oxy_mll=oxy_match.ctdoxy1(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);
    else
        %creat variables for secondary oxy_match
        %disp('Secondary Sensor')
        stn=oxy_match.stnnbr(ii);
        oxy_volt=oxy_match.ctdoxv2(ii);
        oxy_mll=oxy_match.ctdoxy2(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);
    end;
    
    %calculate the saturated oxygen
    %sal_oxy=sw_salt(cond/sw_c3515,temp,prs);
    oxsat=sw_satO2(temp,sal);
    %boxyL=ox_units(boxy,sal,sw_ptmp(sal,temp,prs,0),'mmk2mll');
  
    % Define oxygen coefficients
    Soc=oxy_coeff(1);
    Voffset=oxy_coeff(2);
    A=oxy_coeff(3);
    B=oxy_coeff(4);
    C=oxy_coeff(5);
    E=oxy_coeff(6);
    p1=oxy_coeff(7);
    %stn=repmat(stn,length(oxy_volt),1);
 
    
    
    
    %%%% Calculate new oxygen using new oxygen coefficents
    new_oxy=((Soc.*(oxy_volt+Voffset))+(p1.*stn)).*...
        (1.0+(A.*temp)+(B.*temp.^2)+(C.*temp.^3)).*oxsat.*...
        exp((E.*prs)./(temp+273.15));
    %%%% Convert new_oxy ml/l to umol/kg    
    new_oxmm = ox_units(new_oxy,sal,sw_ptmp(sal,temp,prs,0),'mll2mmk');
    
    
   %return 
    %create structure variables
    oxy_match_cal.stnnbr(ii,1)=oxy_match.stnnbr(ii);
    oxy_match_cal.niskinbtlpos(ii,1)=oxy_match.niskinbtlpos(ii);
    oxy_match_cal.sampnbr(ii,1)=oxy_match.sampnbr(ii);
    oxy_match_cal.btloxy(ii,1)=oxy_match.btloxy(ii);
    oxy_match_cal.quality1(ii,1)=oxy_match.quality1(ii);
    oxy_match_cal.quality2(ii,1)=oxy_match.quality2(ii);
    oxy_match_cal.ctdprs(ii,1)=prs;
    oxy_match_cal.ctdtmp(ii,1)=temp;
    oxy_match_cal.ctdcon(ii,1)=cond;
    oxy_match_cal.ctdsal(ii,1)=sal; 
    oxy_match_cal.ctdoxy(ii,1)=new_oxy; 
    oxy_match_cal.ctdoxyumolkg(ii,1)=new_oxmm; 
   
%end;% end for station for loop 
%return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%               End of Calibration Calculations                           %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

    %Check for structure existance and save if not
    
        %disp('Saving oxy_match_cal structure')
        %save(db_file,'oxy_match_cal','-append');
    

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                           Begin Plotting                                %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %staa=0;
    %stab=58;
    plot_index=find(oxy_match_cal.quality1 == 2 & oxy_match_cal.stnnbr >= staa ...
                    & oxy_match_cal.stnnbr <= stab);
    %Figure 1
    fcnt=0;
    fcnt=fcnt+1;
    fh(fcnt)=figure;
    subplot(2,1,1)
    %%%% Plot Oxygen Difference vs pressure (full profile)
    plot(oxy_match_cal.ctdprs(plot_index),...
        oxy_match_cal.btloxy(plot_index)-oxy_match_cal.ctdoxyumolkg(plot_index),'go');
    
       diff_o2 = oxy_match_cal.btloxy(plot_index)-oxy_match_cal.ctdoxyumolkg(plot_index);
       diff_o2_pct=diff_o2./oxy_match_cal.btloxy(plot_index);
       mumol=nanmean(diff_o2);
       sumol=nanstd(diff_o2);
       grid on;
       ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
       xlabel('Pressure');
       title([' Mean Diff = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
    
     %Figure 2   
     %fcnt=fcnt+1;
     %fh(fcnt)=figure;
     subplot(2,1,2)
     %Find variable values below 500 db
     deep=find(oxy_match_cal.ctdprs>500 & oxy_match.quality1 == 2 ...
                & oxy_match.stnnbr >= staa & oxy_match.stnnbr <= stab);
     %%%% Plot Oxygen Difference vs pressure (<500 db)
     plot(oxy_match_cal.ctdprs(deep),...
         oxy_match_cal.btloxy(deep)-oxy_match_cal.ctdoxyumolkg(deep),'bo');        
        
        diff_o2 = oxy_match_cal.btloxy(deep)-oxy_match_cal.ctdoxyumolkg(deep);
        diff_o2_pct=diff_o2./oxy_match_cal.btloxy(deep);
        mumol=nanmean(diff_o2);
        sumol=nanstd(diff_o2);
        grid on;
        ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
        xlabel('Pressure');
        title([' Mean Diff (Below 500 db) = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
      
    %Figure 5    
    fcnt=fcnt+1;
    fh(fcnt)=figure;
    subplot(2,1,1)
    %%%% Plot Oxygen Difference vs pressure (full profile)
    plot(oxy_match_cal.stnnbr(plot_index),...
        oxy_match_cal.btloxy(plot_index)-oxy_match_cal.ctdoxyumolkg(plot_index),'go');
    
       diff_o2 = oxy_match_cal.btloxy(plot_index)-oxy_match_cal.ctdoxyumolkg(plot_index);
       diff_o2_pct=diff_o2./oxy_match_cal.btloxy(plot_index);
       mumol=nanmean(diff_o2);
       sumol=nanstd(diff_o2);
       grid on;
       ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
       %xlabel('Station');
       title([' Mean Diff = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
    
     %Figure 4   
     %fcnt=fcnt+1;
     %fh(fcnt)=figure;
     subplot(2,1,2)
     %Find variable values below 500 db
     deep=find(oxy_match_cal.ctdprs>500 & oxy_match.quality1 == 2 ...
                & oxy_match.stnnbr >= staa & oxy_match.stnnbr <= stab);
     %%%% Plot Oxygen Difference vs pressure (<500 db)
     plot(oxy_match_cal.stnnbr(deep),...
         oxy_match_cal.btloxy(deep)-oxy_match_cal.ctdoxyumolkg(deep),'bo');        
        
        diff_o2 = oxy_match_cal.btloxy(deep)-oxy_match_cal.ctdoxyumolkg(deep);
        diff_o2_pct=diff_o2./oxy_match_cal.btloxy(deep);
        mumol=nanmean(diff_o2);
        sumol=nanstd(diff_o2);
        grid on;
        ylabel('Bottle -  CTD Oxygen (\mu mol/kg)');
        xlabel('Station');
        title([' Mean Diff (Below 500 db) = ',num2str(mumol,4),'   std.dev = ',num2str(sumol,4)]);
     
return;
