diff --git a/Cantera/src/thermo/Makefile.in b/Cantera/src/thermo/Makefile.in
index 55aa795db..12f9237a1 100644
--- a/Cantera/src/thermo/Makefile.in
+++ b/Cantera/src/thermo/Makefile.in
@@ -18,8 +18,10 @@ do_ranlib = @DO_RANLIB@
CXX_FLAGS = @CXXFLAGS@ $(CXX_OPT)
# Extended Cantera Thermodynamics Object Files
-CATHERMO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o
-CATHERMO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h
+CATHERMO_OBJ = SingleSpeciesTP.o StoichSubstanceSSTP.o \
+ MolalityVPSSTP.o VPStandardStateTP.o
+CATHERMO_H = SingleSpeciesTP.h StoichSubstanceSSTP.h \
+ MolalityVPSSTP.h VPStandardStateTP.h
CXX_INCLUDES = -I.. @CXX_INCLUDES@
LIB = @buildlib@/libcaThermo.a
diff --git a/Cantera/src/thermo/MolalityVPSSTP.cpp b/Cantera/src/thermo/MolalityVPSSTP.cpp
new file mode 100644
index 000000000..4be018f5a
--- /dev/null
+++ b/Cantera/src/thermo/MolalityVPSSTP.cpp
@@ -0,0 +1,505 @@
+/**
+ *
+ * @file MolalityVPSSTP.cpp
+ */
+/*
+ * Copywrite (2005) Sandia Corporation. Under the terms of
+ * Contract DE-AC04-94AL85000 with Sandia Corporation, the
+ * U.S. Government retains certain rights in this software.
+ */
+/*
+ * $Author$
+ * $Date$
+ * $Revision$
+ */
+#ifndef MAX
+#define MAX(x,y) (( (x) > (y) ) ? (x) : (y))
+#endif
+
+
+#include "MolalityVPSSTP.h"
+
+
+namespace Cantera {
+
+ /*
+ * Default constructor.
+ *
+ * This doesn't do much more than initialize constants with
+ * default values for water at 25C.
+ */
+ MolalityVPSSTP::MolalityVPSSTP() :
+ VPStandardStateTP(),
+ m_indexSolvent(0),
+ m_weightSolvent(18.0),
+ m_xmolSolventMIN(0.01),
+ m_Mnaught(18.0E-3)
+ {
+ }
+
+ /**
+ * Copy Constructor:
+ *
+ * Note this stuff will not work until the underlying phase
+ * has a working copy constructor
+ */
+ MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP &b) :
+ VPStandardStateTP(),
+ m_indexSolvent(b.m_indexSolvent),
+ m_xmolSolventMIN(b.m_xmolSolventMIN),
+ m_Mnaught(b.m_Mnaught),
+ m_molalities(b.m_molalities)
+ {
+ throw CanteraError("MolalityVPSSTP::operator=()",
+ "Not Implemented Fully");
+ *this = operator=(b);
+ }
+
+ /*
+ * operator=()
+ *
+ * Note this stuff will not work until the underlying phase
+ * has a working assignment operator
+ */
+ MolalityVPSSTP& MolalityVPSSTP::
+ operator=(const MolalityVPSSTP &b) {
+ if (&b != this) {
+ VPStandardStateTP::operator=(b);
+ m_indexSolvent = b.m_indexSolvent;
+ m_weightSolvent = b.m_weightSolvent;
+ m_xmolSolventMIN = b.m_xmolSolventMIN;
+ m_Mnaught = b.m_Mnaught;
+ m_molalities = b.m_molalities;
+ }
+ throw CanteraError("MolalityVPSSTP::operator=()",
+ "Not Implemented Fully");
+ return *this;
+ }
+
+ /**
+ *
+ * ~MolalityVPSSTP(): (virtual)
+ *
+ * Destructor: does nothing:
+ *
+ */
+ MolalityVPSSTP::~MolalityVPSSTP() {
+ }
+
+ /**
+ * This routine duplicates the current object and returns
+ * a pointer to ThermoPhase.
+ */
+ ThermoPhase*
+ MolalityVPSSTP::duplMyselfAsThermoPhase() {
+ MolalityVPSSTP* mtp = new MolalityVPSSTP(*this);
+ return (ThermoPhase *) mtp;
+ }
+
+ /*
+ * -------------- Utilities -------------------------------
+ */
+
+ /*
+ * setSolvent():
+ * Utilities for Solvent ID and Molality
+ * Here we also calculate and store the molecular weight
+ * of the solvent and the m_Mnaught parameter.
+ */
+ void MolalityVPSSTP::setSolvent(int k) {
+ if (k < 0 || k >= m_kk) {
+ throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
+ }
+ m_indexSolvent = k;
+ m_weightSolvent = molecularWeight(k);
+ m_Mnaught = m_weightSolvent / 1000.;
+ }
+
+ /*
+ * return the solvent id index number.
+ */
+ int MolalityVPSSTP::solventIndex() const {
+ return m_indexSolvent;
+ }
+
+ /*
+ * Sets the minimum mole fraction in the molality formulation
+ */
+ void MolalityVPSSTP::
+ setMoleFSolventMin(doublereal xmolSolventMIN) {
+ if (xmolSolventMIN <= 0.0) {
+ throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
+ } else if (xmolSolventMIN > 0.9) {
+ throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
+ }
+ m_xmolSolventMIN = xmolSolventMIN;
+ }
+
+ /*
+ * Returns the minimum mole fraction in the molality
+ * formulation.
+ */
+ doublereal MolalityVPSSTP::moleFSolventMin() const {
+ return m_xmolSolventMIN;
+ }
+
+ /**
+ * getMolalities():
+ * We calculate the vector of molalities of the species
+ * in the phase
+ * m_i = (n_i) / (1000 * M_o * n_o_p)
+ *
+ * where M_o is the molecular weight of the solvent
+ * n_o is the mole fraction of the solvent
+ * n_i is the mole fraction of the solute.
+ * n_o_p = max (n_o_min, n_o)
+ * n_o_min = minimum mole fraction of solvent allowed
+ * in the denominator.
+ */
+ void MolalityVPSSTP::getMolalities(doublereal * const molal) const {
+ getMoleFractions(molal);
+ double xmolSolvent = molal[m_indexSolvent];
+ if (xmolSolvent < m_xmolSolventMIN) {
+ xmolSolvent = m_xmolSolventMIN;
+ }
+ double denomInv = 1.0/
+ (m_Mnaught * xmolSolvent);
+ for (int k = 0; k < m_kk; k++) {
+ molal[k] *= denomInv;
+ }
+ for (int k = 0; k < m_kk; k++) {
+ m_molalities[k] = molal[k];
+ }
+ }
+
+ /**
+ * setMolalities():
+ * We are supplied with the molalities of all of the
+ * solute species. We then calculate the mole fractions of all
+ * species and update the ThermoPhase object.
+ *
+ * m_i = (n_i) / (W_o/1000 * n_o_p)
+ *
+ * where M_o is the molecular weight of the solvent
+ * n_o is the mole fraction of the solvent
+ * n_i is the mole fraction of the solute.
+ * n_o_p = max (n_o_min, n_o)
+ * n_o_min = minimum mole fraction of solvent allowed
+ * in the denominator.
+ */
+ void MolalityVPSSTP::setMolalities(const doublereal * const molal) {
+
+ double Lsum = 1.0 / m_Mnaught;
+ for (int k = 0; k < m_kk; k++) {
+ if (k != m_indexSolvent) {
+ m_molalities[k] = molal[k];
+ Lsum += molal[k];
+ }
+ }
+ double tmp = 1.0 / Lsum;
+ m_molalities[m_indexSolvent] = tmp / m_Mnaught;
+ double sum = m_molalities[m_indexSolvent];
+ for (int k = 0; k < m_kk; k++) {
+ if (k != m_indexSolvent) {
+ m_molalities[k] = tmp * molal[k];
+ sum += m_molalities[k];
+ }
+ }
+ if (sum != 1.0) {
+ tmp = 1.0 / sum;
+ for (int k = 0; k < m_kk; k++) {
+ m_molalities[k] *= tmp;
+ }
+ }
+ setMoleFractions(m_molalities.begin());
+ /*
+ * Essentially we don't trust the input: We calculate
+ * the molalities from the mole fractions that we
+ * just obtained.
+ */
+ getMolalities(m_molalities.begin());
+ }
+
+ /*
+ * setMolalitiesByName()
+ *
+ * This routine sets the molalities by name
+ * HKM -> Might need to be more complicated here, setting
+ * neutrals so that the existing mole fractions are
+ * preserved.
+ */
+ void MolalityVPSSTP::setMolalitiesByName(compositionMap& mMap) {
+ int kk = nSpecies();
+ doublereal x;
+ /*
+ * Get a vector of mole fractions
+ */
+ vector_fp mf(kk, 0.0);
+ getMoleFractions(mf.begin());
+ double xmolS = mf[m_indexSolvent];
+ double xmolSmin = max(xmolS, m_xmolSolventMIN);
+ compositionMap::iterator p;
+ for (int k = 0; k < kk; k++) {
+ p = mMap.find(speciesName(k));
+ if (p != mMap.end()) {
+ x = mMap[speciesName(k)];
+ if (x > 0.0) {
+ mf[k] = x * m_Mnaught * xmolSmin;
+ }
+ }
+ }
+ /*
+ * check charge neutrality
+ */
+ int largePos = -1;
+ double cPos = 0.0;
+ int largeNeg = -1;
+ double cNeg = 0.0;
+ double sum = 0.0;
+ for (int k = 0; k < kk; k++) {
+ double ch = charge(k);
+ if (mf[k] > 0.0) {
+ if (ch > 0.0) {
+ if (ch * mf[k] > cPos) {
+ largePos = k;
+ cPos = ch * mf[k];
+ }
+ }
+ if (ch < 0.0) {
+ if (fabs(ch) * mf[k] > cNeg) {
+ largeNeg = k;
+ cNeg = fabs(ch) * mf[k];
+ }
+ }
+ }
+ sum += mf[k] * ch;
+ }
+ if (sum != 0.0) {
+ if (sum > 0.0) {
+ if (cPos > sum) {
+ mf[largePos] -= sum / charge(largePos);
+ } else {
+ throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
+ "unbalanced charges");
+ }
+ } else {
+ if (cNeg > (-sum)) {
+ mf[largeNeg] -= (-sum) / fabs(charge(largeNeg));
+ } else {
+ throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
+ "unbalanced charges");
+ }
+ }
+
+ }
+ sum = 0.0;
+ for (int k = 0; k < kk; k++) {
+ sum += mf[k];
+ }
+ sum = 1.0/sum;
+ for (int k = 0; k < kk; k++) {
+ mf[k] *= sum;
+ }
+ setMoleFractions(mf.begin());
+ /*
+ * After we formally set the mole fractions, we
+ * calculate the molalities again and store it in
+ * this object.
+ */
+ getMolalities(m_molalities.begin());
+ }
+
+ /*
+ * setMolalitiesByNames()
+ *
+ * Set the molalities of the solutes by name
+ */
+ void MolalityVPSSTP::setMolalitiesByName(const string& x) {
+ compositionMap xx;
+ int kk = nSpecies();
+ for (int k = 0; k < kk; k++) {
+ xx[speciesName(k)] = -1.0;
+ }
+ parseCompString(x, xx);
+ setMolalitiesByName(xx);
+ }
+
+
+ /*
+ * Update the internal array that contains the molalities of the
+ * species.
+ */
+ void MolalityVPSSTP::updateMolalities() const {
+ getMolalities(m_molalities.begin());
+ }
+
+
+
+ /*
+ * ------------ Molar Thermodynamic Properties ----------------------
+ */
+
+
+ /*
+ * - Activities, Standard States, Activity Concentrations -----------
+ */
+
+ /**
+ * This method returns the activity convention.
+ * Currently, there are two activity conventions
+ * Molar-based activities
+ * Unit activity of species at either a hypothetical pure
+ * solution of the species or at a hypothetical
+ * pure ideal solution at infinite dilution
+ * cAC_CONVENTION_MOLAR 0
+ * - default
+ *
+ * Molality based activities
+ * (unit activity of solutes at a hypothetical 1 molal
+ * solution referenced to infinite dilution at all
+ * pressures and temperatures).
+ * (solvent is still on molar basis).
+ * cAC_CONVENTION_MOLALITY 1
+ *
+ * We set the convention to molality here.
+ */
+ int MolalityVPSSTP::activityConvention() const {
+ return cAC_CONVENTION_MOLALITY;
+ }
+
+ /**
+ * Get the array of non-dimensional activity coefficients at
+ * the current solution temperature, pressure, and
+ * solution concentration.
+ * These are mole fraction based activity coefficients. In this
+ * object, their calculation is based on translating the values
+ * of Molality based activity coefficients.
+ * See Denbigh p. 278 for a thorough discussion.
+ *
+ * Note, the solvent is treated differently. getMolalityActivityCoeff()
+ * returns the molar based solvent activity coefficient already.
+ * Therefore, we do not have to divide by x_s here.
+ */
+ void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const {
+ getMolalityActivityCoefficients(ac);
+ double xmolSolvent = moleFraction(m_indexSolvent);
+ if (xmolSolvent < m_xmolSolventMIN) {
+ xmolSolvent = m_xmolSolventMIN;
+ }
+ for (int k = 0; k < m_kk; k++) {
+ if (k != m_indexSolvent) {
+ ac[k] /= xmolSolvent;
+ }
+ }
+ }
+
+ /**
+ * osmotic coefficient:
+ *
+ * Calculate the osmotic coefficient of the solvent. Note there
+ * are lots of definitions of the osmotic coefficient floating
+ * around. We use the one defined in the Pitzer paper:
+ *
+ * Definition:
+ * - sum(m_i) * M0 * oc = ln(activity_solvent)
+ */
+ doublereal MolalityVPSSTP::osmoticCoefficient() const {
+ vector_fp act(m_kk);
+ getActivities(act.begin());
+ double sum = 0;
+ for (int k = 0; k < m_kk; k++) {
+ if (k != m_indexSolvent) {
+ sum += MAX(m_molalities[k], 0.0);
+ }
+ }
+ double oc = 1.0;
+ double lac = log(act[m_indexSolvent]);
+ if (sum > 1.0E-200) {
+ oc = - lac / (m_Mnaught * sum);
+ }
+ return oc;
+ }
+
+ /*
+ * ------------ Partial Molar Properties of the Solution ------------
+ */
+
+
+ doublereal MolalityVPSSTP::err(string msg) const {
+ throw CanteraError("MolalityVPSSTP","Base class method "
+ +msg+" called. Equation of state type: "+int2str(eosType()));
+ return 0;
+ }
+
+ /**
+ * Returns the units of the standard and general concentrations
+ * Note they have the same units, as their divisor is
+ * defined to be equal to the activity of the kth species
+ * in the solution, which is unitless.
+ *
+ * This routine is used in print out applications where the
+ * units are needed. Usually, MKS units are assumed throughout
+ * the program and in the XML input files.
+ *
+ * On return uA contains the powers of the units (MKS assumed)
+ * of the standard concentrations and generalized concentrations
+ * for the kth species.
+ *
+ * uA[0] = kmol units - default = 1
+ * uA[1] = m units - default = -nDim(), the number of spatial
+ * dimensions in the Phase class.
+ * uA[2] = kg units - default = 0;
+ * uA[3] = Pa(pressure) units - default = 0;
+ * uA[4] = Temperature units - default = 0;
+ * uA[5] = time units - default = 0
+ */
+ void MolalityVPSSTP::getUnitsStandardConc(double *uA, int k, int sizeUA) {
+ for (int i = 0; i < sizeUA; i++) {
+ if (i == 0) uA[0] = 1.0;
+ if (i == 1) uA[1] = -nDim();
+ if (i == 2) uA[2] = 0.0;
+ if (i == 3) uA[3] = 0.0;
+ if (i == 4) uA[4] = 0.0;
+ if (i == 5) uA[5] = 0.0;
+ }
+ }
+
+
+ /*
+ * Set the thermodynamic state.
+ */
+ void MolalityVPSSTP::setStateFromXML(const XML_Node& state) {
+ VPStandardStateTP::setStateFromXML(state);
+ string comp = getString(state,"soluteMolalities");
+ if (comp != "") {
+ setMolalitiesByName(comp);
+ }
+ if (state.hasChild("pressure")) {
+ double p = getFloat(state, "pressure", "pressure");
+ setPressure(p);
+ }
+ }
+
+ /**
+ * @internal Initialize. This method is provided to allow
+ * subclasses to perform any initialization required after all
+ * species have been added. For example, it might be used to
+ * resize internal work arrays that must have an entry for
+ * each species. The base class implementation does nothing,
+ * and subclasses that do not require initialization do not
+ * need to overload this method. When importing a CTML phase
+ * description, this method is called just prior to returning
+ * from function importPhase.
+ *
+ * @see importCTML.cpp
+ */
+ void MolalityVPSSTP::initThermo() {
+ VPStandardStateTP::initThermo();
+ m_molalities.resize(m_kk);
+ }
+
+}
+
+
+
+
diff --git a/Cantera/src/thermo/MolalityVPSSTP.h b/Cantera/src/thermo/MolalityVPSSTP.h
new file mode 100644
index 000000000..9995a9b72
--- /dev/null
+++ b/Cantera/src/thermo/MolalityVPSSTP.h
@@ -0,0 +1,455 @@
+/**
+ * @file MolalityVPSSTP.h
+ *
+ * Header file for a derived class of ThermoPhase that handles
+ * variable pressure standard state methods for calculating
+ * thermodynamic properties that are further based upon
+ * activities based on the molality scale.
+ * These include most of the
+ * methods for calculating liquid electrolyte thermodynamics.
+ */
+/*
+ * Copywrite (2005) Sandia Corporation. Under the terms of
+ * Contract DE-AC04-94AL85000 with Sandia Corporation, the
+ * U.S. Government retains certain rights in this software.
+ */
+/*
+ * $Author$
+ * $Date$
+ * $Revision$
+ */
+
+#ifndef CT_MOLALITYVPSSTP_H
+#define CT_MOLALITYVPSSTP_H
+
+#include "VPStandardStateTP.h"
+
+namespace Cantera {
+
+ /**
+ * @ingroup thermoprops
+ */
+
+ /**
+ * MolalityVPSSTP is a derived class of ThermoPhase that handles
+ * variable pressure standard state methods for calculating
+ * thermodynamic properties that are further based upon
+ * activities based on the molality scale.
+ * These include most of the
+ * methods for calculating liquid electrolyte thermodynamics.
+ */
+ class MolalityVPSSTP : public VPStandardStateTP {
+
+ public:
+
+ /// Constructors
+ MolalityVPSSTP();
+ MolalityVPSSTP(const MolalityVPSSTP &);
+ /// Assignment operator
+ MolalityVPSSTP& operator=(const MolalityVPSSTP&);
+
+ /// Destructor.
+ virtual ~MolalityVPSSTP();
+
+ /**
+ * Duplication routine for objects which inherit from
+ * ThermoPhase.
+ *
+ * This virtual routine can be used to duplicate thermophase objects
+ * inherited from ThermoPhase even if the application only has
+ * a pointer to ThermoPhase to work with.
+ */
+ virtual ThermoPhase *duplMyselfAsThermoPhase();
+
+ /**
+ *
+ * @name Utilities
+ * @{
+ */
+
+ /**
+ * Equation of state type flag. The ThermoPhase base class returns
+ * zero. Subclasses should define this to return a unique
+ * non-zero value. Known constants defined for this purpose are
+ * listed in mix_defs.h. The MolalityVPSSTP class also returns
+ * zero, as it is a non-complete class.
+ */
+ virtual int eosType() const { return 0; }
+
+
+ /**
+ * @}
+ * @name Molar Thermodynamic Properties
+ * @{
+ */
+
+
+ /**
+ * @}
+ * @name Utilities for Solvent ID and Molality
+ * @{
+ */
+
+ /**
+ * This routine sets the index number of the solvent for
+ * the phase.
+ *
+ * Note, having a solvent
+ * is a precursor to many things having to do with molality.
+ *
+ * @param k the solvent index number
+ */
+ void setSolvent(int k);
+
+ /**
+ * Sets the minimum mole fraction in the molality formulation.
+ * Note the molality formulation is singular in the limit that
+ * the solvent mole fraction goes to zero. Numerically, how
+ * this limit is treated and resolved is an ongoing issue within
+ * Cantera.
+ */
+ void setMoleFSolventMin(doublereal xmolSolventMIN);
+
+ /**
+ * Returns the solvent index.
+ */
+ int solventIndex() const;
+
+ /**
+ * Returns the minimum mole fraction in the molality
+ * formulation.
+ */
+ doublereal moleFSolventMin() const;
+
+ /**
+ * getMolalities()
+ * This function will return the molalities of the
+ * species.
+ *
+ */
+ void getMolalities(doublereal * const molal) const;
+
+
+ void setMolalities(const doublereal * const molal);
+ void setMolalitiesByName(compositionMap& xMap);
+ void setMolalitiesByName(const string &);
+ void updateMolalities() const;
+ /**
+ * @}
+ * @name Mechanical Properties
+ * @{
+ */
+
+ /**
+ * @}
+ * @name Potential Energy
+ *
+ * Species may have an additional potential energy due to the
+ * presence of external gravitation or electric fields. These
+ * methods allow specifying a potential energy for individual
+ * species.
+ * @{
+ */
+
+ /**
+ * @}
+ * @name Activities, Standard States, and Activity Concentrations
+ *
+ * The activity \f$a_k\f$ of a species in solution is
+ * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
+ * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
+ * the chemical potential at unit activity, which depends only
+ * on temperature and pressure.
+ * @{
+ */
+
+ /**
+ * This method returns the activity convention.
+ * Currently, there are two activity conventions
+ * Molar-based activities
+ * Unit activity of species at either a hypothetical pure
+ * solution of the species or at a hypothetical
+ * pure ideal solution at infinite dilution
+ * cAC_CONVENTION_MOLAR 0
+ * - default
+ *
+ * Molality based acvtivities
+ * (unit activity of solutes at a hypothetical 1 molal
+ * solution referenced to infinite dilution at all
+ * pressures and temperatures).
+ * cAC_CONVENTION_MOLALITY 1
+ *
+ * We set the convention to molality here.
+ */
+ int activityConvention() const;
+
+ /**
+ * This method returns an array of generalized concentrations
+ * \f$ C_k\f$ that are defined such that
+ * \f$ a_k = C_k / C^0_k, \f$ where \f$ C^0_k \f$
+ * is a standard concentration
+ * defined below. These generalized concentrations are used
+ * by kinetics manager classes to compute the forward and
+ * reverse rates of elementary reactions.
+ *
+ * @param c Array of generalized concentrations. The
+ * units depend upon the implementation of the
+ * reaction rate expressions within the phase.
+ */
+ virtual void getActivityConcentrations(doublereal* c) const {
+ err("getActivityConcentrations");
+ }
+
+ /**
+ * The standard concentration \f$ C^0_k \f$ used to normalize
+ * the generalized concentration. In many cases, this quantity
+ * will be the same for all species in a phase - for example,
+ * for an ideal gas \f$ C^0_k = P/\hat R T \f$. For this
+ * reason, this method returns a single value, instead of an
+ * array. However, for phases in which the standard
+ * concentration is species-specific (e.g. surface species of
+ * different sizes), this method may be called with an
+ * optional parameter indicating the species.
+ */
+ virtual doublereal standardConcentration(int k=0) const {
+ err("standardConcentration");
+ return -1.0;
+ }
+
+ /**
+ * Returns the natural logarithm of the standard
+ * concentration of the kth species
+ */
+ virtual doublereal logStandardConc(int k=0) const {
+ err("logStandardConc");
+ return -1.0;
+ }
+
+ /**
+ * Returns the units of the standard and generalized
+ * concentrations Note they have the same units, as their
+ * ratio is defined to be equal to the activity of the kth
+ * species in the solution, which is unitless.
+ *
+ * This routine is used in print out applications where the
+ * units are needed. Usually, MKS units are assumed throughout
+ * the program and in the XML input files.
+ *
+ * uA[0] = kmol units - default = 1
+ * uA[1] = m units - default = -nDim(), the number of spatial
+ * dimensions in the Phase class.
+ * uA[2] = kg units - default = 0;
+ * uA[3] = Pa(pressure) units - default = 0;
+ * uA[4] = Temperature units - default = 0;
+ * uA[5] = time units - default = 0
+ */
+ virtual void getUnitsStandardConc(double *uA, int k = 0,
+ int sizeUA = 6);
+
+ /**
+ * Get the array of non-dimensional activities (molality
+ * based for this class and classes that derive from it) at
+ * the current solution temperature, pressure, and
+ * solution concentration.
+ */
+ virtual void getActivities(doublereal* ac) const {
+ err("getActivities");
+ }
+
+ /**
+ * Get the array of non-dimensional activity coefficients at
+ * the current solution temperature, pressure, and
+ * solution concentration.
+ * These are mole fraction based activity coefficients. In this
+ * object, their calculation is based on translating the values
+ * of Molality based activity coefficients.
+ * See Denbigh p. 278 for a thorough discussion
+ */
+ void getActivityCoefficients(doublereal* ac) const;
+
+ /**
+ * Get the array of non-dimensional molality based
+ * activity coefficients at the current solution temperature,
+ * pressure, and solution concentration.
+ * See Denbigh p. 278 for a thorough discussion
+ */
+ virtual void getMolalityActivityCoefficients(doublereal *acMolality)
+ const {
+ err("getMolalityActivityCoefficients");
+ }
+
+ /**
+ * Calculate the osmotic coefficient
+ * units = dimensionless
+ */
+ virtual double osmoticCoefficient() const;
+
+ //@}
+ /// @name Partial Molar Properties of the Solution
+ //@{
+
+
+ /**
+ * Get the species electrochemical potentials.
+ * These are partial molar quantities.
+ * This method adds a term \f$ Fz_k \phi_k \f$ to the
+ * to each chemical potential.
+ *
+ * Units: J/kmol
+ */
+ void getElectrochemPotentials(doublereal* mu) const {
+ getChemPotentials(mu);
+ double ve = Faraday * electricPotential();
+ for (int k = 0; k < m_kk; k++) {
+ mu[k] += ve*charge(k);
+ }
+ }
+
+
+ //@}
+ /// @name Properties of the Standard State of the Species in the Solution
+ //@{
+
+
+
+ //@}
+ /// @name Thermodynamic Values for the Species Reference States
+ //@{
+
+
+ ///////////////////////////////////////////////////////
+ //
+ // The methods below are not virtual, and should not
+ // be overloaded.
+ //
+ //////////////////////////////////////////////////////
+
+ /**
+ * @name Specific Properties
+ * @{
+ */
+
+
+ /**
+ * @name Setting the State
+ *
+ * These methods set all or part of the thermodynamic
+ * state.
+ * @{
+ */
+
+ //@}
+
+ /**
+ * @name Chemical Equilibrium
+ * Routines that implement the Chemical equilibrium capability
+ * for a single phase, based on the element-potential method.
+ * @{
+ */
+
+ /**
+ * This method is used by the ChemEquil element-potential
+ * based equilibrium solver.
+ * It sets the state such that the chemical potentials of the
+ * species within the current phase satisfy
+ * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
+ * \left(\frac{\lambda_m} {\hat R T}\right) \f] where
+ * \f$ \lambda_m \f$ is the element potential of element m. The
+ * temperature is unchanged. Any phase (ideal or not) that
+ * implements this method can be equilibrated by ChemEquil.
+ */
+ virtual void setToEquilState(const doublereal* lambda_RT) {
+ err("setToEquilState");
+ }
+
+ // called by function 'equilibrate' in ChemEquil.h to transfer
+ // the element potentials to this object
+ void setElementPotentials(const vector_fp& lambda) {
+ m_lambda = lambda;
+ }
+
+ void getElementPotentials(doublereal* lambda) {
+ copy(m_lambda.begin(), m_lambda.end(), lambda);
+ }
+
+ //@}
+
+
+ /**
+ * Set equation of state parameter values from XML
+ * entries. This method is called by function importPhase in
+ * file importCTML.cpp when processing a phase definition in
+ * an input file. It should be overloaded in subclasses to set
+ * any parameters that are specific to that particular phase
+ * model.
+ *
+ * The MolalityVPSSTP object defines a new method for setting
+ * the concentrations of a phase. The new method is defined by a
+ * block called "soluteMolalities". If this block
+ * is found, the concentrations within that phase are
+ * set to the "name":"molalities pairs found within that
+ * XML block. The solvent concentration is then set
+ * to everything else.
+ *
+ * @param eosdata An XML_Node object corresponding to
+ * the "thermo" entry for this phase in the input file.
+ *
+ */
+ virtual void setStateFromXML(const XML_Node& state);
+
+ /// The following methods are used in the process of constructing
+ /// the phase and setting its parameters from a specification in an
+ /// input file. They are not normally used in application programs.
+ /// To see how they are used, see files importCTML.cpp and
+ /// ThermoFactory.cpp.
+
+ /**
+ * @internal Initialize. This method is provided to allow
+ * subclasses to perform any initialization required after all
+ * species have been added. For example, it might be used to
+ * resize internal work arrays that must have an entry for
+ * each species. The base class implementation does nothing,
+ * and subclasses that do not require initialization do not
+ * need to overload this method. When importing a CTML phase
+ * description, this method is called just prior to returning
+ * from function importPhase.
+ *
+ * @see importCTML.cpp
+ */
+ virtual void initThermo();
+
+ protected:
+
+ int m_indexSolvent;
+ doublereal m_weightSolvent;
+ /*
+ * In any molality implementation, it makes sense to have
+ * a minimum solvent mole fraction requirement, since the
+ * implementation becomes singular in the xmolSolvent=0
+ * limit. The default is to set it to 0.01.
+ * We then modify the molality definition to ensure that
+ * molal_solvent = 0 when xmol_solvent = 0.
+ */
+ doublereal m_xmolSolventMIN;
+ /*
+ * This is the multiplication factor that goes inside
+ * log expressions involving the molalities of species.
+ * Its equal to Wt_0 / 1000.
+ * where Wt_0 = weight of solvent (kg/kmol)
+ */
+ doublereal m_Mnaught;
+
+ mutable vector_fp m_molalities;
+ private:
+ doublereal err(string msg) const;
+
+ };
+
+}
+
+#endif
+
+
+
+
+
diff --git a/Cantera/src/thermo/StoichSubstanceSSTP.h b/Cantera/src/thermo/StoichSubstanceSSTP.h
index bcec4cb3f..60850ed18 100644
--- a/Cantera/src/thermo/StoichSubstanceSSTP.h
+++ b/Cantera/src/thermo/StoichSubstanceSSTP.h
@@ -108,7 +108,6 @@ namespace Cantera {
*/
virtual doublereal thermalExpansionCoeff() const ;
- //@}
/**
* @}
diff --git a/Cantera/src/thermo/VPStandardStateTP.cpp b/Cantera/src/thermo/VPStandardStateTP.cpp
new file mode 100644
index 000000000..b0a285d40
--- /dev/null
+++ b/Cantera/src/thermo/VPStandardStateTP.cpp
@@ -0,0 +1,300 @@
+/**
+ *
+ * @file VPStandardStateTP.cpp
+ */
+/*
+ * Copywrite (2005) Sandia Corporation. Under the terms of
+ * Contract DE-AC04-94AL85000 with Sandia Corporation, the
+ * U.S. Government retains certain rights in this software.
+ */
+/*
+ * $Author$
+ * $Date$
+ * $Revision$
+ */
+
+// turn off warnings under Windows
+#ifdef WIN32
+#pragma warning(disable:4786)
+#pragma warning(disable:4503)
+#endif
+
+#include "VPStandardStateTP.h"
+
+
+namespace Cantera {
+
+ /*
+ * Default constructor
+ */
+ VPStandardStateTP::VPStandardStateTP() :
+ ThermoPhase(),
+ m_tlast(-1.0)
+ {
+ }
+
+ /*
+ * Copy Constructor:
+ *
+ * Note this stuff will not work until the underlying phase
+ * has a working copy constructor.
+ *
+ * The copy constructor just calls the assignment operator
+ * to do the heavy lifting.
+ */
+ VPStandardStateTP::VPStandardStateTP(const VPStandardStateTP &b) :
+ ThermoPhase(),
+ m_tlast(-1.0)
+ {
+ *this = b;
+ }
+
+ /*
+ * operator=()
+ *
+ * Note this stuff will not work until the underlying phase
+ * has a working assignment operator
+ */
+ VPStandardStateTP& VPStandardStateTP::
+ operator=(const VPStandardStateTP &b) {
+ if (&b != this) {
+ /*
+ * Mostly, this is a passthrough to the underlying
+ * assignment operator for the ThermoPhae parent object.
+ */
+ ThermoPhase::operator=(b);
+ /*
+ * However, we have to handle data that we own.
+ */
+ m_tlast = b.m_tlast;
+ m_h0_RT = b.m_h0_RT;
+ m_cp0_R = b.m_cp0_R;
+ m_g0_RT = b.m_g0_RT;
+ m_s0_R = b.m_s0_R;
+ }
+ return *this;
+ }
+
+ /*
+ * ~VPStandardStateTP(): (virtual)
+ *
+ * This destructor does nothing. All of the owned objects
+ * handle themselves.
+ */
+ VPStandardStateTP::~VPStandardStateTP() {
+ }
+
+ /*
+ * Duplication function.
+ * This calls the copy constructor for this object.
+ */
+ ThermoPhase* VPStandardStateTP::duplMyselfAsThermoPhase() {
+ VPStandardStateTP* vptp = new VPStandardStateTP(*this);
+ return (ThermoPhase *) vptp;
+ }
+
+ /*
+ * -------------- Utilities -------------------------------
+ */
+
+
+ /*
+ * ------------Molar Thermodynamic Properties -------------------------
+ */
+
+
+ doublereal VPStandardStateTP::err(string msg) const {
+ throw CanteraError("VPStandardStateTP","Base class method "
+ +msg+" called. Equation of state type: "+int2str(eosType()));
+ return 0;
+ }
+
+ /**
+ * Returns the units of the standard and general concentrations
+ * Note they have the same units, as their divisor is
+ * defined to be equal to the activity of the kth species
+ * in the solution, which is unitless.
+ *
+ * This routine is used in print out applications where the
+ * units are needed. Usually, MKS units are assumed throughout
+ * the program and in the XML input files.
+ *
+ * On return uA contains the powers of the units (MKS assumed)
+ * of the standard concentrations and generalized concentrations
+ * for the kth species.
+ *
+ * uA[0] = kmol units - default = 1
+ * uA[1] = m units - default = -nDim(), the number of spatial
+ * dimensions in the Phase class.
+ * uA[2] = kg units - default = 0;
+ * uA[3] = Pa(pressure) units - default = 0;
+ * uA[4] = Temperature units - default = 0;
+ * uA[5] = time units - default = 0
+ */
+ void VPStandardStateTP::
+ getUnitsStandardConc(double *uA, int k, int sizeUA) {
+ for (int i = 0; i < sizeUA; i++) {
+ if (i == 0) uA[0] = 1.0;
+ if (i == 1) uA[1] = -nDim();
+ if (i == 2) uA[2] = 0.0;
+ if (i == 3) uA[3] = 0.0;
+ if (i == 4) uA[4] = 0.0;
+ if (i == 5) uA[5] = 0.0;
+ }
+ }
+
+ /*
+ * ---- Partial Molar Properties of the Solution -----------------
+ */
+
+ /**
+ * Get the array of non-dimensional species chemical potentials
+ * These are partial molar Gibbs free energies.
+ * \f$ \mu_k / \hat R T \f$.
+ * Units: unitless
+ *
+ * We close the loop on this function, here, calling
+ * getChemPotentials() and then dividing by RT.
+ */
+ void VPStandardStateTP::getChemPotentials_RT(doublereal* muRT) const{
+ getChemPotentials(muRT);
+ doublereal invRT = 1.0 / _RT();
+ for (int k = 0; k < m_kk; k++) {
+ muRT[k] *= invRT;
+ }
+ }
+
+ /*
+ * ----- Thermodynamic Values for the Species Reference States ----
+ */
+
+ /**
+ * Returns the vector of nondimensional
+ * enthalpies of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ */
+ void VPStandardStateTP::getEnthalpy_RT_ref(doublereal *hrt) const {
+ /*
+ * Call the function that makes sure the local copy of
+ * the species reference thermo functions are up to date
+ * for the current temperature.
+ */
+ _updateRefStateThermo();
+ /*
+ * Copy the enthalpy function into return vector.
+ */
+ copy(m_h0_RT.begin(), m_h0_RT.end(), hrt);
+ }
+
+ /**
+ * Returns the vector of nondimensional
+ * enthalpies of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ */
+ void VPStandardStateTP::getGibbs_RT_ref(doublereal *grt) const {
+ /*
+ * Call the function that makes sure the local copy of
+ * the species reference thermo functions are up to date
+ * for the current temperature.
+ */
+ _updateRefStateThermo();
+ /*
+ * Copy the gibbs function into return vector.
+ */
+ copy(m_g0_RT.begin(), m_g0_RT.end(), grt);
+ }
+
+ /**
+ * Returns the vector of the
+ * gibbs function of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ * units = J/kmol
+ *
+ * This is filled in here so that derived classes don't have to
+ * take care of it.
+ */
+ void VPStandardStateTP::getGibbs_ref(doublereal *g) const {
+ getGibbs_RT_ref(g);
+ double RT = _RT();
+ for (int k = 0; k < m_kk; k++) {
+ g[k] *= RT;
+ }
+ }
+
+ /**
+ * Returns the vector of nondimensional
+ * entropies of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ */
+ void VPStandardStateTP::getEntropy_R_ref(doublereal *er) const {
+ /*
+ * Call the function that makes sure the local copy of
+ * the species reference thermo functions are up to date
+ * for the current temperature.
+ */
+ _updateRefStateThermo();
+ /*
+ * Copy the gibbs function into return vector.
+ */
+ copy(m_s0_R.begin(), m_s0_R.end(), er);
+ }
+
+ /**
+ * Returns the vector of nondimensional
+ * constant pressure heat capacities of the reference state
+ * at the current temperature of the solution
+ * and reference pressure for the species.
+ */
+ void VPStandardStateTP::getCp_R_ref(doublereal *cpr) const {
+ /*
+ * Call the function that makes sure the local copy of
+ * the species reference thermo functions are up to date
+ * for the current temperature.
+ */
+ _updateRefStateThermo();
+ /*
+ * Copy the gibbs function into return vector.
+ */
+ copy(m_cp0_R.begin(), m_cp0_R.end(), cpr);
+ }
+
+ /**
+ * Perform initializations after all species have been
+ * added.
+ */
+ void VPStandardStateTP::initThermo() {
+ ThermoPhase::initThermo();
+ m_kk = nSpecies();
+ int leng = m_kk;
+ m_h0_RT.resize(leng);
+ m_g0_RT.resize(leng);
+ m_cp0_R.resize(leng);
+ m_s0_R.resize(leng);
+ }
+
+ /**
+ * void _updateRefStateThermo() (private, const)
+ *
+ * This function gets called for every call to functions in this
+ * class. It checks to see whether the temperature has changed and
+ * thus the reference thermodynamics functions for all of the species
+ * must be recalculated.
+ * If the temperature has changed, the species thermo manager is called
+ * to recalculate G, Cp, H, and S at the current temperature.
+ */
+ void VPStandardStateTP::_updateRefStateThermo() const {
+ doublereal tnow = temperature();
+ if (m_tlast != tnow) {
+ m_spthermo->update(tnow, m_cp0_R.begin(), m_h0_RT.begin(),
+ m_s0_R.begin());
+ m_tlast = tnow;
+ for (int k = 0; k < m_kk; k++) {
+ m_g0_RT[k] = m_h0_RT[k] - m_s0_R[k];
+ }
+ }
+ }
+}
+
+
+
+
diff --git a/Cantera/src/thermo/VPStandardStateTP.h b/Cantera/src/thermo/VPStandardStateTP.h
new file mode 100644
index 000000000..735361192
--- /dev/null
+++ b/Cantera/src/thermo/VPStandardStateTP.h
@@ -0,0 +1,439 @@
+/**
+ * @file VPStandardStateTP.h
+ *
+ * Header file for a derived class of ThermoPhase that handles
+ * variable pressure standard state methods for calculating
+ * thermodynamic properties. These include most of the
+ * methods for calculating liquid electrolyte thermodynamics.
+ */
+/*
+ * Copywrite (2005) Sandia Corporation. Under the terms of
+ * Contract DE-AC04-94AL85000 with Sandia Corporation, the
+ * U.S. Government retains certain rights in this software.
+ */
+/*
+ * $Author$
+ * $Date$
+ * $Revision$
+ */
+
+#ifndef CT_VPSTANDARDSTATETP_H
+#define CT_VPSTANDARDSTATETP_H
+
+#include "ThermoPhase.h"
+
+namespace Cantera {
+
+ class XML_Node;
+
+ /**
+ * @ingroup thermoprops
+ *
+ * This is a filter class for ThermoPhase that implements
+ * a variable pressure standard state for ThermoPhase objects.
+ *
+ * In addition support for the molality unit scale is provided.
+ *
+ * Currently, it really is just a shell. The ThermoPhase object
+ * itself is based around the general concepts of
+ * VPStandardStateTP. Therefore, there really isn't much going
+ * on here.
+ * However, this may change. The ThermoPhase object itself
+ * could change. Additionally, this object may revolve around
+ * the molality unit scale in the near future. We will have to see
+ * how things fare.
+ */
+
+ class VPStandardStateTP : public ThermoPhase {
+
+ public:
+
+ /// Constructor.
+ VPStandardStateTP();
+
+ /// Copy Constructor.
+ VPStandardStateTP(const VPStandardStateTP &);
+
+ /// Assignment operator
+ VPStandardStateTP& operator=(const VPStandardStateTP &);
+ /// Destructor.
+ virtual ~VPStandardStateTP();
+
+ /*
+ * Duplication routine
+ */
+ virtual ThermoPhase *duplMyselfAsThermoPhase();
+
+ /**
+ *
+ * @name Utilities
+ * @{
+ */
+
+ /**
+ * Equation of state type flag. The base class returns
+ * zero. Subclasses should define this to return a unique
+ * non-zero value. Constants defined for this purpose are
+ * listed in mix_defs.h.
+ */
+ virtual int eosType() const { return 0; }
+
+
+ /**
+ * @}
+ * @name Molar Thermodynamic Properties of the Solution
+ * @{
+ */
+
+ /*
+ * These are handled by inherited objects. At this level,
+ * this pass-through routine doesn't add anything to the
+ * ThermoPhase description.
+ */
+
+
+ /**
+ * @}
+ * @name Mechanical Properties
+ * @{
+ */
+
+ /*
+ * These are handled by inherited objects. At this level,
+ * this pass-through routine doesn't add anything to the
+ * ThermoPhase description.
+ */
+
+ /**
+ * @}
+ * @name Electric Potential
+ *
+ * The phase may be at some non-zero electrical
+ * potential. These methods set or get the value of the
+ * electric potential.
+ * @{
+ */
+
+ /*
+ * These are handled by inherited objects. At this level,
+ * this pass-through routine doesn't add anything to the
+ * ThermoPhase description.
+ */
+
+ /**
+ * @}
+ * @name Activities and Activity Concentrations
+ *
+ * The activity \f$a_k\f$ of a species in solution is
+ * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
+ * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T)\f$ is
+ * the chemical potential at unit activity, which depends only
+ * on temperature.
+ * @{
+ */
+
+
+ /**
+ * Returns the units of the standard and generalized
+ * concentrations Note they have the same units, as their
+ * ratio is defined to be equal to the activity of the kth
+ * species in the solution, which is unitless.
+ *
+ * This routine is used in print out applications where the
+ * units are needed. Usually, MKS units are assumed throughout
+ * the program and in the XML input files.
+ *
+ * uA[0] = kmol units - default = 1
+ * uA[1] = m units - default = -nDim(), the number of spatial
+ * dimensions in the Phase class.
+ * uA[2] = kg units - default = 0;
+ * uA[3] = Pa(pressure) units - default = 0;
+ * uA[4] = Temperature units - default = 0;
+ * uA[5] = time units - default = 0
+ */
+ virtual void getUnitsStandardConc(double *uA, int k = 0,
+ int sizeUA = 6);
+
+ //@}
+ /// @name Partial Molar Properties of the Solution
+ //@{
+
+ /**
+ * Get the array of non-dimensional species chemical potentials
+ * These are partial molar Gibbs free energies.
+ * \f$ \mu_k / \hat R T \f$.
+ * Units: unitless
+ *
+ * We close the loop on this function, here, calling
+ * getChemPotentials() and then dividing by RT.
+ */
+ virtual void getChemPotentials_RT(doublereal* mu) const;
+
+
+ //@}
+ /// @name Properties of the Standard State of the Species in the Solution
+ //@{
+
+ /*
+ * These are handled by inherited objects. At this level,
+ * this pass-through routine doesn't add anything to the
+ * ThermoPhase description.
+ *
+ * However, we assume these methods exist for inherited objects.
+ * Therefore, we will bring the error routines up to this object
+ */
+
+ /**
+ * Get the array of chemical potentials at unit activity.
+ * These
+ * are the standard state chemical potentials \f$ \mu^0_k(T,P)
+ * \f$.. The values are evaluated at the current
+ * temperature and pressure.
+ */
+ virtual void getStandardChemPotentials(doublereal* mu) const {
+ err("getStandardChemPotentials");
+ }
+
+ /**
+ * Get the nondimensional Enthalpy functions for the species
+ * at their standard states at the current
+ * T and P of the solution.
+ */
+ virtual void getEnthalpy_RT(doublereal* hrt) const {
+ err("getEnthalpy_RT");
+ }
+
+ /**
+ * Get the array of nondimensional Enthalpy functions for the
+ * standard state species
+ * at the current T and P of the solution.
+ */
+ virtual void getEntropy_R(doublereal* sr) const {
+ err("getEntropy_R");
+ }
+
+ /**
+ * Get the nondimensional Gibbs functions for the species
+ * at their standard states of solution at the current T and P
+ * of the solution.
+ */
+ virtual void getGibbs_RT(doublereal* grt) const {
+ err("getGibbs_RT");
+ }
+
+ /**
+ * Get the nondimensional Gibbs functions for the standard
+ * state of the species at the current T and P.
+ */
+ virtual void getPureGibbs(doublereal* gpure) const {
+ err("getPureGibbs");
+ }
+
+ /**
+ * Returns the vector of nondimensional
+ * internal Energies of the standard state at the current temperature
+ * and pressure of the solution for each species.
+ */
+ virtual void getIntEnergy_RT(doublereal *urt) const {
+ err("getIntEnergy_RT");
+ }
+
+ /**
+ * Get the nondimensional Heat Capacities at constant
+ * pressure for the standard state of the species
+ * at the current T and P.
+ */
+ virtual void getCp_R(doublereal* cpr) const {
+ err("getCp_R");
+ }
+
+ /**
+ * Get the molar volumes of each species in their standard
+ * states at the current
+ * T and P of the solution.
+ * units = m^3 / kmol
+ */
+ virtual void getStandardVolumes(doublereal *vol) const {
+ err("getStandardVolumes");
+ }
+
+ //@}
+ /// @name Thermodynamic Values for the Species Reference States --------------------
+ //@{
+
+ /**
+ * Returns the vector of nondimensional
+ * enthalpies of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ */
+ virtual void getEnthalpy_RT_ref(doublereal *hrt) const;
+
+ /**
+ * Returns the vector of nondimensional
+ * enthalpies of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ */
+ virtual void getGibbs_RT_ref(doublereal *grt) const;
+
+ /**
+ * Returns the vector of the
+ * gibbs function of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ * units = J/kmol
+ */
+ virtual void getGibbs_ref(doublereal *g) const;
+
+ /**
+ * Returns the vector of nondimensional
+ * entropies of the reference state at the current temperature
+ * of the solution and the reference pressure for the species.
+ */
+ virtual void getEntropy_R_ref(doublereal *er) const;
+
+ /**
+ * Returns the vector of nondimensional
+ * constant pressure heat capacities of the reference state
+ * at the current temperature of the solution
+ * and reference pressure for the species.
+ */
+ virtual void getCp_R_ref(doublereal *cprt) const;
+
+ ///////////////////////////////////////////////////////
+ //
+ // The methods below are not virtual, and should not
+ // be overloaded.
+ //
+ //////////////////////////////////////////////////////
+
+ /**
+ * @name Specific Properties
+ * @{
+ */
+
+
+ /**
+ * @name Setting the State
+ *
+ * These methods set all or part of the thermodynamic
+ * state.
+ * @{
+ */
+
+ //@}
+
+ /**
+ * @name Chemical Equilibrium
+ * Chemical equilibrium.
+ * @{
+ */
+
+ //@}
+
+
+ /**
+ * Set equation of state parameter values from XML
+ * entries. This method is called by function importPhase in
+ * file importCTML.cpp when processing a phase definition in
+ * an input file. It should be overloaded in subclasses to set
+ * any parameters that are specific to that particular phase
+ * model.
+ *
+ * @param eosdata An XML_Node object corresponding to
+ * the "thermo" entry for this phase in the input file.
+ */
+ virtual void setParametersFromXML(const XML_Node& eosdata) {}
+
+
+ //---------------------------------------------------------
+ /// @name Critical state properties.
+ /// These methods are only implemented by some subclasses.
+
+ //@{
+
+ //@}
+
+ /// @name Saturation properties.
+ /// These methods are only implemented by subclasses that
+ /// implement full liquid-vapor equations of state.
+ ///
+
+
+ //@}
+
+ /// The following methods are used in the process of constructing
+ /// the phase and setting its parameters from a specification in an
+ /// input file. They are not normally used in application programs.
+ /// To see how they are used, see files importCTML.cpp and
+ /// ThermoFactory.cpp.
+
+ /**
+ * @internal Initialize. This method is provided to allow
+ * subclasses to perform any initialization required after all
+ * species have been added. For example, it might be used to
+ * resize internal work arrays that must have an entry for
+ * each species. The base class implementation does nothing,
+ * and subclasses that do not require initialization do not
+ * need to overload this method. When importing a CTML phase
+ * description, this method is called just prior to returning
+ * from function importPhase.
+ *
+ * @see importCTML.cpp
+ */
+ virtual void initThermo();
+
+ protected:
+ /*
+ * The last temperature at which the reference thermodynamic
+ * properties were calculated at.
+ */
+ mutable doublereal m_tlast;
+ /**
+ * Vector containing the species reference enthalpies at T = m_tlast
+ */
+ mutable vector_fp m_h0_RT;
+
+ /**
+ * Vector containing the species reference constant pressure
+ * heat capacities at T = m_tlast
+ */
+ mutable vector_fp m_cp0_R;
+
+ /**
+ * Vector containing the species reference Gibbs functions
+ * at T = m_tlast
+ */
+ mutable vector_fp m_g0_RT;
+
+ /**
+ * Vector containing the species reference entropies
+ * at T = m_tlast
+ */
+ mutable vector_fp m_s0_R;
+
+ private:
+
+ /**
+ * VPStandardStateTP has its own err routine
+ *
+ */
+ doublereal err(string msg) const;
+
+ /**
+ * This function gets called for every call to functions in this
+ * class. It checks to see whether the temperature has changed and
+ * thus the reference thermodynamics functions for all of the species
+ * must be recalculated.
+ * If the temperature has changed, the species thermo manager is called
+ * to recalculate G, Cp, H, and S at the current temperature.
+ */
+ void _updateRefStateThermo() const;
+ };
+
+}
+
+#endif
+
+
+
+
+