Fix inexact slope solver pH from TA

comparisons/compare_against_master.m

@@ -0,0 +1,98 @@
+% Compares your latest commit with CO2SYS v3.2.0
+
+
+
+%% Run CO2SYS v3.2.0
+disp('Running CO2SYS v3.2.0')
+rehash path % Circumvent caching of CO2SYS function
+system('git checkout v3.2.0'); % Checkout the v3.2.0 tag
+
+%% Set up input conditions
+PARvalues = [2250 2100 8.1 400 405];
+PARTYPEs = 1:5;
+pHSCALEIN_opts = 1:4;
+K1K2CONSTANTS_opts = 1:15;
+KSO4CONSTANTS_opts = 1:4;
+KFCONSTANT_opts = 1;
+SALvalue = 33.1;
+[P1, P2, P1type, P2type, sal, pHscales, K1K2, ~, KSO4, KF, ...
+    BSal, ~, ~] = CO2SYSigen(PARvalues, PARTYPEs, SALvalue, pHSCALEIN_opts, ...
+    K1K2CONSTANTS_opts, KSO4CONSTANTS_opts, KFCONSTANT_opts);
+tempin = 24;
+tempout = 12;
+presin = 1;
+presout = 1647;
+si = 10;
+phos = 1;
+
+%% Determine whether to calculate each input row or not
+% xrow = 1 + 210; % just do one row, or...
+xrow = 1:numel(P1); % ... do all rows (do this for saving output file)
+P1 = P1(xrow);
+P2 = P2(xrow);
+P1type = P1type(xrow);
+P2type = P2type(xrow);
+sal = sal(xrow);
+pHscales = pHscales(xrow);
+K1K2 = K1K2(xrow);
+
+tic
+[DATA_v3, HEADERS_v3] = ...
+    CO2SYS(P1, P2, P1type, P2type, sal, tempin, tempout, presin, ...
+    presout, si, phos, 0, 0, pHscales, K1K2, KSO4, KF, BSal);
+toc
+
+
+
+%% Run current commit CO2SYS
+disp('Running CO2SYS current commit/head')
+rehash path % Circumvent caching of CO2SYS function
+system('git switch -'); % Switch back to your commit
+
+%% Set up input conditions
+PARvalues = [2250 2100 8.1 400 405];
+PARTYPEs = 1:5;
+pHSCALEIN_opts = 1:4;
+K1K2CONSTANTS_opts = 1:15;
+KSO4CONSTANTS_opts = 1:4;
+KFCONSTANT_opts = 1;
+SALvalue = 33.1;
+[P1, P2, P1type, P2type, sal, pHscales, K1K2, ~, KSO4, KF, ...
+    BSal, ~, ~] = CO2SYSigen(PARvalues, PARTYPEs, SALvalue, pHSCALEIN_opts, ...
+    K1K2CONSTANTS_opts, KSO4CONSTANTS_opts, KFCONSTANT_opts);
+tempin = 24;
+tempout = 12;
+presin = 1;
+presout = 1647;
+si = 10;
+phos = 1;
+
+%% Determine whether to calculate each input row or not
+% xrow = 1 + 210; % just do one row, or...
+xrow = 1:numel(P1); % ... do all rows (do this for saving output file)
+P1 = P1(xrow);
+P2 = P2(xrow);
+P1type = P1type(xrow);
+P2type = P2type(xrow);
+sal = sal(xrow);
+pHscales = pHscales(xrow);
+K1K2 = K1K2(xrow);
+
+tic
+[DATA, HEADERS] = ...
+    CO2SYS(P1, P2, P1type, P2type, sal, tempin, tempout, presin, ...
+    presout, si, phos, 0, 0, pHscales, K1K2, KSO4, KF, BSal);
+toc
+
+fprintf("\n\nRelative change vs v3.2.0:\n"); ...
+fprintf("%20s %20s %20s %20s %20s\n", "Variable", "Mean rel. change", "Min rel. change", "Max rel. change", "# of samples"); ...
+for V = 1:length(HEADERS_v3)
+    x = DATA_v3(:,V);
+    y = DATA(:,V);
+    relerr = abs(y - x) ./ abs(x);
+    ix = x ~= -999; % Only compare non fill-in values
+    maxrelerr = max(relerr(idx));
+    minrelerr = min(relerr(idx));
+    meanrelerr = mean(relerr(idx));
+    fprintf("%20s %20.2g %20.2g %20.2g %20i\n", HEADERS_v3{V}, meanrelerr, minrelerr, maxrelerr, length(ix))
+end

