% calibrate_sal_oxygen_secondary
%
% WOCE standard for accuracy is 2e-3 psu (salinity) and 
%   1% of saturation (oxygen), ~0.05 ml/l (larger for T<20).
%
% This version of the code comes up with calibration coefficients for
% secondary sensors. (see calibrate_sal_oxygen_secmary.m for secmary
% sensor analysis).

close all;
clear;

% This version used secondary density to convert bottle oxygen data to mL/L
%load /Users/renellys.perez/Breve/Observations/PNE2017_RB1701/calibration/results/calibration_data_dens2.mat

% This version used primary density to convert bottle oxygen data to mL/L
load /Users/renellys.perez/Breve/Observations/PNE2017_RB1701/calibration/results/calibration_data.mat

% create plots for secondary salinity vs. autosal values 
figure
subplot(2,2,1)
plot(samples(:,4),samples(:,10),'bo','MarkerSize',6,'MarkerFaceColor','b');
hold on;
plot([min(samples(:,4)) max(samples(:,4))], ...
        [min(samples(:,4)) max(samples(:,4))],'k','LineWidth',2);
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('sec sal','Fontsize',12,'Fontname','Palatino'); ylabel('autosal','Fontsize',12,'Fontname','Palatino'); 
title('Salinity','Fontsize',14,'Fontname','Palatino');
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')

% create plots for secondary oxygen vs. oxygen titratation values 
subplot(2,2,3)
plot(samples(:,6),samples(:,11),'bo','MarkerSize',6,'MarkerFaceColor','b');
hold on;
plot([min(samples(:,6)) max(samples(:,6))], ...
        [min(samples(:,6)) max(samples(:,6))],'k','LineWidth',2);
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('sec oxy','Fontsize',12,'Fontname','Palatino'); ylabel('titrated','Fontsize',12,'Fontname','Palatino'); 
title('Oxygen','Fontsize',14,'Fontname','Palatino');
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')

% Drop the obvious outliers as bad autosal runs
% salinity: del = secondary sensor minus autosal
del=samples(:,4)-samples(:,10);
figure
hist(del,50); xlabel('salinity differences, sec-autosal','Fontsize',12,'Fontname','Palatino')
title('dropping outliers','Fontsize',14,'Fontname','Palatino')
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
pause(.5);
q=find(abs(del-nanmean(del))>3*nanstd(del));
while length(q)
    samples(q,10)=NaN;
    del=samples(:,4)-samples(:,10);
    hist(del,50); xlabel('salinity differences, sec-autosal','Fontsize',12,'Fontname','Palatino')
    title('dropping outliers','Fontsize',14,'Fontname','Palatino')
    set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
    q=find(abs(del-nanmean(del))>3*nanstd(del));
    pause(.5);
end;
print -depsc secondary_hist_salinity.eps

% Drop the obvious outliers as bad titrated runs
% oxygen: del=secondary sensor minus titrated value
del=samples(:,6)-samples(:,11);
figure
hist(del,50); xlabel('oxygen differences, sec-titration','Fontsize',12,'Fontname','Palatino')
title('dropping outliers','Fontsize',14,'Fontname','Palatino')
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
pause(.5);
q=find(abs(del-nanmean(del))>3*nanstd(del));
while length(q)
    samples(q,11)=NaN;
    del=samples(:,6)-samples(:,11);
    hist(del,50); xlabel('oxygen differences, sec-titration','Fontsize',12,'Fontname','Palatino')
    title('dropping outliers','Fontsize',14,'Fontname','Palatino')
    set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
    q=find(abs(del-nanmean(del))>3*nanstd(del));
    pause(.5);
end;
print -depsc secondary_hist_oxygen.eps

% create plots for secondary salinity vs. autosal values, outliers removed
figure(1)
subplot(2,2,2)
plot(samples(:,4),samples(:,10),'bo','MarkerSize',6,'MarkerFaceColor','b');
hold on;
plot([min(samples(:,4)) max(samples(:,4))], ...
        [min(samples(:,4)) max(samples(:,4))],'k','LineWidth',2);
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('sec sal','Fontsize',12,'Fontname','Palatino'); ylabel('autosal','Fontsize',12,'Fontname','Palatino'); 
title('Salinity (dropping outliers)','Fontsize',14,'Fontname','Palatino');
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')

% create plots for secondary oxygen vs. oxygen titratation values, outliers removed
subplot(2,2,4)
plot(samples(:,6),samples(:,11),'bo','MarkerSize',6,'MarkerFaceColor','b');
hold on;
plot([min(samples(:,6)) max(samples(:,6))], ...
        [min(samples(:,6)) max(samples(:,6))],'k','LineWidth',2);
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('sec oxy','Fontsize',12,'Fontname','Palatino'); ylabel('titrated','Fontsize',12,'Fontname','Palatino'); 
title('Oxygen (dropping outliers)','Fontsize',14,'Fontname','Palatino');
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
print -depsc secondary_scatterplot.eps


