Continued scale conversion

Author: mvdh7

PyCO2SYS/__init__.py

@@ -19,7 +19,7 @@
     # api,
     # bio,
     # buffers,
-    # constants,
+    constants,
     convert,
     engine,
     equilibria,

PyCO2SYS/convert.py

@@ -1,7 +1,24 @@
 # PyCO2SYS: marine carbonate system calculations in Python.
 # Copyright (C) 2020--2024  Matthew P. Humphreys et al.  (GNU GPLv3)
-"""Convert units and calculate conversion factors."""
+"""
+PyCO2SYS.convert
+================
+Convert units and calculate conversion factors.
+
+pH scale conversions
+--------------------
+There is a function for every variant of pH_<scale1>_to_<scale2>, where <scale1> and
+<scale2> are any of free, total, sws and nbs.  Each function has a different set of
+arguments, depending on what is needed.  For example, to get the conversion factor to
+go from the total to the NBS scale, use:
+
+  >>> factor = pyco2.convert.pH_total_to_nbs(
+        total_fluoride, total_sulfate, k_HF_free, k_HSO4_free, fH
+      )
 
+``factor`` should be multiplied by [H+] or K value(s) on the total scale to convert to
+the NBS scale.
+"""
 import copy
 from jax import numpy as np
 from . import constants

PyCO2SYS/engine/system.py

@@ -1,10 +1,11 @@
 # PyCO2SYS: marine carbonate system calculations in Python.
 # Copyright (C) 2020--2024  Matthew P. Humphreys et al.  (GNU GPLv3)
+from operator import mul
 import itertools
 import networkx as nx
 from jax import numpy as np
 from matplotlib import pyplot as plt
-from . import convert, equilibria, salts
+from .. import constants, convert, equilibria, salts
 
 # Define functions for calculations that depend neither on icase nor opts:
 get_funcs = {
@@ -27,7 +28,11 @@
         k_H2S_total_1atm * total_to_sws_1atm
     ),
     # Pressure correction factors for equilibrium constants
+    "factor_k_H2S": equilibria.pcx.factor_k_H2S,
     # Equilibrium constants at pressure and on the seawater pH scale
+    "k_H2S_sws": lambda k_H2S_sws_1atm, factor_k_H2S: (
+        k_H2S_sws_1atm, factor_k_H2S
+    )
     # Equilibrium constants at pressure and on the requested pH scale (YES, HERE!)
 }
 
@@ -48,6 +53,11 @@
 # Define functions for calculations that depend on opts:
 # (unlike in previous versions, each opt may only affect one parameter)
 get_funcs_opts = {}
