a revised version below:
using NonlinearSolve
using ModelingToolkit
# !!! used UNITS !!!
# eq/liter
# Pa
# constants
mmHG = 133.322 # Pa
Kc = 2.46e-11/mmHG
Ka = 3e-7
@variables HA A H HCO3
@parameters pCO2 Atot SID
eqs = [
H * A ~ Ka * HA
H * HCO3 ~ Kc * pCO2
Atot ~ HA + A
SID ~ HCO3 + A
]
@mtkbuild ns_stewart = NonlinearSystem(eqs)
display(ns_stewart)
println()
# Initial values (guesses)
guesses_stewart = (
HCO3 => 24.9e-3,
H => 3.95e-8,
A => 17.1e-3,
HA => 2e-3
)
ps = (
Atot => 19e-3,
SID => 41.7e-3,
pCO2 => 5.3e3
)
prob = NonlinearProblem(ns_stewart, guesses_stewart, ps)
sol = solve(prob)
@show sol
@show H = -log10(sol[1])
println("HCO3 = $(sol[2]*1000) meq/liter")
NEW Results:
julia> include("stewart1.jl")
Model ns_stewart with 2 equations
Unknowns (2):
H
HCO3
Parameters (3):
pCO2
Atot
SID
sol = [3.9275240698625016e-8, 0.02489948137605699]
H = -(log10(sol[1])) = 7.405881144702646
HCO3 = 24.89948137605699 meq/liter