clear all

figure(1); clf; paper

%% First, mlr of t s o2 on pH
% running will save file mlr_pred
% load predict;	% t s o2 ph and predicted TCO2 TA 
% X = [t s o2]; [beta stats] = mlr(X,ph,1); 
load mlr_pred
rmse = sqrt(mean((y-pred).^2));
disp(['Root Mean Square Error:  ' num2str(rmse)]);
axes('position',[.25 .3 .5 .4]);
plot(y,pred,'k.');
mn = min([y;pred]); mx = max([y;pred]);
axis([mn mx mn mx]); hold on; plot([mn mx],[mn mx],'k');
set(gca,'fontweight','bold','fontsize',12);
xlabel('measured pH'); ylabel('predicted pH');
% title('pH prediction from measured T, S, O_2 (r^2 = 0.99)');
t=text(7.6,7.9,'r^2 = 0.99'); set(t,'fontsize',12);
t=text(7.6,7.85,'rmse = 0.02'); set(t,'fontsize',12);

return

%% Next, mlr of TCO2 & TA used to compute predicted pH
load predict;	% t s o2 ph and predicted TCO2 TA 
load predict_out.asc	% results of computing pH from predicted TCO2/TA
correct=0;

%% Columns 7:9 are pHT, pHF, and pHSWS
scale = 9;
% pHpred = [predict_out_1(:,scale); predict_out_2(:,scale)];
pHpred = [predict_out(:,scale)];

%% Convert measured pH to SWS 
%% ("measured" pH is based on NBS buffer calibration, and NBS buffers
%%  are on the molality scale (mol/Kg water) whereas you want the
%%  seawater scale (mol/Kg solution).  This will reduce "measured"
%%  pH).  M.E.Q. Pilson provides equations for conversion
%% At the salinities you are working with, the average conversion factor 
%% is .9661 and the stdev is .0003; 
% phsws = (1 - .00100575*s).*ph;
%% Wrong: you are not looking at concentration scales, you just have
%% to account for the effects of density
%% NBS is mol/kg of water ~ mol/l
%% SWS is mol/kg solution  
%% So, divide predicted by density in kg/l
if correct
 display(['Adjusting pH scale']);
 addpath '/u/ryjo/matlab/CSIRO'
 lat=repmat(36.8,length(ph),1);
 pres = sw_pres(Depth,lat);      % [db]
 dens = sw_dens(s,t,pres);	% kg/m3
 dens = dens/1000;		% kg/l
 phsws = ph./dens;
 figure(1); clf; bigfig;
 plot(ph,phsws); xlabel('pH'); ylabel('pH_{SWS}');
 ph=phsws;
end

idx = find(Depth >0);

subplot(212);
plot(ph(idx),pHpred(idx),'ko');
mn = min([ph(idx);pHpred(idx)]); mx = max([ph(idx);pHpred(idx)]);
axis([mn mx mn mx]); hold on; plot([mn mx],[mn mx],'k');
set(gca,'fontweight','bold','fontsize',12);
xlabel('measured pH'); ylabel('predicted pH');
title('pH prediction via C_T and A_T predicted from measured T, S, O_2 ');

rmse = sqrt(mean((ph(idx)-pHpred(idx)).^2));
disp(['Root Mean Square Error:  ' num2str(rmse)]);

figure(2); clf; paper
subplot(321); plot(s,TA,'*'); xlabel('salinity'); ylabel('TA');
 axis([min(s) max(s) min(TA) max(TA)]);
subplot(322); plot(s,TCO2,'*'); xlabel('salinity'); ylabel('TCO_2');
 axis([min(s) max(s) min(TCO2) max(TCO2)]);
subplot(323); plot(t,TA,'*'); xlabel('temperature'); ylabel('TA');
 axis([min(t) max(t) min(TA) max(TA)]);
subplot(324); plot(t,TCO2,'*'); xlabel('temperature'); ylabel('TCO_2');
 axis([min(t) max(t) min(TCO2) max(TCO2)]);
subplot(325); plot(o2,TA,'*'); xlabel('oxygen'); ylabel('TA');
 axis([min(o2) max(o2) min(TA) max(TA)]);
subplot(326); plot(o2,TCO2,'*'); xlabel('oxygen'); ylabel('TCO_2');
 axis([min(o2) max(o2) min(TCO2) max(TCO2)]);

frame
t=text(.5,.95,'Predictions from AUV T, S, O2');
set(t,'fontsize',12,'horizontalalignment','center');

