% calibrate_sal_oxygen
%
% WOCE standard for accuracy is 2e-3 psu (salinity) and 
%   1% of saturation (oxygen), ~0.05 ml/l (larger for T<20).
%
%  NOTE: the values for "btl().data" are different from cruise
%   to cruise.  CHECK CAREFULLY before applying this by examining
%   btl().names
   
%{
clear all;  close all;

if (1)  % load all data, create matrix "samples" containing
        % simultaneous sensor and bottle values
        % (set to zero to load already-saved "samples", e.g. if
        %  tweaking code to adjust sensor values to match bottle values).
        
        % samples(:,1)  -   cast number
        % samples(:,2)  -   niskin number
        % samples(:,3)  -   primary salinity (sensor value at bottle depth)
        %   (:,4)       -   secondary salinity (sensor)
        %   (:,5)       -   primary oxygen (sensor value at bottle depth)
        %   (:,6)       -   secondary oxygen  (both oxy in ml/l)
        %   (:,7)       -   primary temp
        %   (:,8)       -   secondary temp
        %   (:,9)       -   pressure
        %   (:,10)      -   salinity value from autosal (bottle value)
        %   (:,11)      -   oxygen value from titration (bottle value)
%}


%load('calibration_data.mat');
load('/Users/gregory.foltz/pne_cruises/pne17/proc_data/calibration_data.mat');

% create plots
figure
subplot(2,2,1)
plot(samples(:,3),samples(:,10),'.');
hold on;
plot([min(samples(:,3)) max(samples(:,3))], ...
        [min(samples(:,3)) max(samples(:,3))],'k');
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('pri sal'); ylabel('autosal');

subplot(2,2,2);
plot(samples(:,4),samples(:,10),'.');
hold on;
plot([min(samples(:,4)) max(samples(:,4))], ...
        [min(samples(:,4)) max(samples(:,4))],'k');
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('sec sal'); ylabel('autosal');

subplot(2,2,3);
plot(samples(:,5),samples(:,11),'.');
hold on;
plot([min(samples(:,5)) max(samples(:,5))], ...
        [min(samples(:,5)) max(samples(:,5))],'k');
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('pri oxy'); ylabel('titrated');

subplot(2,2,4);
plot(samples(:,6),samples(:,11),'.');
hold on;
plot([min(samples(:,6)) max(samples(:,6))], ...
        [min(samples(:,6)) max(samples(:,6))],'k');
hold off;
set(gca,'dataaspectratio',[1 1 1]);
xlabel('sec oxy'); ylabel('titrated');
 

% drop the obvious outliers as bad titration/autosal runs

if (1)
% salinity: del = primary sensor minus autosal
del=samples(:,3)-samples(:,10);
figure
hist(del,50); xlabel('salinity differences, pri-autosal')
title('dropping outliers')
pause(.5);
q=find(abs(del-nanmean(del))>3*nanstd(del));
while length(q)
    samples(q,10)=NaN;
    del=samples(:,3)-samples(:,10);
    hist(del,50); xlabel('salinity differences, pri-autosal')
    title('dropping outliers')
    q=find(abs(del-nanmean(del))>3*nanstd(del));
    pause(.5);
end;

end;

% oxygen: del=primary sensor minus titrated value
del=samples(:,5)-samples(:,11);

% PNE11: bad values for titrated flasks beneath the ML in casts
% 1-4
%q=find(samples(:,11)>4.4 & ...
%    samples(:,9)>150 & samples(:,9)<800);    % pressure
%samples(q,11)=NaN;

figure
hist(del,50); xlabel('oxygen differences, pri-titration')
title('dropping outliers')
pause(.5);
q=find(abs(del-nanmean(del))>3*nanstd(del));
while length(q)
    samples(q,11)=NaN;
    del=samples(:,5)-samples(:,11);
    hist(del,50); xlabel('oxygen differences, pri-titration')
    title('dropping outliers')
    q=find(abs(del-nanmean(del))>3*nanstd(del));
    pause(.5);
end;

%save calibrations samples

%else
%    load calibrations;
%end;


% salinity: del=primary sensor minus autosal value (true value)
%   so true value = sensor - del
del=samples(:,3)-samples(:,10);
goodind=find(isfinite(del));
del=del(goodind);
% 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(7)*castnum^2
%   + x(8) * sensor value + x(9) * sensor^2
%  del = A x + n
As=NaN*ones(length(del),1);
As(:,1)=1;    % mean offset
As(:,2)=samples(goodind,7);  % primary temp
As(:,3)=As(:,2).^2;
As(:,4)=samples(goodind,9);  % pressure
As(:,5)=As(:,4).^2;
%keyboard
As(:,6)=samples(goodind,1);  % cast num
As(:,7)=As(:,6).^2;
As(:,8)=samples(goodind,3);  % pri sal
As(:,9)=As(:,8).^2;

xs=pinv(As)*del;
figure
subplot(211)
plot(1:length(del),As*xs,1:length(del),del,'.')
title('(primary S sensor - autosal value) and fit')
ylabel('psu')
subplot(223);
resid=As*xs-del;
plot(1:length(del),resid,'.');
title('residuals')
ylabel('psu');
subplot(224);
hist(resid,25);
title('histogram of resids')
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;


% oxygen: del=primary sensor minus titrated value (true value)
%   so true value = sensor - del
del=samples(:,5)-samples(:,11);
goodind=find(isfinite(del));
del=del(goodind);
% 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(7)*castnum^2
%   + x(8) * sensor value + x(9) * sensor^2
%  del = A x + n
Ao=NaN*ones(length(del),1);

Ao(:,1)=1;    % mean offset
Ao(:,2)=samples(goodind,7);  % primary 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,5);  % pri oxygen
Ao(:,9)=Ao(:,8).^2; 

xo=pinv(Ao)*del;
figure
subplot(211)
plot(1:length(del),Ao*xo,1:length(del),del,'.')
title('(primary O2 sensor - titrated value) and fit')
ylabel('ml / l')
subplot(223);
resid=Ao*xo-del;
plot(1:length(del),resid,'.');
title('residuals')
ylabel('ml / l');
subplot(224);
hist(resid,25);
title('histogram of resids')
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;

save('/Users/gregory.foltz/pne_cruises/pne17/proc_data/calibration_coeff.mat','xs','xo','resid_oxy','resid_sal');
%save calibration_coeff xs xo resid_*

%use convert_data_pne17.m to produce final CTD profiles