main/CO2SYS.m

@@ -339,7 +339,6 @@
 % NOTHING BELOW THIS SHOULD REQUIRE EDITING BY USER!
 %**************************************************************************
 
-
 % Declare global variables
 global pHScale WhichKs WhoseKSO4 WhoseKF WhoseTB Pbar
 global Sal sqrSal TempK logTempK TempCi TempCo Pdbari Pdbaro;
@@ -1903,34 +1902,56 @@ function Constants(TempC,Pdbar)
 nF=(abs(deltapH) > pHTol);
 while any(nF)
     H         = 10.^(-pH);
+    dH_dpHx = -ln10 .* H;
     Denom     = (H.*H + K1F.*H + K1F.*K2F);
+    dDenom_dH = K1F + 2 * H;
     CAlk      = TCi.*K1F.*(H + 2.*K2F)./Denom;
+    dCAlk_dH = (K1F .* TCi - dDenom_dH .* CAlk) ./ Denom;
     BAlk      = TBF.*KBF./(KBF + H);
+    dBAlk_dH = -BAlk ./ (KBF + H);
     OH        = KWF./H;
+    dOH_dH = -KWF ./ (H .* H);
     PhosTop   = KP1F.*KP2F.*H + 2.*KP1F.*KP2F.*KP3F - H.*H.*H;
+    dPhosTop_dH = KP1F .* KP2F - 3 * H .* H;
     PhosBot   = H.*H.*H + KP1F.*H.*H + KP1F.*KP2F.*H + KP1F.*KP2F.*KP3F;
+    dPhosBot_dH = 3 * H .* H + KP1F .* KP2F + 2 * H .* KP1F;
     PAlk      = TPF.*PhosTop./PhosBot;
+    dPAlkex_dPhosTop = TPF ./ PhosBot;
+    dPAlkex_dPhosBot = -PhosTop .* TPF ./ (PhosBot .* PhosBot);
+    dPAlk_dH = dPAlkex_dPhosTop .* dPhosTop_dH + dPAlkex_dPhosBot .* dPhosBot_dH;
     SiAlk     = TSiF.*KSiF./(KSiF + H);
+    dSiAlk_dH = -SiAlk ./ (H + KSiF);
     AmmAlk    = TNH4F.*KNH4F./(KNH4F + H);
+    dAmmAlk_dH = -AmmAlk ./ (H + KNH4F);
     HSAlk     = TH2SF.*KH2SF./(KH2SF + H);
-    [~,~,pHfree,~] = FindpHOnAllScales(pH); % this converts pH to pHfree no matter the scale
-    Hfree     = 10.^-pHfree; % this converts pHfree to Hfree
+    dHSAlk_dH = -HSAlk ./ (H + KH2SF);
+    % Changing from pH to pHfree is just a tranlsation (± C)
+    % and changing from H to Hfree is just a scaling by C
+    % so all we need to do is grab that constant to get dH_scaleA/dH_scaleB = 1/C
+    [~,~,pHfree,~,~,~,scaling_constant] = FindpHOnAllScales(pH); % this converts pH to pHfree no matter the scale
+    Hfree     = 10.^(-pHfree); % this converts pHfree to Hfree
+    dHfree_dH = 1 ./ scaling_constant; %
     HSO4      = TSF./(1 + KSF./Hfree); % since KS is on the free scale
+    dHSO4_dHfree = KSF ./ (Hfree .* Hfree) .* HSO4 ./ (1 + KSF ./ Hfree);
+    dHSO4_dH = dHSO4_dHfree .* dHfree_dH;
     HF        = TFF./(1 + KFF./Hfree); % since KF is on the free scale
