Homec4science
Diffusion SpecMiCP / ReactMiCP (reactivetransport) Browse

specmicp/src/specmicp/equilibrium_data.cpp
583c1a0f5574reactivetransport

Last Change
Raw File

equilibrium_data.cpp
View Options

/*-------------------------------------------------------
- Module : specmicp
- File : equilibrium_data.cpp
- Author : Fabien Georget
Copyright (c) 2014, Fabien Georget, Princeton University
---------------------------------------------------------*/
#include "equilibrium_data.hpp"
#include "database/module.hpp"
namespace specmicp {
//! return log(IAP) for mineral m
double EquilibriumState::logIAP(int m)
{
double logiap = 0.0;
for (int i=1; i<m_data->nb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral(m, i);
return logiap;
}
//! return log(IAP) for mineral m
double EquilibriumState::logIAP_kinetic(int m)
{
double logiap = 0.0;
for (int i=1; i<m_data->nb_component; ++i) logiap += m_main_variables(i)*m_data->nu_mineral_kinetic(m, i);
return logiap;
}
//! Return total concentration of aqueous species for component i
double EquilibriumState::total_aqueous_concentration_component(int i)
{
double totconc = POW10(m_main_variables(i));
for (int j=0; j<m_data->nb_aqueous; ++j)
{
if (m_data->nu_aqueous(j,i) != 0)
totconc += m_data->nu_aqueous(j,i)*m_concentration_secondary(j);
}
return totconc;
}
void EquilibriumState::total_aqueous_concentrations(Eigen::VectorXd& totconc)
{
assert(totconc.rows() == m_data->nb_component);
totconc(0) = m_main_variables(0);
for (int i=1; i<m_data->nb_component; ++i) totconc(i) = total_aqueous_concentration_component(i);
}
double EquilibriumState::pH()
{
database::DatabaseModule namer(m_data);
double id = namer.component_label_to_id("H[+]");
if (id != database::DatabaseModule::no_species) return m_main_variables(id);
else
{
id = namer.aqueous_label_to_id("H[+]");
assert(id != database::DatabaseModule::no_species);
return -std::log10(m_concentration_secondary(id));
}
}
} // end namespace specmicp
c4science · Help