+get_funcs_opts["opt_gas_constant"] = {
+    1: dict(gas_constant=lambda: constants.RGasConstant_DOEv2),
+    2: dict(gas_constant=lambda: constants.RGasConstant_DOEv3),
+    3: dict(gas_constant=lambda: constants.RGasConstant_CODATA2018),
+}
 get_funcs_opts["opt_fH"] = {
     1: dict(fH=convert.fH_TWB82),
     2: dict(fH=convert.fH_PTBO87),
@@ -128,6 +138,7 @@
 }
 
 default_opts = {
+    "opt_gas_constant": 3,
     "opt_fH": 1,
     "opt_k_HF": 1,
     "opt_k_HSO4": 1,

PyCO2SYS/equilibria/p1atm.py

@@ -26,6 +26,7 @@
     with the corrections of WMW14.  Used when opt_k_HSO4 = 3.
 """
 from jax import numpy as np
+from .. import convert
 
 
 def kCO2_W74(temperature, salinity):
@@ -68,7 +69,7 @@ def k_HSO4_free_D90a(temperature, salinity, ionic_strength):
     # TYPO on p. 121: the constant e9 should be e8.
     # Output KS is on the free pH scale in mol/kg-sw.
     # This is from eqs 22 and 23 on p. 123, and Table 4 on p 121:
-    TempK = temperature + 273.15
+    TempK = convert.celsius_to_kelvin(temperature)
     logTempK = np.log(TempK)
     lnk_HSO4 = (
         -4276.1 / TempK
@@ -137,7 +138,7 @@ def k_HSO4_free_WM13(temperature, salinity):
     float
         HSO4 dissociation constant.
     """
-    TempK = temperature + 273.15
+    TempK = convert.celsius_to_kelvin(temperature)
     logKS0 = (
         562.69486
         - 102.5154 * np.log(TempK)
@@ -866,7 +867,7 @@ def k_H2S_total_YM95(temperature, salinity):
     # they agree with Millero 1988 which are on Total Scale.
     # Also, calculations agree at high H2S with AquaEnv when assuming it is on
     # Total Scale.
-    TempK = temperature + 273.15
+    TempK = convert.celsius_to_kelvin(temperature)
     lnkH2S = (
         225.838
         - 13275.3 / TempK

PyCO2SYS/equilibria/pcx.py

@@ -6,9 +6,30 @@
 from .. import convert
 
 
-def Kfac(deltaV, Kappa, Pbar, TempK, RGas):
-    """Calculate pressure correction factor for equilibrium constants."""
-    return np.exp((-deltaV + 0.5 * Kappa * Pbar) * Pbar / (RGas * TempK))
+def pressure_factor(deltaV, kappa, pressure, temperature, gas_constant):
+    """Calculate pressure-correction factors for equilibrium constants.
+    
+    Parameters
+    ----------
+    deltaV : float
+        The deltaV constant for the pressure correction.
+    kappa : float
+        The kappa constant for the pressure correction.
+    pressure : float
+        Hydrostatic pressure in dbar.
+    temperature : float
+        Temperature in °C.
+    gas_constant : float
+        The universal gas constant in ml / (bar * K * mol).
+    
+    Returns
+    -------
+    float
+        The correction factor, to be multiplied by the K value to correct it.
+    """
+    Pbar = convert.decibar_to_bar(pressure)
+    TempK = convert.celsius_to_kelvin(temperature)
+    return np.exp((-deltaV + 0.5 * kappa * Pbar) * Pbar / (gas_constant * TempK))
 
 
 def KSO4fac(TempK, Pbar, RGas):
@@ -134,16 +155,15 @@ def KSifac(TempK, Pbar, RGas):
     return Kfac(deltaV, Kappa, Pbar, TempK, RGas)
 
 
-def KH2Sfac(TempK, Pbar, RGas):
-    """Calculate pressure correction factor for KH2S."""
+def factor_k_H2S(temperature, pressure, gas_constant):
+    """Calculate pressure correction factor for k_H2S."""
     # === CO2SYS.m comments: =======
     # Millero 1995 gives values for deltaV in fresh water instead of SW.
     # Millero 1995 gives -b0 as -2.89 instead of 2.89.
     # Millero 1983 is correct for both.
-    TempC = convert.kelvin_to_celsius(TempK)
-    deltaV = -11.07 - 0.009 * TempC - 0.000942 * TempC**2
-    Kappa = (-2.89 + 0.054 * TempC) / 1000
-    return Kfac(deltaV, Kappa, Pbar, TempK, RGas)
+    deltaV = -11.07 - 0.009 * temperature - 0.000942 * temperature**2
+    kappa = (-2.89 + 0.054 * temperature) / 1000
+    return pressure_factor(deltaV, kappa, pressure, temperature, gas_constant)
 
 
 def KNH3fac(TempK, Pbar, RGas):

newclass_v2.py

@@ -1,3 +1,4 @@
+import PyCO2SYS as pyco2
 from PyCO2SYS import CO2System, system
 import networkx as nx
 
@@ -15,6 +16,8 @@
         # "total_to_sws_1atm",
         # "nbs_to_sws",
         "k_H2S_sws_1atm",
+        # 'gas_constant',
+        "factor_k_H2S",
     )
 )
 sys.plot_graph(