-    Residual  = TAi - CAlk - BAlk - OH - PAlk - SiAlk  - AmmAlk - HSAlk + Hfree + HSO4 + HF;
+    dHF_dHfree = KFF ./ (Hfree .* Hfree) .* HF ./ (1 + KFF ./ Hfree);
+    dHF_dH = dHF_dHfree .* dHfree_dH;
+    Residual  = TAi - CAlk - BAlk - OH - PAlk - SiAlk - AmmAlk - HSAlk + Hfree + HSO4 + HF;
     % find Slope dTA/dpH;
-    % (this is not exact, but keeps all important terms);
-    Slope     = ln10.*(TCi.*K1F.*H.*(H.*H + K1F.*K2F + 4.*H.*K2F)./Denom./Denom + BAlk.*H./(KBF + H) + OH + H);
-    deltapH   = Residual./Slope; %' this is Newton's method
-    % ' to keep the jump from being too big:
+    % (this is now exact!)
+    Slope     = dH_dpHx .* (-dCAlk_dH - dBAlk_dH - dOH_dH - dPAlk_dH ...
+        - dSiAlk_dH - dAmmAlk_dH - dHSAlk_dH + dHfree_dH + dHSO4_dH + dHF_dH);
+    deltapH   = -Residual./Slope; % this is Newton's method
+    % to keep the jump from being too big:
     while any(abs(deltapH) > 1)
         FF=abs(deltapH)>1; deltapH(FF)=deltapH(FF)./2;
     end
     pH(nF) = pH(nF) + deltapH(nF);
     nF     = abs(deltapH) > pHTol;
     loopc=loopc+1;
- 
-    if loopc>10000
+
+    if loopc > 100 % BP: changed max iterations to 100
         Fr=find(abs(deltapH) > pHTol);
         pH(Fr)=NaN;  disp(['pH value did not converge for data on row(s): ' num2str((Fr)')]);
         deltapH=pHTol*0.9;
@@ -2762,7 +2783,7 @@ function Constants(TempC,Pdbar)
 end % end nested function
 
 function varargout=FindpHOnAllScales(pH)
-global pHScale K T TS KS TF KF fH F ntps;
+global pHScale TS KS TF KF fH F
 % ' SUB FindpHOnAllScales, version 01.02, 01-08-97, written by Ernie Lewis.
 % ' Inputs: pHScale%, pH, K(), T(), fH
 % ' Outputs: pHNBS, pHfree, pHTot, pHSWS
@@ -2776,17 +2797,25 @@ function Constants(TempC,Pdbar)
 nF=pHScale(F)==1;  %'"pHtot"
 factor(nF) = 0;
 nF=pHScale(F)==2; % '"pHsws"
-factor(nF) = -log(SWStoTOT(nF))./log(0.1);
+factor(nF) = log10(SWStoTOT(nF));
 nF=pHScale(F)==3; % '"pHfree"
-factor(nF) = -log(FREEtoTOT(nF))./log(0.1);
+factor(nF) = log10(FREEtoTOT(nF));
 nF=pHScale(F)==4;  %'"pHNBS"
-factor(nF) = -log(SWStoTOT(nF))./log(0.1) + log(fHx(nF))./log(0.1);
+factor(nF) = log10(SWStoTOT(nF)) - log10(fHx(nF));
 pHtot  = pH    - factor;    % ' pH comes into this sub on the given scale
-pHNBS  = pHtot - log(SWStoTOT) ./log(0.1) + log(fHx)./log(0.1);
-pHfree = pHtot - log(FREEtoTOT)./log(0.1);
-pHsws  = pHtot - log(SWStoTOT) ./log(0.1);
+pHNBS  = pHtot + log10(SWStoTOT ./ fHx);
+pHfree = pHtot + log10(FREEtoTOT);
+pHsws  = pHtot + log10(SWStoTOT);
 varargout{1}=pHtot;
 varargout{2}=pHsws;
 varargout{3}=pHfree;
 varargout{4}=pHNBS;
+% If requested, return the scaling constant
+% (useful for derivative of H(scale in) wrt H(scale out))
+if nargout > 4
+    varargout{5}=factor;
+    varargout{6}=factor - log10(SWStoTOT ./ fHx);
+    varargout{7}=factor - log10(FREEtoTOT);
+    varargout{8}=factor - log10(SWStoTOT);
+end
 end % end nested function