% get regression for secondary salinity residuals, using model for del: 
% del=mean(del) + x(1) (constant offset)
%   + x(2) * temp  + x(3) * temp^2 
%   + x(4) * press + x(5) * press^2
%   + x(6) * castnum + x(8)*castnum^2
%   + x(8) * sensor value + x(9) * sensor^2
%  del = A x + n

% salinity: del=secondary sensor minus autosal value (true value)
%   so true value = sensor - del
del=samples(:,4)-samples(:,10);
goodind=find(isfinite(del));
del=del(goodind);
As=NaN*ones(length(del),1);
As(:,1)=1;    % mean offset
As(:,2)=samples(goodind,8);  % secondary temp
As(:,3)=As(:,2).^2;
As(:,4)=samples(goodind,9);  % pressure
As(:,5)=As(:,4).^2;
As(:,6)=samples(goodind,1);  % cast num
As(:,7)=As(:,6).^2;
As(:,8)=samples(goodind,4);  % sec sal
As(:,9)=As(:,8).^2;

xs=pinv(As)*del;
figure
subplot(211)
plot(1:length(del),As*xs,'LineWidth',1);
hold on;
plot(1:length(del),del,'.','MarkerSize',8);
plot([1 length(del)],[0 0],'k--','LineWidth',1);
set(gca,'Xlim',[1 length(del)]);
title('(secondary S sensor - autosal value) and fit','Fontsize',14,'Fontname','Palatino')
ylabel('psu','Fontsize',12,'Fontname','Palatino')
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
subplot(223);
resid=As*xs-del;
plot(1:length(del),resid,'b.','MarkerSize',8);
axis([1 length(del) -max(abs(resid)) +max(abs(resid))])
title('residuals','Fontsize',14,'Fontname','Palatino')
ylabel('psu','Fontsize',12,'Fontname','Palatino');
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
subplot(224);
hist(resid,25);
ax=max(abs(resid)); axis([-ax ax 0 150])
title('histogram of resids','Fontsize',14,'Fontname','Palatino')
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
fprintf(['Error for calibrated salinity: ',num2str(1e3*std(resid),'%3.2f'),'e-3 psu.\n']);
disp([' (WOCE standard: 2e-3 psu)'])
resid_sal=NaN*ones(size(samples,1),1);
resid_sal(goodind)=resid;
print -depsc secondary_modelfit_salinity.eps


% get regression for secondary oxygen residuals, using model for del: 
% model for del: del=mean(del) + x(1) (constant offset)
%   + x(2) * temp  + x(3) * temp^2 
%   + x(4) * press + x(5) * press^2
%   + x(6) * castnum + x(8)*castnum^2
%   + x(8) * sensor value + x(9) * sensor^2
%  del = A x + n

% oxygen: del=secondary sensor minus titrated value (true value)
%   so true value = sensor - del
del=samples(:,6)-samples(:,11);
goodind=find(isfinite(del));
del=del(goodind);
Ao=NaN*ones(length(del),1);
Ao(:,1)=1;    % mean offset
Ao(:,2)=samples(goodind,8);  % secondary temp
Ao(:,3)=Ao(:,2).^2;
Ao(:,4)=samples(goodind,9);  % pressure
Ao(:,5)=Ao(:,4).^2;
Ao(:,6)=samples(goodind,1);  % cast num
Ao(:,7)=Ao(:,6).^2;
Ao(:,8)=samples(goodind,6);  % sec oxygen
Ao(:,9)=Ao(:,8).^2; 

xo=pinv(Ao)*del;
figure
subplot(211)
plot(1:length(del),Ao*xo,'LineWidth',1);
hold on;
plot(1:length(del),del,'.','MarkerSize',8);
plot([1 length(del)],[0 0],'k--','LineWidth',1);
set(gca,'Xlim',[1 length(del)]);
title('(secondary O2 sensor - titrated value) and fit','Fontsize',14,'Fontname','Palatino')
ylabel('ml / l','Fontsize',12,'Fontname','Palatino')
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
subplot(223);
resid=Ao*xo-del;
plot(1:length(del),resid,'b.','MarkerSize',8);
axis([1 length(del) -max(abs(resid)) +max(abs(resid))])
title('residuals','Fontsize',14,'Fontname','Palatino')
ylabel('ml / l','Fontsize',12,'Fontname','Palatino');
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
subplot(224);
hist(resid,25);
ax=max(abs(resid)); axis([-ax ax 0 150])
title('histogram of resids','Fontsize',14,'Fontname','Palatino')
set(gca,'Fontsize',12,'Fontname','Palatino','TickDir','out')
fprintf(['Error for calibrated oxygen: ',num2str(100*std(resid),'%3.2f'),'e-2 ml/l.\n']);
disp([' (WOCE standard: ???)'])
resid_oxy=NaN*ones(size(samples,1),1);
resid_oxy(goodind)=resid;
print -depsc secondary_modelfit_oxygen.eps

save calibration_coeff_secondary xs xo resid_*

%use convert_data.m to produce final CTD profiles