HMWSoln.h

Go to the documentation of this file.
00001 /**
00002  *  @file HMWSoln.h
00003  *    Headers for the %HMWSoln ThermoPhase object, which models concentrated
00004  *    electrolyte solutions
00005  *    (see \ref thermoprops and \link Cantera::HMWSoln HMWSoln \endlink) .
00006  *
00007  * Class %HMWSoln represents a concentrated liquid electrolyte phase which
00008  * obeys the Pitzer formulation for nonideality using molality-based
00009  * standard states.
00010  */
00011 /*
00012  * Copywrite (2006) Sandia Corporation. Under the terms of 
00013  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
00014  * U.S. Government retains certain rights in this software.
00015  */
00016 /*
00017  *  $Id: HMWSoln.h 386 2010-01-17 17:55:28Z hkmoffa $
00018  */
00019 
00020 #ifndef CT_HMWSOLN_H
00021 #define CT_HMWSOLN_H
00022 
00023 #include "MolalityVPSSTP.h"
00024 #include "electrolytes.h"
00025 
00026 namespace Cantera {
00027 
00028 
00029   /**
00030    * Major Parameters:
00031    *   The form of the Pitzer expression refers to the 
00032    *   form of the Gibbs free energy expression. The temperature
00033    *   dependence of the Pitzer coefficients are handled by
00034    *   another parameter.
00035    *
00036    * m_formPitzer = Form of the Pitzer expression
00037    *
00038    *  PITZERFORM_BASE = 0
00039    *
00040    *   Only one form is supported atm. This parameter is included for
00041    *   future expansion.
00042    *  
00043    */
00044 #define PITZERFORM_BASE 0
00045 
00046 
00047   /*!
00048    * @name Temperature Dependence of the Pitzer Coefficients
00049    *
00050    *  Note, the temperature dependence of the
00051    * Gibbs free energy also depends on the temperature dependence
00052    * of the standard state and the temperature dependence of the
00053    * Debye-Huckel constant, which includes the dielectric constant
00054    * and the density. Therefore, this expression defines only part
00055    * of the temperature dependence for the mixture thermodynamic
00056    * functions.
00057    *
00058    *  PITZER_TEMP_CONSTANT
00059    *     All coefficients are considered constant wrt temperature
00060    *  PITZER_TEMP_LINEAR
00061    *     All coefficients are assumed to have a linear dependence
00062    *     wrt to temperature.
00063    *  PITZER_TEMP_COMPLEX1
00064    *     All coefficients are assumed to have a complex functional
00065    *     based dependence wrt temperature;  See:
00066    *    (Silvester, Pitzer, J. Phys. Chem. 81, 19 1822 (1977)).
00067    *
00068    *       beta0 = q0 + q3(1/T - 1/Tr) + q4(ln(T/Tr)) +
00069    *               q1(T - Tr) + q2(T**2 - Tr**2)
00070    * 
00071    */
00072   //@{
00073 #define PITZER_TEMP_CONSTANT   0
00074 #define PITZER_TEMP_LINEAR     1
00075 #define PITZER_TEMP_COMPLEX1   2
00076   //@}
00077 
00078   /*
00079    *  @name ways to calculate the value of A_Debye
00080    *
00081    *  These defines determine the way A_Debye is calculated
00082    */
00083  //@{
00084 #define    A_DEBYE_CONST  0
00085 #define    A_DEBYE_WATER  1
00086   //@}
00087 
00088   class WaterProps;
00089   class PDSS_Water;
00090 
00091   /**
00092    * Class %HMWSoln represents a dilute or concentrated liquid electrolyte
00093    * phase which obeys the Pitzer formulation for nonideality.
00094    *
00095    * As a prerequisite to the specification of thermodynamic quantities,
00096    * The concentrations of the ionic species are assumed to obey the 
00097    * electroneutrality condition. 
00098    *
00099    * <HR>
00100    * <H2> Specification of Species Standard %State Properties </H2>
00101    * <HR>
00102    *
00103    * The solvent is assumed to be liquid water. A real model for liquid
00104    * water (IAPWS 1995 formulation) is used as its standard state.
00105    * All standard state properties for the solvent are based on
00106    * this real model for water, and involve function calls
00107    * to the object that handles the real water model, #Cantera::WaterPropsIAPWS.
00108    *
00109    * The standard states for solutes are on the unit molality basis. 
00110    * Therefore, in the documentation below, the normal \f$ o \f$ 
00111    * superscript is replaced with
00112    * the \f$ \triangle \f$ symbol. The reference state symbol is now
00113    *  \f$ \triangle, ref \f$.
00114    * 
00115    *  
00116    *  It is assumed that the reference state thermodynamics may be
00117    *  obtained by a pointer to a populated species thermodynamic property
00118    *  manager class (see ThermoPhase::m_spthermo). How to relate pressure
00119    *  changes to the reference state thermodynamics is resolved at this level.
00120    *
00121    *  For solutes that rely on ThermoPhase::m_spthermo, are assumed to
00122    *  have an incompressible standard state mechanical property.
00123    *  In other words, the molar volumes are independent of temperature
00124    *  and pressure.
00125    *
00126    *  For these incompressible,
00127    *  standard states, the molar internal energy is
00128    *  independent of pressure. Since the thermodynamic properties
00129    *  are specified by giving the standard-state enthalpy, the
00130    *  term \f$ P_0 \hat v\f$ is subtracted from the specified molar
00131    *  enthalpy to compute the molar internal energy. The entropy is
00132    *  assumed to be independent of the pressure.
00133    *
00134    * The enthalpy function is given by the following relation.
00135    *
00136    *       \f[
00137    *   \raggedright  h^\triangle_k(T,P) = h^{\triangle,ref}_k(T) 
00138    *         + \tilde{v}_k \left( P - P_{ref} \right) 
00139    *       \f]
00140    *
00141    * For an incompressible,
00142    * stoichiometric substance, the molar internal energy is
00143    * independent of pressure. Since the thermodynamic properties
00144    * are specified by giving the standard-state enthalpy, the
00145    * term \f$ P_{ref} \tilde v\f$ is subtracted from the specified reference molar
00146    * enthalpy to compute the molar internal energy.
00147    *
00148    *       \f[
00149    *            u^\triangle_k(T,P) = h^{\triangle,ref}_k(T) - P_{ref} \tilde{v}_k
00150    *       \f]
00151    *
00152    *
00153    * The solute standard state heat capacity and entropy are independent
00154    * of pressure. The solute standard state gibbs free energy is obtained
00155    * from the enthalpy and entropy functions.
00156    *
00157    * The vector Constituents::m_speciesSize[] is used to hold the
00158    * base values of species sizes. These are defined as the 
00159    * molar volumes of species at infinite dilution at 300 K and 1 atm
00160    * of water. m_speciesSize are calculated during the initialization of the
00161    * %HMWSoln object and are then not touched.
00162    *
00163    * The current model assumes that an incompressible molar volume for
00164    * all solutes. The molar volume for the water solvent, however,
00165    * is obtained from a pure water equation of state, waterSS.
00166    * Therefore, the water standard state varies with both T and P.
00167    * It is an error to request standard state water properties  at a T and P
00168    * where the water phase is not a stable phase, i.e., beyond its
00169    * spinodal curve.
00170    *   
00171    * <HR>
00172    * <H2> Specification of Solution Thermodynamic Properties </H2>
00173    * <HR>
00174    *
00175    * Chemical potentials
00176    * of the solutes,  \f$ \mu_k \f$, and the solvent, \f$ \mu_o \f$, which are based 
00177    * on the molality form, have the following general format:
00178    *
00179    * \f[
00180    *    \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} \frac{m_k}{m^\triangle}) 
00181    * \f]
00182    * \f[
00183    *    \mu_o = \mu^o_o(T,P) + RT ln(a_o) 
00184    * \f]
00185    *
00186    * where \f$ \gamma_k^{\triangle} \f$ is the molality based activity coefficient for species
00187    * \f$k\f$.
00188    * 
00189    * Individual activity coefficients of ions can not be independently measured. Instead,
00190    * only binary pairs forming electroneutral solutions can be measured. This problem
00191    * leads to a redundancy in the evaluation of species standard state properties.
00192    * The redundancy issue is resolved by setting the standard state chemical potential
00193    * enthalpy, entropy, and volume for the hydrogen ion, H+, to zero, for every temperature
00194    * and pressure. After this convention is applied, all other standard state 
00195    * properties of ionic species contain meaningfull information.
00196    *
00197    *
00198    *  <H3> Ionic Strength </H3>
00199    *
00200    *  Most of the parameterizations within the model use the ionic strength 
00201    *  as a key variable. The ionic strength, \f$ I\f$ is defined as follows
00202    *
00203    *  \f[
00204    *    I = \frac{1}{2} \sum_k{m_k  z_k^2}
00205    *  \f]
00206    *
00207    *
00208    * \f$ m_k \f$ is the molality of the kth species. \f$ z_k \f$ is the charge
00209    * of the kth species. Note, the ionic strength is a defined units quantity.
00210    * The molality has defined units of gmol kg-1, and therefore the ionic
00211    * strength has units of sqrt( gmol kg<SUP>-1</SUP>).
00212    *
00213    * In some instances, from some authors, a different 
00214    * formulation is used for the ionic strength in the equations below. The different
00215    * formulation is due to the possibility of the existence of weak acids and how
00216    * association wrt to the weak acid equilibrium relation affects the calculation 
00217    * of the activity coefficients via the assumed value of the ionic strength.
00218    *
00219    * If we are to assume that the association reaction doesn't have an effect
00220    * on the ionic strength, then we will want to consider the associated weak
00221    * acid as in effect being fully dissociated, when we calculate an effective
00222    * value for the ionic strength. We will call this calculated value, the
00223    * stoichiometric ionic strength, \f$ I_s \f$, putting a subscript s to denote
00224    * it from the more straightforward calculation of \f$ I \f$.
00225    *
00226    *  \f[
00227    *    I_s = \frac{1}{2} \sum_k{m_k^s  z_k^2}
00228    *  \f]
00229    *
00230    *  Here, \f$ m_k^s \f$ is the value of the molalities calculated assuming that
00231    *  all weak acid-base pairs are in their fully dissociated states. This calculation may
00232    * be simplified by considering that the weakly associated acid may be made up of two
00233    * charged species, k1 and k2, each with their own charges, obeying the following relationship:
00234    *
00235    *   \f[
00236    *      z_k = z_{k1} +  z_{k2}
00237    *   \f]
00238    *  Then, we may only need to specify one charge value, say, \f$  z_{k1}\f$, 
00239    *  the cation charge number,
00240    *  in order to get both numbers, since we have already specified \f$ z_k \f$ 
00241    *  in the definition of original species.
00242    *  Then, the stoichiometric ionic strength may be calculated via the following formula.
00243    *
00244    *  \f[
00245    *    I_s = \frac{1}{2} \left(\sum_{k,ions}{m_k  z_k^2}+ 
00246    *               \sum_{k,weak_assoc}(m_k  z_{k1}^2 + m_k  z_{k2}^2) \right)
00247    *  \f]
00248    *
00249    *  The specification of which species are weakly associated acids is made in the input 
00250    *  file via the
00251    *  <TT> stoichIsMods </TT> XML block, where the charge for k1 is also specified. 
00252    *  An example is given below:
00253    * 
00254    * @code
00255    *          <stoichIsMods>
00256    *                NaCl(aq):-1.0
00257    *          </stoichIsMods>
00258    * @endcode
00259    *
00260    *
00261    *  Because we need the concept of a weakly associated acid in order to calculated 
00262    *  \f$ I_s \f$ we need to 
00263    *  catalog all species in the phase. This is done using the following categories:
00264    *
00265    *  -  <B>cEST_solvent</B>    :           Solvent species (neutral)
00266    *  -  <B>cEST_chargedSpecies</B>         Charged species (charged)
00267    *  -  <B>cEST_weakAcidAssociated</B>     Species which can break apart into charged species.
00268    *                                        It may or may not be charged.  These may or 
00269    *                                        may not be be included in the
00270    *                                        species solution vector.
00271    *  -  <B>cEST_strongAcidAssociated</B>   Species which always breaksapart into charged species.
00272    *                                        It may or may not be charged. Normally, these
00273    *                                        aren't included in the speciation vector.
00274    *  -  <B>cEST_polarNeutral </B>          Polar neutral species
00275    *  -  <B>cEST_nonpolarNeutral</B>        Non poloar neutral species
00276    *
00277    *  Polar and non-polar neutral species are differentiated, because some additions 
00278    *  to the activity 
00279    *  coefficient expressions distinguish between these two types of solutes. 
00280    *  This is the so-called salt-out effect.
00281    *
00282    * The type of species is specified in the <TT>electrolyteSpeciesType</TT> XML block.
00283    * Note, this is not
00284    * considered a part of the specification of the standard state for the species, 
00285    * at this time. Therefore,
00286    * this information is put under the <TT>activityCoefficient</TT> XML block. An example 
00287    * is given below
00288    *
00289    * @code
00290    *         <electrolyteSpeciesType>
00291    *                H2L(L):solvent
00292    *                H+:chargedSpecies
00293    *                NaOH(aq):weakAcidAssociated
00294    *                NaCl(aq):strongAcidAssociated
00295    *                NH3(aq):polarNeutral
00296    *                O2(aq):nonpolarNeutral
00297    *         </electrolyteSpeciesType>
00298    * @endcode
00299    *
00300    *
00301    *  Much of the species electrolyte type information is infered from other information in the
00302    *  input file. For example, as species which is charged is given the "chargedSpecies" default
00303    *  category. A neutral solute species is put into the "nonpolarNeutral" category by default.
00304    *
00305    *
00306    *  <H3> Specification of the Excess Gibbs Free Energy </H3>
00307    *
00308    *  Pitzer's formulation may best be represented as a specification of the excess gibbs
00309    *  free energy, \f$ G^{ex} \f$, defined as the deviation of the total gibbs free energy from 
00310    *  that of an ideal molal solution. 
00311    *     \f[
00312    *         G = G^{id} + G^{ex}
00313    *     \f]
00314    *
00315    *  The ideal molal solution contribution, not equal to an ideal solution contribution
00316    *  and in fact containing a singularity at the zero solvent mole fraction limit, is
00317    *  given below.
00318    *     \f[
00319    *         G^{id} = n_o \mu^o_o + \sum_{k\ne o} n_k \mu_k^{\triangle}
00320    *               + \tilde{M}_o n_o ( RT (\sum{m_i(\ln(m_i)-1)}))
00321    *     \f]
00322    *
00323    *  From the excess Gibbs free energy formulation, the activity coefficient expression
00324    *  and the osmotic coefficient expression for the solvent may be defined, by 
00325    *  taking the appropriate derivatives. Using this approach garranties that the
00326    *  entire system will obey the Gibbs-Duhem relations.
00327    *
00328    *  Pitzer employs the following general expression for the excess Gibbs free energy
00329    *
00330    *  \f[    
00331    *    \begin{array}{cclc}
00332    *     \frac{G^{ex}}{\tilde{M}_o  n_o RT} &= &
00333    *          \left( \frac{4A_{Debye}I}{3b} \right) \ln(1 + b \sqrt{I})
00334    *        +   2 \sum_c \sum_a m_c m_a B_{ca} 
00335    *        +     \sum_c \sum_a m_c m_a Z C_{ca}
00336    *          \\&&
00337    *        +   \sum_{c < c'} \sum m_c m_{c'} \left[ 2 \Phi_{c{c'}} + \sum_a m_a \Psi_{c{c'}a} \right]
00338    *        +   \sum_{a < a'} \sum m_a m_{a'} \left[ 2 \Phi_{a{a'}} + \sum_c m_c \Psi_{a{a'}c} \right]
00339    *          \\&&
00340    *        + 2 \sum_n \sum_c m_n m_c \lambda_{nc} + 2 \sum_n \sum_a m_n m_a \lambda_{na} 
00341    *        + 2 \sum_{n < n'} \sum m_n m_{n'} \lambda_{n{n'}}
00342    *        +  \sum_n m^2_n \lambda_{nn}         
00343    *    \end{array}
00344    *  \f]
00345    *
00346    *  <I>a</I> is a subscribt over all anions, <I>c</I> is a subscript extending over all
00347    *  cations, and  <I>i</I> is a subscrit that extends over all anions and cations.
00348    *  <I>n</I> is a subscript that extends only over neutral solute molecules. 
00349    *  The second line contains cross terms where cations affect 
00350    *  cations and/or cation/anion pairs,
00351    *  and anions affect anions or cation/anion pairs. Note part of the coefficients,
00352    *  \f$ \Phi_{c{c'}} \f$ and  \f$ \Phi_{a{a'}} \f$  stem from the theory
00353    *  of unsymmetrical mixing of electrolytes with different charges. This
00354    *  theory depends on the total ionic stregnth of the solution, and therefore,
00355    *  \f$ \Phi_{c{c'}} \f$ and  \f$ \Phi_{a{a'}} \f$  will depend on <I>I</I>, the
00356    *  ionic strength.  \f$ B_{ca}\f$ is a strong function of the 
00357    *  total ionic strength,  <I>I</I>,
00358    *  of the electrolyte. The rest of the coefficients are assumed to be independent of the
00359    *  molalities or ionic strengths. However, all coefficients are potentially functions
00360    *  of the temperature and pressure of the solution.
00361    *  
00362    *   <I>A</I> is the Debye-Huckel constant. It's specification is described in its own
00363    *   section below. 
00364    *
00365    *   \f$ I\f$ is the ionic strength of the solution, and is given by:
00366    *
00367    *   \f[
00368    *       I = \frac{1}{2} \sum_k{m_k  z_k^2}
00369    *   \f]
00370    *
00371    *   In contrast to several other Debye-Huckel implementations (see \ref DebyeHuckel), the 
00372    *   parameter \f$ b\f$ in the above equation is a constant that 
00373    *   doesn not vary with respect to ion idenity. This is an important simplification
00374    *   as it avoids troubles with satisfaction of the Gibbs-Duhem analysis.  
00375    *
00376    *   The function \f$ Z \f$ is given by
00377    *
00378    *   \f[
00379    *       Z = \sum_i m_i \left| z_i \right|
00380    *   \f]
00381    *
00382    *   The value of \f$ B_{ca}\f$ is given by the following function
00383    *
00384    *   \f[
00385    *      B_{ca} = \beta^{(0)}_{ca} + \beta^{(1)}_{ca} g(\alpha^{(1)}_{ca} \sqrt{I}) 
00386    *             + \beta^{(2)}_{ca} g(\alpha^{(2)}_{ca} \sqrt{I}) 
00387    *   \f]
00388    *
00389    *   where
00390    *
00391    *   \f[
00392    *       g(x) = 2 \frac{(1 - (1 + x)\exp[-x])}{x^2}
00393    *   \f]
00394    *
00395    *   The formulation for \f$ B_{ca}\f$ combined with the formulation of the
00396    *   Debye-Huckel term in the eqn. for the excess Gibbs free energy stems
00397    *   essentially from an empirical fit to the ionic strength dependent data
00398    *   based over a wide sampling of binary electroyte systems. \f$ C_{ca} \f$,
00399    *   \f$ \lambda_{nc} \f$, \f$ \lambda_{na} \f$, \f$ \lambda_{nn} \f$,
00400    *   \f$ \Psi_{c{c'}a} \f$, \f$ \Psi_{a{a'}c} \f$ are experimentally derived
00401    *   coefficients that may have pressure and/or temperature dependencies.
00402    *    The \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$ formulations are 
00403    *   slightly more complicated. \f$ b \f$ is a univeral
00404    *   constant defined to be equal to \f$ 1.2\ kg^{1/2}\ gmol^{-1/2} \f$. The exponential
00405    *   coefficient \f$ \alpha^{(1)}_{ca} \f$ is usually 
00406    *   fixed at \f$ \alpha^{(1)}_{ca} = 2.0\ kg^{1/2} gmol^{-1/2}\f$
00407    *   except for 2-2 electrolytes, while other parameters were fit to experimental
00408    *   data. For 2-2 electrolytes, \f$ \alpha^{(1)}_{ca} = 1.4\ kg^{1/2}\ gmol^{-1/2}\f$
00409    *    is used in combination with either \f$ \alpha^{(2)}_{ca} = 12\ kg^{1/2}\ gmol^{-1/2}\f$
00410    *   or \f$ \alpha^{(2)}_{ca} = k A_\psi \f$, where <I>k</I> is a constant. For electrolytes other
00411    *   than 2-2 electrolytes the \f$ \beta^{(2)}_{ca} g(\alpha^{(2)}_{ca} \sqrt{I}) \f$  term
00412    *   is not used in the fitting procedure; it is only used for divalent metal
00413    *   solfates and other high-valence electrolytes which exhibit significant 
00414    *   association at low ionic strengths.
00415    *
00416    *   The \f$ \beta^{(0)}_{ca} \f$,  \f$ \beta^{(1)}_{ca}\f$,  \f$ \beta^{(2)}_{ca} \f$,
00417    *   and \f$ C_{ca} \f$ binary coefficients are referred to as ion-interaction or
00418    *   Pitzer parameters. These Pitzer parameters may vary with temperature and pressure
00419    *   but they do not depend on the ionic strength. Their values and temperature
00420    *   derivatives of their values have been tabulated for a range of electrolytes
00421    *   
00422    *   The \f$ \Phi_{c{c'}} \f$ and \f$ \Phi_{a{a'}} \f$ contributions, which
00423    *   capture cation-cation and anion-anion interactions, also have an
00424    *   ionic strength dependence.
00425    *
00426    *   Ternary contributions \f$ \Psi_{c{c'}a} \f$ and \f$ \Psi_{a{a'}c} \f$ 
00427    *   have been measured also for some systems. The success of the Pitzer
00428    *   method lies in its ability to model nonlinear activity coefficients
00429    *   of complex multicomponent systems with just binary and minor
00430    *   ternary contributions, which can be independently measured in
00431    *   binary or ternary subsystems.
00432    *
00433    *
00434    *  <H3> Multicomponent Activity Coefficients for Solutes </H3>
00435    *
00436    *    The formulas for activity coefficients of solutes may be obtained by taking the 
00437    *    following derivative of the excess Gibbs Free Energy formulation described above:
00438    *
00439    *    \f[
00440    *     \ln(\gamma_k^\triangle) = \frac{d\left( \frac{G^{ex}}{M_o n_o RT} \right)}{d(m_k)}\Bigg|_{n_i}
00441    *    \f]
00442    *
00443    *    In the formulas below the following conventions are used. The subscript <I>M</I> refers
00444    *    to a particular cation. The subscript X refers to a particular anion, whose
00445    *    activity is being currently evaluated. the subscript <I>a</I> refers to a summation
00446    *    over all anions in the solution, while the subscript <I>c</I> refers to a summation
00447    *    over all cations in the solutions.
00448    *
00449    *     The activity coefficient for a particular cation <I>M</I> is given by
00450    *
00451    *   \f[
00452    *      \ln(\gamma_M^\triangle) = -z_M^2(F) + \sum_a m_a \left( 2 B_{Ma} + Z C_{Ma} \right)
00453    *      + z_M   \left( \sum_a  \sum_c m_a m_c C_{ca} \right)
00454    *             + \sum_c m_c \left[ 2 \Phi_{Mc} + \sum_a m_a \Psi_{Mca} \right]
00455    *             + \sum_{a < a'} \sum m_a m_{a'} \Psi_{Ma{a'}}
00456    *             +  2 \sum_n m_n \lambda_{nM}
00457    *   \f]
00458    *
00459    *     The activity coefficient for a particular anion <I>X</I> is given by
00460    *
00461    *   \f[
00462    *      \ln(\gamma_X^\triangle) = -z_X^2(F) + \sum_a m_c \left( 2 B_{cX} + Z C_{cX} \right)
00463    *      + \left|z_X \right|  \left( \sum_a  \sum_c m_a m_c C_{ca} \right)
00464    *             + \sum_a m_a \left[ 2 \Phi_{Xa} + \sum_c m_c \Psi_{cXa} \right]
00465    *             + \sum_{c < c'} \sum m_c m_{c'} \Psi_{c{c'}X}
00466    *             +  2 \sum_n m_n \lambda_{nM}
00467    *   \f]
00468    *              where the function \f$ F \f$ is given by
00469    *
00470    *   \f[
00471    *       F = - A_{\phi} \left[ \frac{\sqrt{I}}{1 + b \sqrt{I}} 
00472    *                 + \frac{2}{b} \ln{\left(1 + b\sqrt{I}\right)} \right]
00473    *                 + \sum_a \sum_c m_a m_c B'_{ca}
00474    *                 + \sum_{c < c'} \sum m_c m_{c'} \Phi'_{c{c'}}
00475    *                 + \sum_{a < a'} \sum m_a m_{a'} \Phi'_{a{a'}}
00476    *   \f]
00477    *
00478    *   We have employed the definition of \f$ A_{\phi} \f$, also used by Pitzer
00479    *   which is equal to
00480    *
00481    *   \f[
00482    *     A_{\phi} = \frac{A_{Debye}}{3}
00483    *   \f]
00484    *
00485    *   In the above formulas, \f$ \Phi'_{c{c'}} \f$  and \f$ \Phi'_{a{a'}} \f$ are the
00486    *   ionic strength derivatives of \f$ \Phi_{c{c'}} \f$  and \f$  \Phi_{a{a'}} \f$,
00487    *   respectively.
00488    *
00489    *   The function \f$ B'_{MX} \f$ is defined as:
00490    *
00491    *   \f[
00492    *       B'_{MX} = \left( \frac{\beta^{(1)}_{MX} h(\alpha^{(1)}_{MX} \sqrt{I})}{I}  \right) 
00493    *                 \left( \frac{\beta^{(2)}_{MX} h(\alpha^{(2)}_{MX} \sqrt{I})}{I}  \right)
00494    *   \f]
00495    *
00496    *  where \f$ h(x) \f$ is defined as
00497    *
00498    *   \f[
00499    *       h(x) = g'(x) \frac{x}{2} = 
00500    *        \frac{2\left(1 - \left(1 + x + \frac{x^2}{2} \right)\exp(-x) \right)}{x^2}
00501    *   \f]
00502    *
00503    *   The activity coefficient for neutral species <I>N</I> is given by 
00504    *
00505    *   \f[
00506    *       \ln(\gamma_N^\triangle) = 2 \left( \sum_i m_i \lambda_{iN}\right)
00507    *   \f]
00508    *
00509    *
00510    *  <H3> Activity of the Water Solvent </H3> 
00511    *
00512    *  The activity for the solvent water,\f$ a_o \f$, is not independent and must be 
00513    *  determined either from the Gibbs-Duhem relation or from taking the appropriate derivative
00514    *  of the same excess Gibbs free energy function as was used to formulate
00515    *  the solvent activity coefficients. Pitzer's description follows the later approach to
00516    *  derive a formula for the osmotic coefficient, \f$ \phi \f$.
00517    *
00518    *  \f[
00519    *       \phi - 1 = - \left( \frac{d\left(\frac{G^{ex}}{RT} \right)}{d(\tilde{M}_o n_o)}  \right)
00520    *                \frac{1}{\sum_{i \ne 0} m_i}
00521    *  \f]
00522    *
00523    *  The osmotic coefficient may be related to the water activity by the following relation:
00524    * 
00525    *  \f[
00526    *       \phi = - \frac{1}{\tilde{M}_o \sum_{i \neq o} m_i} \ln(a_o) 
00527    *           = - \frac{n_o}{\sum_{i \neq o}n_i} \ln(a_o)
00528    *  \f]
00529    *
00530    *
00531    *  The result is the following
00532    *
00533    *  \f[    
00534    *    \begin{array}{ccclc}
00535    *      \phi - 1 &= &
00536    *          \frac{2}{\sum_{i \ne 0} m_i}
00537    *           \bigg[ &
00538    *        -  A_{\phi} \frac{I^{3/2}}{1 + b \sqrt{I}}  
00539    *        +   \sum_c  \sum_a m_c m_a \left( B^{\phi}_{ca} + Z C_{ca}\right)
00540    *          \\&&&
00541    *        +   \sum_{c < c'} \sum m_c m_{c'} \left[ \Phi^{\phi}_{c{c'}} + \sum_a m_a \Psi_{c{c'}a} \right]
00542    *        +   \sum_{a < a'} \sum m_a m_{a'} \left[ \Phi^{\phi}_{a{a'}} + \sum_c m_c \Psi_{a{a'}c} \right]
00543    *          \\&&&
00544    *        + \sum_n \sum_c m_n m_c \lambda_{nc} +  \sum_n \sum_a m_n m_a \lambda_{na} 
00545    *        + \sum_{n < n'} \sum m_n m_{n'} \lambda_{n{n'}}
00546    *        + \frac{1}{2} \left( \sum_n m^2_n \lambda_{nn}\right)         
00547    *          \bigg]
00548    *    \end{array}
00549    *  \f]
00550    *
00551    *  It can be shown that the expression 
00552    *
00553    *
00554    *  \f[
00555    *     B^{\phi}_{ca} = \beta^{(0)}_{ca} + \beta^{(1)}_{ca} \exp{(- \alpha^{(1)}_{ca} \sqrt{I})} 
00556    *             + \beta^{(2)}_{ca} \exp{(- \alpha^{(2)}_{ca} \sqrt{I} )} 
00557    *  \f]
00558    *     
00559    *  is consistent with the expression \f$ B_{ca} \f$ in the \f$ G^{ex} \f$ expression
00560    *  after carrying out the derivative wrt \f$ m_M \f$.
00561    *
00562    *  Also taking into account that  \f$ {\Phi}_{c{c'}} \f$ and
00563    *  \f$ {\Phi}_{a{a'}} \f$ has an ionic strength dependence.
00564    *
00565    *  \f[
00566    *    \Phi^{\phi}_{c{c'}} = {\Phi}_{c{c'}} + I \frac{d{\Phi}_{c{c'}}}{dI}
00567    *  \f]
00568    *
00569    *  \f[
00570    *    \Phi^{\phi}_{a{a'}} = \Phi_{a{a'}} + I \frac{d\Phi_{a{a'}}}{dI}
00571    *  \f]
00572    *
00573    *
00574    *  <H3> Temperature and Pressure Dependence of the Pitzer Parameters </H3>
00575    *
00576    *  In general most of the coefficients introduced in the previous section may
00577    *  have a temperature and pressure dependence. The temperature and pressure
00578    *  dependence of these coefficients strongly influence the value of the
00579    *  excess Enthalpy and excess Volumes of Pitzer solutions. Therefore, these
00580    *  are readily measurable quantities.
00581    *  HMWSoln provides several
00582    *  different methods for putting these dependencies into the coefficients.
00583    *  HMWSoln has an implementation described by Silverter and Pitzer (1977),
00584    *  which was used to fit experimental data for NaCl over an extensive range,
00585    *  below the critical temperature of water.
00586    *  They found a temperature funcdtional form for fitting the 3 following
00587    *  coefficients that describe the Pitzer parameterization for a single salt
00588    *  to be adequate to describe how the excess gibbs free energy values for
00589    *  the binary salt changes with respect to temperature.
00590    *  The following functional form
00591    *  was used to fit the temperature dependence of the Pitzer Coefficients
00592    *  for each cation - anion pair, M X.
00593    *
00594    *  \f[
00595    *      \beta^{(0)}_{MX} = q^{b0}_0 
00596    *                       + q^{b0}_1 \left( T - T_r \right)
00597    *                       + q^{b0}_2 \left( T^2 - T_r^2 \right)
00598    *                       + q^{b0}_3 \left( \frac{1}{T} - \frac{1}{T_r}\right)
00599    *                       + q^{b0}_4 \ln \left( \frac{T}{T_r} \right)
00600    *  \f]
00601    *  \f[
00602    *      \beta^{(1)}_{MX} = q^{b1}_0  + q^{b1}_1 \left( T - T_r \right)
00603    *                       + q^{b1}_{2} \left( T^2 - T_r^2 \right)
00604    *  \f]
00605    *  \f[
00606    *     C^{\phi}_{MX} = q^{Cphi}_0 
00607    *                   + q^{Cphi}_1 \left( T - T_r \right)
00608    *                   + q^{Cphi}_2 \left( T^2 - T_r^2 \right)
00609    *                   + q^{Cphi}_3 \left( \frac{1}{T} - \frac{1}{T_r}\right)
00610    *                   + q^{Cphi}_4 \ln \left( \frac{T}{T_r} \right)
00611    *  \f]
00612    *
00613    *  where
00614    *
00615    *  \f[ 
00616    *      C^{\phi}_{MX} =  2 {\left| z_M z_X \right|}^{1/2} C_{MX}
00617    *  \f]
00618    *
00619    *  In later papers, Pitzer has added additional temperature dependencies
00620    *  to all of the other remaining second and third order virial coefficients.
00621    *  Some of these dependencies are justified and motivated by theory. 
00622    *  Therefore,
00623    *  a formalism wherein all of the coefficients in the base theory have
00624    *  temperature dependencies associated with them has been implemented 
00625    *  within the  %HMWSoln object. Much of the formalism, however,
00626    *  has been unexercised.
00627    *
00628    *  In the %HMWSoln object, the temperature dependence of the Pitzer 
00629    *  parameters are specified in the following way.
00630    *
00631    *    - PIZTER_TEMP_CONSTANT      - string name "CONSTANT"
00632    *      - Assumes that all coefficients are independent of temperature
00633    *        and pressure
00634    *    - PIZTER_TEMP_COMPLEX1     - string name "COMPLEX" or "COMPLEX1"
00635    *      - Uses the full temperature dependence for the
00636    *        \f$\beta^{(0)}_{MX} \f$ (5 coeffs), 
00637    *        the  \f$\beta^{(1)}_{MX} \f$ (3 coeffs),
00638    *        and \f$ C^{\phi}_{MX} \f$ (5 coeffs) parameters described above.
00639    *    - PITZER_TEMP_LINEAR        - string name "LINEAR"
00640    *      - Uses just the temperature dependence for the
00641    *        \f$\beta^{(0)}_{MX} \f$, the \f$\beta^{(1)}_{MX} \f$,
00642    *        and \f$ C^{\phi}_{MX} \f$ coefficients described above. 
00643    *        There are 2 coefficients for each term.
00644    *
00645    *  The temperature dependence is specified in an attributes field
00646    *  in the <TT>  activityCoefficients </TT>  XML block, 
00647    *  called <TT> TempModel </TT>. Permissible values for that 
00648    *  attribute are <TT> CONSTANT, COMPLEX1</TT>, and <TT> LINEAR.</TT>
00649    *
00650    *  The specification of the binary interaction between a cation and
00651    *  an anion is given by the coefficients, \f$ B_{MX}\f$ and 
00652    *  \f$ C_{MX}\f$
00653    *  The specification of \f$ B_{MX}\f$ is a function of 
00654    *  \f$\beta^{(0)}_{MX} \f$, \f$\beta^{(1)}_{MX} \f$, 
00655    *  \f$\beta^{(2)}_{MX} \f$,  \f$\alpha^{(1)}_{MX} \f$, and
00656    *  \f$\alpha^{(2)}_{MX} \f$.
00657    *  \f$ C_{MX}\f$ is calculated from \f$C^{\phi}_{MX} \f$
00658    *  from the formula above.
00659    *  All of the underlying coeficients are specified in the
00660    *  XML element block <TT> binarySaltParameters </TT>, which
00661    *  has the attribute <TT> cation </TT> and <TT> anion </TT>
00662    *  to identify the interaction. XML elements named
00663    *  <TT> beta0, beta1, beta2, Cphi, Alpha1, Alpha2 </TT> 
00664    *  within each <TT> binarySaltParameters </TT> block
00665    *  specify the parameters. Within each of these blocks 
00666    *  multiple parameters describing temperature or pressure
00667    *  dependence are serially listed in the order that they
00668    *  appear in the equation in this document. An example of
00669    *  the <TT> beta0 </TT> block that fits the <TT> COMPLEX1 </TT> temperature
00670    *  dependence given above is
00671    *
00672    * @code
00673       <binarySaltParameters cation="Na+" anion="OH-"> 
00674         <beta0> q0, q1, q2, q3, q4  </beta0>
00675       </binarySaltParameters>
00676      @endcode
00677    *
00678    *  The parameters for \f$ \beta^{(0)}\f$ fit the following equation:
00679    *
00680    *  \f[
00681    *        \beta^{(0)} = q_0^{{\beta}0} + q_1^{{\beta}0} \left( T - T_r \right) 
00682    *               + q_2^{{\beta}0} \left( T^2 - T_r^2 \right)  
00683    *               + q_3^{{\beta}0} \left( \frac{1}{T} - \frac{1}{T_r} \right)
00684    *               + q_4^{{\beta}0} \ln \left( \frac{T}{T_r} \right)
00685    *  \f]
00686    *
00687    *  This same COMPLEX1 </TT> temperature
00688    *  dependence given above is used for the following parameters:
00689    *  \f$ \beta^{(0)}_{MX} \f$, \f$ \beta^{(1)}_{MX} \f$, 
00690    *  \f$ \beta^{(2)}_{MX} \f$, \f$ \Theta_{cc'} \f$,  \f$\Theta_{aa'} \f$,
00691    *  \f$ \Psi_{c{c'}a} \f$ and \f$ \Psi_{ca{a'}} \f$.
00692    *
00693    *
00694    *  <H3> Like-Charged Binary Ion Parameters and the Mixing Parameters </H3>
00695    *
00696    *  The previous section contained the functions, \f$ \Phi_{c{c'}} \f$,
00697    *  \f$ \Phi_{a{a'}} \f$ and their derivatives wrt the
00698    *  ionic strength, \f$ \Phi'_{c{c'}} \f$ and \f$ \Phi'_{a{a'}} \f$.
00699    *  Part of these terms come from theory.
00700    *
00701    *  Since like charged ions repel each other and are generally
00702    *  not near each other, the virial coefficients for same-charged ions
00703    *  are small. However, Pitzer doesn't ignore these in his 
00704    *  formulation. Relatively larger and longer range terms between
00705    *  like-charged ions exist however, which appear only for 
00706    *  unsymmetrical mixing of same-sign charged ions with different
00707    *  charges. \f$ \Phi_{ij} \f$, where \f$ ij \f$ is either \f$ a{a'} \f$
00708    *  or \f$ c{c'} \f$ is given by
00709    *
00710    *
00711    *
00712    *  \f[
00713    *      {\Phi}_{ij} = \Theta_{ij} + \,^E \Theta_{ij}(I)
00714    *  \f]
00715    *
00716    *
00717    *  \f$ \Theta_{ij} \f$ is the small virial coefficient expansion term.
00718    *  Dependent in general on temperature and pressure, it's ionic
00719    *  strength dependence is ignored in Pitzer's approach. 
00720    *  \f$ \,^E\Theta_{ij}(I) \f$ accounts for the electrostatic
00721    *  unsymmetrical mixing effects and is dependent only on the
00722    *  charges of the ions i, j, the total ionic strength and on
00723    *  the dielectric constant and density of the solvent.
00724    *  This seems to be a relatively well-documented part of the theory.
00725    *  They theory below comes from Pitzer summation (Pitzer) in the
00726    *  appendix. It's also mentioned in Bethke's book (Bethke), and
00727    *  the equations are summarized in Harvie & Weare (1980).
00728    *  Within the code, \f$ \,^E\Theta_{ij}(I) \f$ is evaluated according
00729    *  to the algorithm described in Appendix B [Pitzer] as
00730    *
00731    *  \f[ 
00732    *     \,^E\Theta_{ij}(I) = \left( \frac{z_i z_j}{4I} \right)
00733    *        \left( J(x_{ij})  - \frac{1}{2} J(x_{ii})  
00734    *                          - \frac{1}{2} J(x_{jj})  \right)
00735    *  \f]
00736    *
00737    *  where \f$ x_{ij} = 6 z_i z_j A_{\phi} \sqrt{I} \f$ and
00738    *
00739    *   \f[ 
00740    *      J(x) = \frac{1}{x} \int_0^{\infty}{\left( 1 + q +
00741    *             \frac{1}{2} q^2 - e^q \right) y^2 dy}
00742    *   \f]
00743    *
00744    *  and \f$ q = - (\frac{x}{y}) e^{-y} \f$. \f$ J(x) \f$ is evaluated by 
00745    *  numerical integration.
00746    *
00747    *  The \f$  \Theta_{ij} \f$ term is a constant that is specified
00748    *  by the XML element <TT> thetaCation </TT> and 
00749    *  <TT> thetaAnion </TT>,  which
00750    *  has the attribute <TT> cation1 </TT>, <TT> cation2 </TT> and
00751    *  <TT> anion1 </TT>, <TT> anion2 </TT> respectively
00752    *  to identify the interaction. No temperature or 
00753    *  pressure dependence of this parameter is currently allowed.
00754    *  An example of the block is presented below.
00755    *
00756    * @code
00757        <thetaCation cation1="Na+" cation2="H+">
00758                   <Theta> 0.036 </Theta>
00759        </thetaCation>
00760      @endcode
00761    *
00762    *
00763    *  <H3> Ternary Pitzer Parameters </H3>
00764    *
00765    *  The \f$  \Psi_{c{c'}a} \f$ and \f$  \Psi_{ca{a'}} \f$ terms
00766    *  represent ternary interactions between two cations and
00767    *  an anion and two anions and a cation, respectively.
00768    *  In Pitzer's implementation these terms are usually small
00769    *  in absolute size. Currently these parameters do not have
00770    *  any dependence on temperature, pressure, or ionic strength.
00771    *  
00772    *   Their values are input using the XML element
00773    *   <TT> psiCommonCation </TT> and  <TT> psiCommonAnion </TT>.
00774    *   The species id's are specified in attribute fields in
00775    *   the XML element. The fields  <TT>cation</TT>, 
00776    *   <TT> anion1</TT>, and <TT> anion2</TT>
00777    *   are used for  <TT>psiCommonCation</TT>. The fields <TT> anion</TT>,
00778    *   <TT>cation1</TT> and  <TT>cation2</TT> are used for 
00779    *   <TT> psiCommonAnion</TT>. An example block is given below.
00780    *   The <TT> Theta </TT> field below is a duplicate of the
00781    *   <TT> thetaAnion </TT> field mentioned above. The two fields
00782    *   are input into the same block for convenience, and because
00783    *   their data are highly correlated, in practice.
00784    *   It is an error for the
00785    *   two blocks to specify different information about 
00786    *   thetaAnion (or thetaCation) in different blocks. It's
00787    *   ok to specify duplicate but consistent information
00788    *   in multiple blocks.
00789    *
00790    * @code
00791        <psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
00792              <Theta> -0.05 </Theta>
00793              <Psi> -0.006 </Psi>
00794        </psiCommonCation>
00795       @endcode
00796    *
00797    *  <H3> Treatment of Neutral Species </H3>
00798    *
00799    *  Binary virial-coefficient-like interactions between two neutral 
00800    *  species may be specified in the \f$ \lambda_{mn} \f$ terms
00801    *  that appear in the formulas above.
00802    *  Currently these interactions are independent of temperature,
00803    *  pressure, and ionic strength. Also, currently, the neutrality
00804    *  of the species are not checked. Therefore, this interaction
00805    *  may involve charged species in the solution as well.
00806    *  The identity of the species is specified by the 
00807    *  <TT>species1</TT> and <TT>species2</TT> attributes to the XML
00808    *  <TT>lambdaNeutral</TT> node. These terms are symmetrical;
00809    *  <TT>species1</TT> and <TT>species2</TT> may be reversed and
00810    *  the term will be the same. An example is given below.
00811    *
00812    * @code
00813        <lambdaNeutral species1="CO2" species2="CH4">
00814              <lambda> 0.05 </lambda>
00815        </lambdaNeutral>
00816      @endcode
00817    *
00818    *  <H3> Example of the Specification of Parameters for the Activity
00819    *       Coefficients </H3>
00820    *  
00821    * An example is given below.
00822    *
00823    * An example <TT> activityCoefficients </TT> XML block for this 
00824    * formulation is supplied below
00825    *
00826    * @verbatim
00827     <activityCoefficients model="Pitzer" TempModel="complex1">
00828                 <!-- Pitzer Coefficients
00829                      These coefficients are from Pitzer's main
00830                      paper, in his book.
00831                   -->
00832                 <A_Debye model="water" />
00833                 <ionicRadius default="3.042843"  units="Angstroms">
00834                 </ionicRadius>
00835                 <binarySaltParameters cation="Na+" anion="Cl-">
00836                   <beta0> 0.0765, 0.008946, -3.3158E-6,
00837                           -777.03, -4.4706
00838                   </beta0>
00839                   <beta1> 0.2664, 6.1608E-5, 1.0715E-6, 0.0, 0.0 </beta1>
00840                   <beta2> 0.0,   0.0, 0.0, 0.0, 0.0  </beta2>
00841                   <Cphi> 0.00127, -4.655E-5, 0.0,
00842                          33.317, 0.09421
00843                   </Cphi>
00844                   <Alpha1> 2.0 </Alpha1>
00845                 </binarySaltParameters>
00846 
00847                 <binarySaltParameters cation="H+" anion="Cl-">
00848                   <beta0> 0.1775, 0.0, 0.0, 0.0, 0.0 </beta0>
00849                   <beta1> 0.2945, 0.0, 0.0, 0.0, 0.0 </beta1>
00850                   <beta2> 0.0,    0.0, 0.0, 0.0, 0.0 </beta2>
00851                   <Cphi> 0.0008, 0.0, 0.0, 0.0, 0.0 </Cphi>
00852                   <Alpha1> 2.0 </Alpha1>
00853                 </binarySaltParameters>
00854 
00855                 <binarySaltParameters cation="Na+" anion="OH-">
00856                   <beta0> 0.0864, 0.0, 0.0, 0.0, 0.0 </beta0>
00857                   <beta1> 0.253,  0.0, 0.0  0.0, 0.0 </beta1>
00858                   <beta2> 0.0     0.0, 0.0, 0.0, 0.0 </beta2>
00859                   <Cphi> 0.0044,  0.0, 0.0, 0.0, 0.0 </Cphi>
00860                   <Alpha1> 2.0 </Alpha1>
00861                 </binarySaltParameters>
00862 
00863                 <thetaAnion anion1="Cl-" anion2="OH-">
00864                   <Theta> -0.05,  0.0, 0.0, 0.0, 0.0 </Theta>
00865                 </thetaAnion>
00866 
00867                 <psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
00868                   <Theta> -0.05,  0.0, 0.0, 0.0, 0.0 </Theta>
00869                   <Psi> -0.006 </Psi>
00870                 </psiCommonCation>
00871 
00872                 <thetaCation cation1="Na+" cation2="H+">
00873                   <Theta> 0.036,  0.0, 0.0, 0.0, 0.0 </Theta>
00874                 </thetaCation>
00875 
00876                 <psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
00877                   <Theta> 0.036,  0.0, 0.0, 0.0, 0.0 </Theta>
00878                   <Psi> -0.004 </Psi>
00879                 </psiCommonAnion>
00880 
00881        </activityCoefficients>
00882     @endverbatim
00883    *
00884    *
00885    * <H3> Specification of the Debye-Huckel Constant </H3>
00886    *
00887    *  In the equations above, the formula for  \f$  A_{Debye} \f$ 
00888    *  is needed. The %HMWSoln object uses two methods for specifying these quantities.
00889    *  The default method is to assume that \f$  A_{Debye} \f$  is a constant, given
00890    *  in the initialization process, and storred in the
00891    *  member double, m_A_Debye. Optionally, a full water treatment may be employed that makes
00892    *  \f$ A_{Debye} \f$ a full function of <I>T</I> and <I>P</I> and creates nontrivial entries for
00893    *  the excess heat capacity, enthalpy, and excess volumes of solution.
00894    *
00895    *   \f[
00896    *      A_{Debye} = \frac{F e B_{Debye}}{8 \pi \epsilon R T} {\left( C_o \tilde{M}_o \right)}^{1/2}
00897    *   \f]
00898    * where
00899    * 
00900    *  \f[
00901    *      B_{Debye} = \frac{F} {{(\frac{\epsilon R T}{2})}^{1/2}} 
00902    *  \f]
00903    *  Therefore:
00904    *  \f[
00905    *      A_{Debye} = \frac{1}{8 \pi} 
00906    *                  {\left(\frac{2 N_a \rho_o}{1000}\right)}^{1/2}
00907    *                  {\left(\frac{N_a e^2}{\epsilon R T }\right)}^{3/2}
00908    *  \f]
00909    *
00910    *            Units = sqrt(kg/gmol)
00911    *
00912    *     where
00913    *      - \f$ N_a \f$ is Avrogadro's number
00914    *      - \f$ \rho_w \f$ is the density of water
00915    *      - \f$ e \f$ is the electronic charge
00916    *      - \f$ \epsilon = K \epsilon_o \f$ is the permitivity of water
00917    *           where \f$ K \f$ is the dielectric constant of water,
00918    *           and  \f$ \epsilon_o \f$ is the permitivity of free space.
00919    *      - \f$ \rho_o \f$ is the density of the solvent in its standard state.
00920    *
00921    *            Nominal value at 298 K and 1 atm = 1.172576 (kg/gmol)<SUP>1/2</SUP>
00922    *                  based on:
00923    *                 -   \f$ \epsilon / \epsilon_0 \f$ = 78.54
00924    *                           (water at 25C)
00925    *                 -   \f$ \epsilon_0 \f$= 8.854187817E-12 C<SUP>2</SUP> N<SUP>-1</SUP> m<SUP>-2</SUP>
00926    *                 -   e = 1.60217653E-19 C
00927    *                 -   F = 9.6485309E7 C kmol<SUP>-1</SUP>
00928    *                 -   R = 8.314472E3 kg m<SUP>2</SUP> s<SUP>-2</SUP> kmol<SUP>-1</SUP> K<SUP>-1</SUP>
00929    *                 -   T = 298.15 K
00930    *                 -   B_Debye = 3.28640E9 (kg/gmol)<SUP>1/2</SUP> m<SUP>-1</SUP>
00931    *                 -   \f$N_a\f$ = 6.0221415E26 kmol<SUP>-1</SUP>
00932    *
00933    * An example of a fixed value implementation is given below.
00934    * @code
00935    *   <activityCoefficients model="Pitzer">
00936    *         <!-- A_Debye units = sqrt(kg/gmol)  -->
00937    *         <A_Debye> 1.172576 </A_Debye>
00938    *         <!-- object description continues -->
00939    *   </activityCoefficients>
00940    * @endcode
00941    *
00942    * An example of a variable value implementation within the %HMWSoln object is given below.
00943    * The model attribute, "water", triggers the full implementation.
00944    *
00945    * @code
00946    *   <activityCoefficients model="Pitzer">
00947    *         <!-- A_Debye units = sqrt(kg/gmol)  -->
00948    *         <A_Debye model="water" /> 
00949    *         <!-- object description continues -->
00950    *   </activityCoefficients>
00951    * @endcode
00952    *
00953    *
00954    * <H3>  Temperature and Pressure Dependence of the Activity Coefficients </H3>
00955    *
00956    *  Temperature dependence of the activity coefficients leads to nonzero terms
00957    *  for the excess enthalpy and entropy of solution. This means that the
00958    *  partial molar enthalpies, entropies, and heat capacities are all
00959    *  non-trivial to compute. The following formulas are used.
00960    *
00961    *  The partial molar enthalpy, \f$ \bar s_k(T,P) \f$:
00962    *  
00963    *  \f[
00964    * \bar h_k(T,P) = h^{\triangle}_k(T,P) 
00965    *               - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
00966    * \f]
00967    * The solvent partial molar enthalpy is equal to 
00968    *  \f[
00969    * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o)}{dT}
00970    *       = h^{o}_o(T,P) 
00971    *       + R T^2 (\sum_{k \neq o} m_k)  \tilde{M_o} (\frac{d \phi}{dT})
00972    * \f]
00973    *
00974    *  The partial molar entropy, \f$ \bar s_k(T,P) \f$:
00975    *
00976    *  \f[
00977    *     \bar s_k(T,P) =  s^{\triangle}_k(T,P) 
00978    *             - R \ln( \gamma^{\triangle}_k \frac{m_k}{m^{\triangle}}))
00979    *                    - R T \frac{d \ln(\gamma^{\triangle}_k) }{dT}
00980    * \f]
00981    * \f[
00982    *      \bar s_o(T,P) = s^o_o(T,P) - R \ln(a_o)
00983    *                    - R T \frac{d \ln(a_o)}{dT}
00984    * \f]
00985    *
00986    * The partial molar heat capacity, \f$ C_{p,k}(T,P)\f$:
00987    *
00988    *  \f[
00989    *     \bar C_{p,k}(T,P) =  C^{\triangle}_{p,k}(T,P) 
00990    *             - 2 R T \frac{d \ln( \gamma^{\triangle}_k)}{dT}
00991    *                    - R T^2 \frac{d^2 \ln(\gamma^{\triangle}_k) }{{dT}^2}
00992    * \f]
00993    * \f[
00994    *      \bar C_{p,o}(T,P) = C^o_{p,o}(T,P) 
00995    *                   - 2 R T \frac{d \ln(a_o)}{dT}
00996    *                    - R T^2 \frac{d^2 \ln(a_o)}{{dT}^2}
00997    * \f]
00998    *
00999    *  The pressure dependence of the activity coefficients leads to non-zero terms
01000    *  for the excess Volume of the solution.
01001    *  Therefore, the partial molar volumes are functions
01002    *  of the pressure derivatives of the activity coefficients.
01003    * \f[
01004    *     \bar V_k(T,P) =  V^{\triangle}_k(T,P) 
01005    *                    + R T \frac{d \ln(\gamma^{\triangle}_k) }{dP}
01006    * \f]
01007    * \f[
01008    *      \bar V_o(T,P) = V^o_o(T,P) 
01009    *                    + R T \frac{d \ln(a_o)}{dP}
01010    * \f]
01011    *
01012    *  The majority of work for these functions take place in the internal
01013    *  routines that calculate the first and second derivatives of the log
01014    *  of the activity coefficients wrt temperature,
01015    *  s_update_dlnMolalityActCoeff_dT(), s_update_d2lnMolalityActCoeff_dT2(),
01016    *  and the first 
01017    *  derivative of the log activity coefficients wrt pressure,
01018    *  s_update_dlnMolalityActCoeff_dP().
01019    *
01020    * <HR>
01021    * <H2> %Application within %Kinetics Managers </H2>
01022    * <HR>
01023    *
01024    * For the time being, we have set the standard concentration for all solute 
01025    * species in
01026    * this phase equal to the default concentration of the solvent at the system temperature
01027    * and pressure multiplied by Mnaught (kg solvent / gmol solvent). The solvent
01028    * standard concentration is just equal to its standard state concentration.
01029    *
01030    * 
01031    * This means that the
01032    * kinetics operator essentially works on an generalized concentration basis (kmol / m3),
01033    * with units for the kinetic rate constant specified
01034    * as if all reactants (solvent or solute) are on a concentration basis (kmol /m3).
01035    * The concentration will be modified by the activity coefficients.
01036    *
01037    * For example, a bulk-phase binary reaction between liquid solute species 
01038    * <I>j</I> and <I>k</I>, producing
01039    * a new liquid solute species <I>l</I> would have the
01040    * following equation for its rate of progress variable, \f$ R^1 \f$, which has
01041    * units of kmol m-3 s-1.
01042    *
01043    *   \f[
01044    *    R^1 = k^1 C_j^a C_k^a =  k^1 (C^o_o \tilde{M}_o a_j) (C^o_o \tilde{M}_o a_k) 
01045    *   \f]
01046    *
01047    * where
01048    *
01049    *   \f[
01050    *      C_j^a = C^o_o \tilde{M}_o a_j \quad and \quad C_k^a = C^o_o \tilde{M}_o a_k
01051    *   \f]
01052    *   
01053    *  \f$ C_j^a \f$ is the activity concentration of species <I>j</I>, and 
01054    *  \f$ C_k^a \f$ is the activity concentration of species <I>k</I>. \f$ C^o_o \f$
01055    *  is the concentration of water at 298 K and 1 atm. \f$ \tilde{M}_o \f$ 
01056    *  has units of kg solvent per gmol solvent and is equal to
01057    * 
01058    * \f[
01059    *     \tilde{M}_o = \frac{M_o}{1000} 
01060    * \f]
01061    *
01062    * \f$ a_j \f$ is
01063    *  the activity of species <I>j</I> at the current temperature and pressure
01064    *  and concentration of the liquid phase is given by the molality based
01065    *  activity coefficient multiplied by the molality of the jth species.
01066    *
01067    * \f[
01068    *      a_j  =  \gamma_j^\triangle m_j = \gamma_j^\triangle \frac{n_j}{\tilde{M}_o n_o}
01069    * \f]
01070    *
01071    * \f$k^1 \f$ has units of m<SUP>3</SUP>  kmol<SUP>-1</SUP> s<SUP>-1</SUP>.
01072    *
01073    *  Therefore the generalized activity concentration of a solute species has the following form
01074    *
01075    *  \f[
01076    *      C_j^a = C^o_o \frac{\gamma_j^\triangle n_j}{n_o}
01077    *  \f]
01078    *
01079    *  The generalized activity concentration of the solvent has the same units, but its a simpler form
01080    *   
01081    *  \f[
01082    *      C_o^a = C^o_o a_o
01083    *  \f]
01084    *
01085    *  The reverse rate constant can then be obtained from the law of microscopic reversibility
01086    *  and the equilibrium expression for the system.
01087    *
01088    *   \f[
01089    *      \frac{a_j a_k}{ a_l} = K^{o,1} = \exp(\frac{\mu^o_l - \mu^o_j - \mu^o_k}{R T} )
01090    *   \f]
01091    *
01092    *  \f$ K^{o,1} \f$ is the dimensionless form of the equilibrium constant.
01093    *  
01094    *   \f[
01095    *       R^{-1} = k^{-1} C_l^a =  k^{-1} (C_o  \tilde{M}_o a_l)
01096    *   \f]
01097    *
01098    *  where
01099    *
01100    *   \f[
01101    *       k^{-1} =  k^1 K^{o,1} C_o \tilde{M}_o
01102    *   \f]
01103    *
01104    *  \f$ k^{-1} \f$ has units of s<SUP>-1</SUP>.
01105    * 
01106    *  Note, this treatment may be modified in the future, as events dictate.
01107    *
01108    * <HR>
01109    * <H2> Instantiation of the Class </H2>
01110    * <HR>
01111    *
01112    * The constructor for this phase is now located in the default ThermoFactory
01113    * for %Cantera. The following code snipet may be used to initialize the phase
01114    * using the default construction technique within %Cantera.
01115    *
01116    * @code
01117    *      ThermoPhase *HMW = newPhase("HMW_NaCl.xml", "NaCl_electrolyte");
01118    * @endcode
01119    *
01120    *
01121    * A new %HMWSoln object may be created by  the following code snippets:
01122    *
01123    * @code
01124    *      HMWSoln *HMW = new HMWSoln("HMW_NaCl.xml", "NaCl_electrolyte");
01125    * @endcode
01126    *
01127    * or
01128    *
01129    * @code
01130    *    char iFile[80], file_ID[80];
01131    *    strcpy(iFile, "HMW_NaCl.xml");
01132    *    sprintf(file_ID,"%s#NaCl_electrolyte", iFile);
01133    *    XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
01134    *    HMWSoln *dh = new HMWSoln(*xm);
01135    * @endcode
01136    *
01137    * or by the following call to importPhase():
01138    *
01139    * @code
01140    *    char iFile[80], file_ID[80];
01141    *    strcpy(iFile, "HMW_NaCl.xml");
01142    *    sprintf(file_ID,"%s#NaCl_electrolyte", iFile);
01143    *    XML_Node *xm = get_XML_NameID("phase", file_ID, 0);
01144    *    HMWSoln dhphase;
01145    *    importPhase(*xm, &dhphase);
01146    * @endcode
01147    *
01148    *
01149    * <HR>
01150    * <H2> XML Example </H2>
01151    * <HR>
01152    *
01153    * The phase model name for this is called StoichSubstance. It must be supplied
01154    * as the model attribute of the thermo XML element entry.
01155    * Within the phase XML block,
01156    * the density of the phase must be specified. An example of an XML file
01157    * this phase is given below. 
01158    * 
01159    * @verbatim
01160    <phase id="NaCl_electrolyte" dim="3">
01161     <speciesArray datasrc="#species_waterSolution">
01162                H2O(L) Na+ Cl- H+ OH-
01163     </speciesArray>
01164     <state>
01165       <temperature units="K"> 300  </temperature>
01166       <pressure units="Pa">101325.0</pressure>
01167       <soluteMolalities>
01168              Na+:3.0
01169              Cl-:3.0
01170              H+:1.0499E-8
01171              OH-:1.3765E-6
01172       </soluteMolalities>
01173     </state>
01174     <!-- thermo model identifies the inherited class
01175          from ThermoPhase that will handle the thermodynamics.
01176       -->
01177     <thermo model="HMW">
01178        <standardConc model="solvent_volume" />
01179      <activityCoefficients model="Pitzer" TempModel="complex1">
01180                 <!-- Pitzer Coefficients
01181                      These coefficients are from Pitzer's main
01182                      paper, in his book.
01183                   -->
01184                 <A_Debye model="water" />
01185                 <ionicRadius default="3.042843"  units="Angstroms">
01186                 </ionicRadius>
01187                 <binarySaltParameters cation="Na+" anion="Cl-">
01188                   <beta0> 0.0765, 0.008946, -3.3158E-6,
01189                           -777.03, -4.4706
01190                   </beta0>
01191                   <beta1> 0.2664, 6.1608E-5, 1.0715E-6 </beta1>
01192                   <beta2> 0.0    </beta2>
01193                   <Cphi> 0.00127, -4.655E-5, 0.0,
01194                          33.317, 0.09421
01195                   </Cphi>
01196                   <Alpha1> 2.0 </Alpha1>
01197                 </binarySaltParameters>
01198 
01199                 <binarySaltParameters cation="H+" anion="Cl-">
01200                   <beta0> 0.1775, 0.0, 0.0, 0.0, 0.0</beta0>
01201                   <beta1> 0.2945, 0.0, 0.0 </beta1>
01202                   <beta2> 0.0    </beta2>
01203                   <Cphi> 0.0008, 0.0, 0.0, 0.0, 0.0 </Cphi>
01204                   <Alpha1> 2.0 </Alpha1>
01205                 </binarySaltParameters>
01206 
01207                 <binarySaltParameters cation="Na+" anion="OH-">
01208                   <beta0> 0.0864, 0.0, 0.0, 0.0, 0.0 </beta0>
01209                   <beta1> 0.253, 0.0, 0.0 </beta1>
01210                   <beta2> 0.0    </beta2>
01211                   <Cphi> 0.0044, 0.0, 0.0, 0.0, 0.0 </Cphi>
01212                   <Alpha1> 2.0 </Alpha1>
01213                 </binarySaltParameters>
01214 
01215                 <thetaAnion anion1="Cl-" anion2="OH-">
01216                   <Theta> -0.05 </Theta>
01217                 </thetaAnion>
01218 
01219                 <psiCommonCation cation="Na+" anion1="Cl-" anion2="OH-">
01220                   <Theta> -0.05 </Theta>
01221                   <Psi> -0.006 </Psi>
01222                 </psiCommonCation>
01223 
01224                 <thetaCation cation1="Na+" cation2="H+">
01225                   <Theta> 0.036 </Theta>
01226                 </thetaCation>
01227 
01228                 <psiCommonAnion anion="Cl-" cation1="Na+" cation2="H+">
01229                   <Theta> 0.036 </Theta>
01230                   <Psi> -0.004 </Psi>
01231                 </psiCommonAnion>
01232 
01233        </activityCoefficients>
01234 
01235        <solvent> H2O(L) </solvent>
01236     </thermo>
01237     <elementArray datasrc="elements.xml"> O H Na Cl </elementArray>
01238     <kinetics model="none" >
01239     </kinetics>
01240   </phase> 
01241   @endverbatim
01242    *
01243    *
01244    *
01245    * @ingroup thermoprops
01246    *
01247    */
01248   class HMWSoln : public MolalityVPSSTP {
01249 
01250   public:
01251         
01252     //! Default Constructor 
01253     HMWSoln();
01254 
01255     //! Construct and initialize an HMWSoln ThermoPhase object 
01256     //! directly from an asci input file
01257     /*!
01258      * Working constructors
01259      *
01260      *  The two constructors below are the normal way
01261      *  the phase initializes itself. They are shells that call
01262      *  the routine initThermo(), with a reference to the
01263      *  XML database to get the info for the phase.
01264      *
01265      * @param inputFile Name of the input file containing the phase XML data
01266      *                  to set up the object
01267      * @param id        ID of the phase in the input file. Defaults to the
01268      *                  empty string.
01269      */
01270     HMWSoln(std::string inputFile, std::string id = "");
01271 
01272     //! Construct and initialize an HMWSoln ThermoPhase object 
01273     //! directly from an XML database
01274     /*!
01275      *  @param phaseRef XML phase node containing the description of the phase
01276      *  @param id     id attribute containing the name of the phase. 
01277      *                (default is the empty string)
01278      */
01279     HMWSoln(XML_Node& phaseRef, std::string id = "");
01280 
01281 
01282     //! Copy Constructor
01283     /*!
01284      * Copy constructor for the object. Constructed
01285      * object will be a clone of this object, but will
01286      * also own all of its data.
01287      * This is a wrapper around the assignment operator
01288      * 
01289      * @param right Object to be copied.
01290      */
01291     HMWSoln(const HMWSoln &right);
01292 
01293     //! Asignment operator
01294     /*!
01295      * Assignment operator for the object. Constructed
01296      * object will be a clone of this object, but will
01297      * also own all of its data.
01298      * 
01299      * @param right Object to be copied.
01300      */
01301     HMWSoln& operator=(const HMWSoln& right);
01302 
01303     
01304     //!  This is a special constructor, used to replicate test problems
01305     //!  during the initial verification of the object
01306     /*!  
01307      *
01308      *
01309      *  test problems:
01310      *  1 = NaCl problem - 5 species -
01311      *   the thermo is read in from an XML file
01312      *
01313      * speci   molality                        charge
01314      *  Cl-     6.0954          6.0997E+00      -1
01315      *  H+      1.0000E-08      2.1628E-09      1
01316      *  Na+     6.0954E+00      6.0997E+00      1
01317      *  OH-     7.5982E-07      1.3977E-06     -1
01318      *  HMW_params____beta0MX__beta1MX__beta2MX__CphiMX_____alphaMX__thetaij
01319      * 10
01320      * 1  2          0.1775  0.2945   0.0      0.00080    2.0      0.0
01321      * 1  3          0.0765  0.2664   0.0      0.00127    2.0      0.0
01322      * 1  4          0.0     0.0      0.0      0.0        0.0     -0.050
01323      * 2  3          0.0     0.0      0.0      0.0        0.0      0.036
01324      * 2  4          0.0     0.0      0.0      0.0        0.0      0.0
01325      * 3  4          0.0864  0.253    0.0      0.0044     2.0      0.0
01326      * Triplet_interaction_parameters_psiaa'_or_psicc'
01327      * 2
01328      * 1  2  3   -0.004
01329      * 1  3  4   -0.006
01330      *
01331      * @param testProb Hard -coded test problem to instantiate.
01332      *                 Current valid values are 1.
01333      */
01334     HMWSoln(int testProb);
01335 
01336     //! Destructor. 
01337     virtual ~HMWSoln();
01338 
01339     //! Duplicator from the ThermoPhase parent class
01340     /*!
01341      * Given a pointer to a ThermoPhase object, this function will
01342      * duplicate the ThermoPhase object and all underlying structures.
01343      * This is basically a wrapper around the copy constructor.
01344      *
01345      * @return returns a pointer to a ThermoPhase
01346      */
01347     ThermoPhase *duplMyselfAsThermoPhase() const;
01348 
01349     /** 
01350      * @name  Utilities  
01351      * @{
01352      */
01353 
01354     /** 
01355      * Equation of state type flag. The base class returns
01356      * zero. Subclasses should define this to return a unique
01357      * non-zero value. Constants defined for this purpose are
01358      * listed in mix_defs.h.
01359      */
01360     virtual int eosType() const;
01361 
01362     /**
01363      * @} 
01364      * @name  Molar Thermodynamic Properties of the Solution --------------
01365      * @{
01366      */
01367 
01368     /// Molar enthalpy. Units: J/kmol. 
01369     /**
01370      * Molar enthalpy of the solution. Units: J/kmol.
01371      *      (HKM -> Bump up to Parent object)
01372      */
01373     virtual doublereal enthalpy_mole() const;
01374 
01375     /**
01376      * Excess molar enthalpy of the solution from 
01377      * the mixing process. Units: J/ kmol.
01378      *
01379      * Note this is kmol of the total solution.
01380      */
01381     virtual doublereal relative_enthalpy() const;
01382 
01383     /**
01384      * Excess molar enthalpy of the solution from 
01385      * the mixing process on a molality basis.
01386      *  Units: J/ (kmol add salt).
01387      *
01388      * Note this is kmol of the guessed at salt composition
01389      */
01390     virtual doublereal relative_molal_enthalpy() const;
01391 
01392 
01393     /// Molar internal energy. Units: J/kmol. 
01394     /**
01395      * Molar internal energy of the solution. Units: J/kmol.
01396      *      (HKM -> Bump up to Parent object)
01397      */
01398     virtual doublereal intEnergy_mole() const;
01399 
01400     /// Molar entropy. Units: J/kmol/K. 
01401     /**
01402      * Molar entropy of the solution. Units: J/kmol/K.
01403      * For an ideal, constant partial molar volume solution mixture with
01404      * pure species phases which exhibit zero volume expansivity:
01405      * \f[
01406      * \hat s(T, P, X_k) = \sum_k X_k \hat s^0_k(T)  
01407      *      - \hat R  \sum_k X_k log(X_k)
01408      * \f]
01409      * The reference-state pure-species entropies 
01410      * \f$ \hat s^0_k(T,p_{ref}) \f$ are computed by the
01411      *  species thermodynamic 
01412      * property manager. The pure species entropies are independent of 
01413      * temperature since the volume expansivities are equal to zero.
01414      * @see SpeciesThermo
01415      *
01416      *      (HKM -> Bump up to Parent object)
01417      */
01418     virtual doublereal entropy_mole() const;
01419 
01420     /// Molar Gibbs function. Units: J/kmol. 
01421     /*
01422      *      (HKM -> Bump up to Parent object)
01423      */
01424     virtual doublereal gibbs_mole() const;
01425 
01426     /// Molar heat capacity at constant pressure. Units: J/kmol/K. 
01427     /*
01428      *     
01429      */
01430     virtual doublereal cp_mole() const;
01431 
01432     /// Molar heat capacity at constant volume. Units: J/kmol/K.
01433     /*
01434      *      (HKM -> Bump up to Parent object)
01435      */
01436     virtual doublereal cv_mole() const;
01437 
01438     //@}
01439     /// @name Mechanical Equation of State Properties ---------------------
01440     //@{
01441     /**
01442      *   In this equation of state implementation, the density is a 
01443      *   function only of the mole fractions. Therefore, it can't be 
01444      *   an independent variable. Instead, the pressure is used as the
01445      *   independent variable. Functions which try to set the thermodynamic
01446      *   state by calling setDensity() may cause an exception to be
01447      *   thrown.  
01448      */
01449 
01450     /**
01451      * Pressure. Units: Pa.
01452      * For this incompressible system, we return the internally storred
01453      * independent value of the pressure.
01454      */ 
01455     virtual doublereal pressure() const;
01456 
01457     //! Set the internally storred pressure (Pa) at constant
01458     //! temperature and composition
01459     /*!
01460      *  This method sets the pressure within the object.
01461      *  The water model is a completely compressible model.
01462      *  Also, the dielectric constant is pressure dependent.
01463      *
01464      *  @param p input Pressure (Pa)
01465      *
01466      * @todo Implement a variable pressure capability
01467      */
01468     virtual void setPressure(doublereal p);
01469 
01470   protected:
01471     /**
01472      * Calculate the density of the mixture using the partial 
01473      * molar volumes and mole fractions as input
01474      *
01475      * The formula for this is
01476      *
01477      * \f[ 
01478      * \rho = \frac{\sum_k{X_k W_k}}{\sum_k{X_k V_k}} 
01479      * \f]
01480      *
01481      * where \f$X_k\f$ are the mole fractions, \f$W_k\f$ are
01482      * the molecular weights, and \f$V_k\f$ are the pure species
01483      * molar volumes.
01484      *
01485      * Note, the basis behind this formula is that in an ideal
01486      * solution the partial molar volumes are equal to the pure
01487      * species molar volumes. We have additionally specified
01488      * in this class that the pure species molar volumes are
01489      * independent of temperature and pressure.
01490      *
01491      * NOTE: This is a non-virtual function, which is not a 
01492      *       member of the ThermoPhase base class. 
01493      */
01494     void calcDensity();
01495 
01496   public:
01497     //! Returns the current value of the density
01498     /*!
01499      *  @return value of the density. Units: kg/m^3
01500      */
01501     virtual doublereal density() const;
01502 
01503     //! Set the internally storred density (kg/m^3) of the phase.
01504     /*!
01505      * Overwritten setDensity() function is necessary because of
01506      * the underlying water model.
01507      *
01508      * @todo Now have a compressible ss equation for liquid water.
01509      *       Therefore, this phase is compressible. May still
01510      *       want to change the independent variable however. 
01511      *
01512      *  NOTE: This is an overwritten function from the State.h
01513      *        class
01514      * 
01515      * @param rho Input density (kg/m^3).
01516      */
01517     void setDensity(const doublereal rho);
01518 
01519     //! Set the internally storred molar density (kmol/m^3) for the phase.
01520     /**
01521      * Overwritten setMolarDensity() function is necessary because of the
01522      * underlying water model.
01523      *
01524      * This function will now throw an error condition if the input
01525      * isn't exactly equal to the current molar density.
01526      *
01527      *  NOTE: This is a virtual function overwritten from the State.h
01528      *        class
01529      *
01530      * @param conc   Input molar density (kmol/m^3).
01531      */
01532     void setMolarDensity(const doublereal conc);
01533 
01534     //! Set the temperature (K)
01535     /*!
01536      * Overwritten setTemperature(double) from State.h. This
01537      * function sets the temperature, and makes sure that
01538      * the value propagates to underlying objects, such as
01539      * the water standard state model.
01540      *
01541      * @todo Make State::setTemperature a virtual function
01542      *
01543      * @param temp Temperature in kelvin
01544      */
01545     virtual void setTemperature(const doublereal temp);
01546 
01547     //! Set the temperature (K) and pressure (Pa)
01548     /*!
01549      *  Set the temperature and pressure.
01550      *
01551      * @param t    Temperature (K)
01552      * @param p    Pressure (Pa)
01553      */
01554     virtual void setState_TP(doublereal t, doublereal p);
01555 
01556 
01557     /**
01558      * The isothermal compressibility. Units: 1/Pa.
01559      * The isothermal compressibility is defined as
01560      * \f[
01561      * \kappa_T = -\frac{1}{v}\left(\frac{\partial v}{\partial P}\right)_T
01562      * \f]
01563      */
01564     virtual doublereal isothermalCompressibility() const;
01565 
01566     /**
01567      * The thermal expansion coefficient. Units: 1/K.
01568      * The thermal expansion coefficient is defined as
01569      *
01570      * \f[
01571      * \beta = \frac{1}{v}\left(\frac{\partial v}{\partial T}\right)_P
01572      * \f]
01573      */
01574     virtual doublereal thermalExpansionCoeff() const;
01575 
01576     /**
01577      * @} 
01578      * @name Potential Energy
01579      * 
01580      * Species may have an additional potential energy due to the
01581      * presence of external gravitation or electric fields. These
01582      * methods allow specifying a potential energy for individual
01583      * species.
01584      * @{
01585      */
01586 
01587 
01588 
01589     /**
01590      * @}
01591      * @name Activities, Standard States,  and Activity Concentrations
01592      *
01593      * The activity \f$a_k\f$ of a species in solution is
01594      * related to the chemical potential by \f[ \mu_k = \mu_k^0(T)
01595      * + \hat R T \log a_k. \f] The quantity \f$\mu_k^0(T,P)\f$ is
01596      * the chemical potential at unit activity, which depends only
01597      * on temperature and the pressure.
01598      * Activity is assumed to be molality-based here.
01599      * @{
01600      */
01601 
01602     
01603     //! This method returns an array of generalized activity concentrations
01604     /*!
01605      * The generalized activity concentrations, \f$ C_k^a\f$,  are defined such that 
01606      * \f$ a_k = C^a_k / C^0_k, \f$ where \f$ C^0_k \f$ 
01607      * is a standard concentration
01608      * defined below.  These generalized concentrations are used
01609      * by kinetics manager classes to compute the forward and
01610      * reverse rates of elementary reactions. 
01611      *
01612      *  The generalized activity concentration of a solute species has the following form
01613      *
01614      *  \f[
01615      *      C_j^a = C^o_o \frac{\gamma_j^\triangle n_j}{n_o}
01616      *  \f]
01617      *
01618      *  The generalized activity concentration of the solvent has the same units, but its a simpler form
01619      *   
01620      *  \f[
01621      *      C_o^a = C^o_o a_o
01622      *  \f]
01623      *
01624      *
01625      * @param c Array of generalized concentrations. The 
01626      *          units are kmol m-3 for both the solvent and the solute species
01627      */
01628     virtual void getActivityConcentrations(doublereal* c) const;
01629 
01630     //! Return the standard concentration for the kth species
01631     /*!
01632      * The standard concentration \f$ C^0_k \f$ used to normalize
01633      * the activity (i.e., generalized) concentration for use
01634      *
01635      * We have set the standard concentration for all solute  species in
01636      * this phase equal to the default concentration of the solvent at the system temperature
01637      * and pressure multiplied by Mnaught (kg solvent / gmol solvent). The solvent
01638      * standard concentration is just equal to its standard state concentration.
01639      *
01640      *   \f[
01641      *      C_j^0 = C^o_o \tilde{M}_o \quad and  C_o^0 = C^o_o 
01642      *   \f]
01643      *
01644      * The consequence of this is that the standard concentrations have unequal units
01645      * between the solvent and the solute. However, both the solvent and the solute
01646      * activity concentrations will have the same units of kmol kg<SUP>-3</SUP>.
01647      *
01648      * This means that the
01649      * kinetics operator essentially works on an generalized concentration basis (kmol / m3),
01650      * with units for the kinetic rate constant specified
01651      * as if all reactants (solvent or solute) are on a concentration basis (kmol /m3).
01652      * The concentration will be modified by the activity coefficients.
01653      *
01654      * For example, a bulk-phase binary reaction between liquid solute species 
01655      * <I>j</I> and <I>k</I>, producing
01656      * a new liquid solute species <I>l</I> would have the
01657      * following equation for its rate of progress variable, \f$ R^1 \f$, which has
01658      * units of kmol m-3 s-1.
01659      *
01660      *   \f[
01661      *    R^1 = k^1 C_j^a C_k^a =  k^1 (C^o_o \tilde{M}_o a_j) (C^o_o \tilde{M}_o a_k) 
01662      *   \f]
01663      *
01664      * where
01665      *
01666      *   \f[
01667      *      C_j^a = C^o_o \tilde{M}_o a_j \quad and \quad C_k^a = C^o_o \tilde{M}_o a_k
01668      *   \f]
01669      *   
01670      * \f$ C_j^a \f$ is the activity concentration of species <I>j</I>, and 
01671      * \f$ C_k^a \f$ is the activity concentration of species <I>k</I>. \f$ C^o_o \f$
01672      * is the concentration of water at 298 K and 1 atm. \f$ \tilde{M}_o \f$ 
01673      * has units of kg solvent per gmol solvent and is equal to
01674      * 
01675      * \f[
01676      *     \tilde{M}_o = \frac{M_o}{1000} 
01677      * \f]
01678      *
01679      * \f$ a_j \f$ is
01680      *  the activity of species <I>j</I> at the current temperature and pressure
01681      *  and concentration of the liquid phase is given by the molality based
01682      *  activity coefficient multiplied by the molality of the jth species.
01683      *
01684      * \f[
01685      *      a_j  =  \gamma_j^\triangle m_j = \gamma_j^\triangle \frac{n_j}{\tilde{M}_o n_o}
01686      * \f]
01687      *
01688      * \f$k^1 \f$ has units of m<SUP>3</SUP>  kmol<SUP>-1</SUP> s<SUP>-1</SUP>.
01689      *
01690      *  Therefore the generalized activity concentration of a solute species has the following form
01691      *
01692      *  \f[
01693      *      C_j^a = C^o_o \frac{\gamma_j^\triangle n_j}{n_o}
01694      *  \f]
01695      *
01696      *  The generalized activity concentration of the solvent has the same units, but its a simpler form
01697      *   
01698      *  \f[
01699      *      C_o^a = C^o_o a_o
01700      *  \f]
01701      *
01702      *
01703      * @param k Optional parameter indicating the species. The default
01704      *         is to assume this refers to species 0.
01705      * @return 
01706      *   Returns the standard Concentration in units of 
01707      *   m<SUP>3</SUP> kmol<SUP>-1</SUP>.
01708      *
01709      * @param k Species index
01710      */
01711     virtual doublereal standardConcentration(int k=0) const;
01712 
01713     
01714     //! Returns the natural logarithm of the standard 
01715     //! concentration of the kth species
01716     /*!
01717      * @param k Species index
01718      */
01719     virtual doublereal logStandardConc(int k=0) const;
01720 
01721     //! Returns the units of the standard and generalized concentrations.
01722     /*!
01723      * Note they have the same units, as their
01724      * ratio is defined to be equal to the activity of the kth
01725      * species in the solution, which is unitless.
01726      *
01727      * This routine is used in print out applications where the
01728      * units are needed. Usually, MKS units are assumed throughout
01729      * the program and in the XML input files.
01730      *
01731      * The base %ThermoPhase class assigns the default quantities
01732      * of (kmol/m3) for all species.
01733      * Inherited classes are responsible for overriding the default 
01734      * values if necessary.
01735      *
01736      * @param uA Output vector containing the units
01737      *  uA[0] = kmol units - default  = 1
01738      *  uA[1] = m    units - default  = -nDim(), the number of spatial
01739      *                                dimensions in the Phase class.
01740      *  uA[2] = kg   units - default  = 0;
01741      *  uA[3] = Pa(pressure) units - default = 0;
01742      *  uA[4] = Temperature units - default = 0;
01743      *  uA[5] = time units - default = 0
01744      * @param k species index. Defaults to 0.
01745      * @param sizeUA output int containing the size of the vector.
01746      *        Currently, this is equal to 6.
01747      */
01748     virtual void getUnitsStandardConc(double *uA, int k = 0,
01749                                       int sizeUA = 6) const;
01750 
01751     //! Get the array of non-dimensional activities at
01752     //! the current solution temperature, pressure, and solution concentration.
01753     /*!
01754      *
01755      * We resolve this function at this level by calling
01756      * on the activityConcentration function. However, 
01757      * derived classes may want to override this default
01758      * implementation.
01759      *
01760      * (note solvent is on molar scale).
01761      *
01762      * @param ac  Output vector of activities. Length: m_kk.
01763      */
01764     virtual void getActivities(doublereal* ac) const;
01765 
01766     //@}
01767     /// @name  Partial Molar Properties of the Solution -----------------
01768     //@{
01769 
01770     //! Get the species chemical potentials. Units: J/kmol.
01771     /*!
01772      *
01773      * This function returns a vector of chemical potentials of the 
01774      * species in solution.
01775      *
01776      * \f[
01777      *    \mu_k = \mu^{\triangle}_k(T,P) + R T ln(\gamma_k^{\triangle} m_k)
01778      * \f] 
01779      *
01780      * @param mu  Output vector of species chemical 
01781      *            potentials. Length: m_kk. Units: J/kmol
01782      */
01783     virtual void getChemPotentials(doublereal* mu) const;
01784    
01785     //! Returns an array of partial molar enthalpies for the species
01786     //! in the mixture. Units (J/kmol)
01787     /*!
01788      * For this phase, the partial molar enthalpies are equal to the
01789      * standard state enthalpies modified by the derivative of the
01790      * molality-based activity coefficent wrt temperature
01791      *
01792      *  \f[
01793      * \bar h_k(T,P) = h^{\triangle}_k(T,P) 
01794      *               - R T^2 \frac{d \ln(\gamma_k^\triangle)}{dT}
01795      * \f]
01796      * The solvent partial molar enthalpy is equal to 
01797      *  \f[
01798      * \bar h_o(T,P) = h^{o}_o(T,P) - R T^2 \frac{d \ln(a_o)}{dT}
01799      *       = h^{o}_o(T,P) 
01800      *       + R T^2 (\sum_{k \neq o} m_k)  \tilde{M_o} (\frac{d \phi}{dT})
01801      * \f]
01802      *
01803      *
01804      * @param hbar    Output vector of species partial molar enthalpies.
01805      *                Length: m_kk. units are J/kmol.
01806      */
01807     virtual void getPartialMolarEnthalpies(doublereal* hbar) const;
01808 
01809   
01810     //! Returns an array of partial molar entropies of the species in the
01811     //! solution. Units: J/kmol/K.
01812     /*!
01813      * Maxwell's equations provide an answer for how calculate this
01814      *   (p.215 Smith and Van Ness)
01815      *
01816      *      d(chemPot_i)/dT = -sbar_i
01817      *      
01818      * For this phase, the partial molar entropies are equal to the
01819      * SS species entropies plus the ideal solution contribution
01820      * plus complicated functions of the
01821      * temperature derivative of the activity coefficents.
01822      *
01823      *  \f[
01824      *     \bar s_k(T,P) =  s^{\triangle}_k(T,P) 
01825      *             - R \ln( \gamma^{\triangle}_k \frac{m_k}{m^{\triangle}}))
01826      *                    - R T \frac{d \ln(\gamma^{\triangle}_k) }{dT}
01827      * \f]
01828      * \f[
01829      *      \bar s_o(T,P) = s^o_o(T,P) - R \ln(a_o)
01830      *                    - R T \frac{d \ln(a_o)}{dT}
01831      * \f]
01832      *
01833      *  @param sbar    Output vector of species partial molar entropies.
01834      *                 Length = m_kk. units are J/kmol/K.
01835      */
01836     virtual void getPartialMolarEntropies(doublereal* sbar) const;
01837      
01838     //! Return an array of partial molar volumes for the
01839     //! species in the mixture. Units: m^3/kmol.
01840     /*!
01841      *  For this solution, the partial molar volumes are functions
01842      *  of the pressure derivatives of the activity coefficients.
01843      *
01844      *  \f[
01845      *     \bar V_k(T,P)  = V^{\triangle}_k(T,P) 
01846      *                    + R T \frac{d \ln(\gamma^{\triangle}_k) }{dP}
01847      * \f]
01848      * \f[
01849      *      \bar V_o(T,P) = V^o_o(T,P) 
01850      *                    + R T \frac{d \ln(a_o)}{dP}
01851      * \f]
01852      *
01853      *  @param vbar   Output vector of speciar partial molar volumes.
01854      *                Length = m_kk. units are m^3/kmol.
01855      */
01856     virtual void getPartialMolarVolumes(doublereal* vbar) const;
01857 
01858     //! Return an array of partial molar heat capacities for the
01859     //! species in the mixture.  Units: J/kmol/K
01860     /*!
01861      *  The following formulas are implemented within the code.
01862      *
01863      *  \f[
01864      *     \bar C_{p,k}(T,P) =  C^{\triangle}_{p,k}(T,P) 
01865      *             - 2 R T \frac{d \ln( \gamma^{\triangle}_k)}{dT}
01866      *                    - R T^2 \frac{d^2 \ln(\gamma^{\triangle}_k) }{{dT}^2}
01867      * \f]
01868      * \f[
01869      *      \bar C_{p,o}(T,P) = C^o_{p,o}(T,P) 
01870      *                   - 2 R T \frac{d \ln(a_o)}{dT}
01871      *                    - R T^2 \frac{d^2 \ln(a_o)}{{dT}^2}
01872      * \f]
01873      *
01874      *
01875      * @param cpbar   Output vector of species partial molar heat 
01876      *                capacities at constant pressure.
01877      *                Length = m_kk. units are J/kmol/K.
01878      */
01879     virtual void getPartialMolarCp(doublereal* cpbar) const;
01880 
01881     //@}
01882 
01883   
01884   
01885   protected:
01886 
01887     //! Updates the standard state thermodynamic functions at the current T and P of the solution.
01888     /*!
01889      * @internal
01890      *
01891      * This function gets called for every call to a public function in this
01892      * class. It checks to see whether the temperature or pressure has changed and
01893      * thus whether the ss thermodynamics functions must be recalculated.
01894      *
01895      * @param pres  Pressure at which to evaluate the standard states.
01896      *              The default, indicated by a -1.0, is to use the current pressure
01897      */                      
01898     //virtual void _updateStandardStateThermo() const;
01899 
01900     //@}
01901     /// @name Thermodynamic Values for the Species Reference States ---
01902     //@{
01903 
01904 
01905     ///////////////////////////////////////////////////////
01906     //
01907     //  The methods below are not virtual, and should not
01908     //  be overloaded.
01909     //
01910     //////////////////////////////////////////////////////
01911 
01912     /**
01913      * @name Specific Properties
01914      * @{
01915      */
01916 
01917 
01918     /**
01919      * @name Setting the State
01920      *
01921      * These methods set all or part of the thermodynamic
01922      * state.
01923      * @{
01924      */
01925 
01926     //@}
01927 
01928     /**
01929      * @name Chemical Equilibrium
01930      * Chemical equilibrium.
01931      * @{
01932      */
01933   public:
01934 
01935     //!This method is used by the ChemEquil equilibrium solver.
01936     /*!
01937      * It sets the state such that the chemical potentials satisfy
01938      * \f[ \frac{\mu_k}{\hat R T} = \sum_m A_{k,m}
01939      * \left(\frac{\lambda_m} {\hat R T}\right) \f] where 
01940      * \f$ \lambda_m \f$ is the element potential of element m. The
01941      * temperature is unchanged.  Any phase (ideal or not) that
01942      * implements this method can be equilibrated by ChemEquil.
01943      *
01944      * @param lambda_RT Input vector of dimensionless element potentials
01945      *                  The length is equal to nElements().
01946      */ 
01947     virtual void setToEquilState(const doublereal* lambda_RT) {
01948       updateStandardStateThermo();
01949       err("setToEquilState");
01950     }
01951 
01952     //@}
01953 
01954 
01955     //! Set the equation of state parameters
01956     /*!
01957      * @internal
01958      *  The number and meaning of these depends on the subclass. 
01959      *
01960      * @param n number of parameters
01961      * @param c array of \a n coefficients
01962      */
01963     virtual void setParameters(int n, doublereal* const c);
01964 
01965     //! Get the equation of state parameters in a vector
01966     /*!
01967      * @internal
01968      * The number and meaning of these depends on the subclass. 
01969      *
01970      * @param n number of parameters
01971      * @param c array of \a n coefficients
01972      */
01973     virtual void getParameters(int &n, doublereal * const c) const;
01974     
01975     //! Set equation of state parameter values from XML
01976     //! entries.
01977     /*!
01978      * This method is called by function importPhase in
01979      * file importCTML.cpp when processing a phase definition in
01980      * an input file. It should be overloaded in subclasses to set
01981      * any parameters that are specific to that particular phase
01982      * model. 
01983      *   
01984      * @param eosdata An XML_Node object corresponding to
01985      *                the "thermo" entry for this phase in the input file.
01986      */
01987     virtual void setParametersFromXML(const XML_Node& eosdata);
01988  
01989     //---------------------------------------------------------
01990     /// @name Critical state properties.
01991     /// These methods are only implemented by some subclasses.
01992         
01993     //@{
01994         
01995     /// Critical temperature (K).
01996     virtual doublereal critTemperature() const {
01997       err("critTemperature"); return -1.0;
01998     }
01999         
02000     /// Critical pressure (Pa).
02001     virtual doublereal critPressure() const {
02002       err("critPressure"); return -1.0;
02003     }
02004         
02005     /// Critical density (kg/m3).
02006     virtual doublereal critDensity() const {
02007       err("critDensity"); return -1.0;
02008     }                
02009         
02010     //@}
02011         
02012     /// @name Saturation properties.
02013     /// These methods are only implemented by subclasses that 
02014     /// implement full liquid-vapor equations of state.
02015     ///
02016     virtual doublereal satTemperature(doublereal p) const {
02017       err("satTemperature"); return -1.0;
02018     }
02019           
02020     //! Get the saturation pressure for a given temperature. 
02021     /*!
02022      * Note the limitations of this function. Stability considerations
02023      * concernting multiphase equilibrium are ignored in this 
02024      * calculation. Therefore, the call is made directly to the SS of 
02025      * water underneath. The object is put back into its original
02026      * state at the end of the call.
02027      *
02028      * @todo This is probably not implemented correctly. The stability
02029      *       of the salt should be added into this calculation. The
02030      *       underlying water model may be called to get the stability
02031      *       of the pure water solution, if needed.
02032      *
02033      * @param T  Temperature (kelvin)
02034      */
02035     virtual doublereal satPressure(doublereal T) const;
02036         
02037     virtual doublereal vaporFraction() const {
02038       err("vaprFraction"); return -1.0;
02039     }
02040         
02041     virtual void setState_Tsat(doublereal t, doublereal x) {
02042       err("setState_sat"); 
02043     }
02044 
02045     virtual void setState_Psat(doublereal p, doublereal x) {
02046       err("setState_sat"); 
02047     }
02048 
02049     //@}
02050 
02051 
02052     /*
02053      *  -------------- Utilities -------------------------------
02054      */
02055 
02056     /**
02057      * Return a reference to the species thermodynamic property
02058      * manager. 
02059      *
02060      * @todo This method will fail if no species thermo
02061      * manager has been installed.
02062      */
02063     SpeciesThermo& speciesThermo() { return *m_spthermo; }
02064 
02065     //! Initialization of a HMWSoln phase using an xml file
02066     /*!
02067      * This routine is a precursor to initThermo(XML_Node*)
02068      * routine, which does most of the work.
02069      *
02070      * @param inputFile XML file containing the description of the
02071      *        phase
02072      *
02073      * @param id  Optional parameter identifying the name of the
02074      *            phase. If none is given, the first XML
02075      *            phase element will be used.
02076      */
02077     void constructPhaseFile(std::string inputFile, std::string id);
02078 
02079     //!   Import and initialize a HMWSoln phase 
02080     //!   specification in an XML tree into the current object.
02081     /*!
02082      *   Here we read an XML description of the phase.
02083      *   We import descriptions of the elements that make up the
02084      *   species in a phase.
02085      *   We import information about the species, including their
02086      *   reference state thermodynamic polynomials. We then freeze
02087      *   the state of the species.
02088      *
02089      *   Then, we read the species molar volumes from the xml 
02090      *   tree to finish the initialization.
02091      *
02092      * @param phaseNode This object must be the phase node of a
02093      *             complete XML tree
02094      *             description of the phase, including all of the
02095      *             species data. In other words while "phase" must
02096      *             point to an XML phase object, it must have
02097      *             sibling nodes "speciesData" that describe
02098      *             the species in the phase.
02099      *
02100      * @param id   ID of the phase. If nonnull, a check is done
02101      *             to see if phaseNode is pointing to the phase
02102      *             with the correct id. 
02103      */
02104     void constructPhaseXML(XML_Node& phaseNode, std::string id);
02105 
02106     //! Internal initialization required after all species have
02107     //! been added
02108     /*!
02109      * @internal Initialize. This method is provided to allow
02110      * subclasses to perform any initialization required after all
02111      * species have been added. For example, it might be used to
02112      * resize internal work arrays that must have an entry for
02113      * each species.  The base class implementation does nothing,
02114      * and subclasses that do not require initialization do not
02115      * need to overload this method.  When importing a CTML phase
02116      * description, this method is called just prior to returning
02117      * from function importPhase.
02118      *
02119      * @see importCTML.cpp
02120      */
02121     virtual void initThermo();
02122 
02123     //! Initialize the phase parameters from an XML file.
02124     /*!
02125      *  initThermoXML()                 (virtual from ThermoPhase)
02126      *
02127      *  This gets called from importPhase(). It processes the XML file
02128      *  after the species are set up. This is the main routine for
02129      *  reading in activity coefficient parameters.
02130      *
02131      * @param phaseNode This object must be the phase node of a
02132      *             complete XML tree
02133      *             description of the phase, including all of the
02134      *             species data. In other words while "phase" must
02135      *             point to an XML phase object, it must have
02136      *             sibling nodes "speciesData" that describe
02137      *             the species in the phase.
02138      * @param id   ID of the phase. If nonnull, a check is done
02139      *             to see if phaseNode is pointing to the phase
02140      *             with the correct id.
02141      */
02142     virtual void initThermoXML(XML_Node& phaseNode, std::string id);
02143 
02144     //!  Report the molar volume of species k
02145     /*!
02146      *
02147      * units - \f$ m^3 kmol^-1 \f$
02148      *
02149      * @param k species index
02150      *
02151      * @deprecated
02152      *   The getPartialMolarVolumes() expression is more precise.
02153      */
02154     double speciesMolarVolume(int k) const;
02155 
02156   
02157     //! Value of the Debye Huckel constant as a function of temperature
02158     //! and pressure.
02159     /*!
02160      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02161      *
02162      *            Units = sqrt(kg/gmol)
02163      *
02164      * @param temperature  Temperature of the derivative calculation
02165      *                     or -1 to indicate the current temperature
02166      *
02167      * @param pressure    Pressure of the derivative calcualtion
02168      *                    or -1 to indicate the current pressure
02169      */
02170     virtual double A_Debye_TP(double temperature = -1.0, 
02171                               double pressure = -1.0) const;
02172     
02173     //! Value of the derivative of the Debye Huckel constant with 
02174     //! respect to temperature as a function of temperature
02175     //! and pressure.
02176     /*!
02177      *
02178      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02179      *
02180      *            Units = sqrt(kg/gmol)
02181      *
02182      * @param temperature  Temperature of the derivative calculation
02183      *                     or -1 to indicate the current temperature
02184      *
02185      * @param pressure    Pressure of the derivative calcualtion
02186      *                    or -1 to indicate the current pressure
02187      */
02188     virtual double dA_DebyedT_TP(double temperature = -1.0, 
02189                                  double pressure = -1.0) const;
02190 
02191     /**
02192      * Value of the derivative of the Debye Huckel constant with 
02193      * respect to pressure, as a function of temperature
02194      * and pressure.
02195      *
02196      *      A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02197      *
02198      *  Units = sqrt(kg/gmol)
02199      *
02200      * @param temperature  Temperature of the derivative calculation
02201      *                     or -1 to indicate the current temperature
02202      *
02203      * @param pressure    Pressure of the derivative calcualtion
02204      *                    or -1 to indicate the current pressure
02205      */
02206     virtual double dA_DebyedP_TP(double temperature = -1.0, 
02207                                  double pressure = -1.0) const;
02208 
02209     /**
02210      * Return Pitzer's definition of A_L. This is basically the
02211      * derivative of the A_phi multiplied by 4 R T**2
02212      *
02213      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02214      *            dA_phidT = d(A_Debye)/dT / 3.0
02215      *            A_L = dA_phidT * (4 * R * T * T)
02216      *
02217      *            Units = sqrt(kg/gmol) (RT)
02218      * 
02219      *
02220      * @param temperature  Temperature of the derivative calculation
02221      *                     or -1 to indicate the current temperature
02222      *
02223      * @param pressure    Pressure of the derivative calcualtion
02224      *                    or -1 to indicate the current pressure
02225      */
02226     double ADebye_L(double temperature = -1.0,
02227                     double pressure = -1.0) const;
02228 
02229     /**
02230      * Return Pitzer's definition of A_J. This is basically the
02231      * temperature derivative of A_L, and the second derivative
02232      * of A_phi
02233      *
02234      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02235      *            dA_phidT = d(A_Debye)/dT / 3.0
02236      *            A_J = 2 A_L/T + 4 * R * T * T * d2(A_phi)/dT2
02237      *
02238      *            Units = sqrt(kg/gmol) (R)
02239      *
02240      * @param temperature  Temperature of the derivative calculation
02241      *                     or -1 to indicate the current temperature
02242      *
02243      * @param pressure    Pressure of the derivative calcualtion
02244      *                    or -1 to indicate the current pressure
02245      */
02246     double ADebye_J(double temperature = -1.0,
02247                     double pressure = -1.0) const;
02248     /**
02249      * Return Pitzer's definition of A_V. This is the
02250      * derivative wrt pressure of A_phi multiplied by - 4 R T
02251      *
02252      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02253      *            dA_phidT = d(A_Debye)/dP / 3.0
02254      *            A_V = - dA_phidP * (4 * R * T)
02255      *
02256      *            Units = sqrt(kg/gmol) (RT) / Pascal
02257      *
02258      * @param temperature  Temperature of the derivative calculation
02259      *                     or -1 to indicate the current temperature
02260      *
02261      * @param pressure    Pressure of the derivative calcualtion
02262      *                    or -1 to indicate the current pressure
02263      * 
02264      */
02265     double ADebye_V(double temperature = -1.0,
02266                     double pressure = -1.0) const;
02267     
02268     //! Value of the 2nd derivative of the Debye Huckel constant with 
02269     //! respect to temperature as a function of temperature
02270     //! and pressure.
02271     /*!
02272      *
02273      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02274      *
02275      *            Units = sqrt(kg/gmol)
02276      *
02277      * @param temperature  Temperature of the derivative calculation
02278      *                     or -1 to indicate the current temperature
02279      *
02280      * @param pressure    Pressure of the derivative calcualtion
02281      *                    or -1 to indicate the current pressure
02282      */
02283     virtual double d2A_DebyedT2_TP(double temperature = -1.0, 
02284                                    double pressure = -1.0) const;
02285 
02286     //! Reports the ionic radius of the kth species
02287     /*!
02288      *  @param k Species index 
02289      */
02290     double AionicRadius(int k = 0) const;
02291 
02292     /**
02293      *
02294      * formPitzer():
02295      *
02296      *  Returns the form of the Pitzer parameterization used
02297      */
02298     int formPitzer() const { return m_formPitzer; }
02299 
02300     /**
02301      * Print out all of the input coefficients.
02302      */
02303     void printCoeffs () const;
02304 
02305     //!  Get the array of unscaled non-dimensional molality based 
02306     //!  activity coefficients at the current solution temperature, 
02307     //!  pressure, and solution concentration.
02308     /*!
02309      *  See Denbigh p. 278 for a thorough discussion. This class must be overwritten in
02310      *  classes which derive from %MolalityVPSSTP. This function takes over from the
02311      *  molar-based activity coefficient calculation, getActivityCoefficients(), in
02312      *  derived classes.
02313      *
02314      * @param acMolality Output vector containing the molality based activity coefficients.
02315      *                   length: m_kk.
02316      */
02317     void getUnscaledMolalityActivityCoefficients(doublereal *acMolality) const;
02318 
02319   private:
02320 
02321     //! Apply the current phScale to a set of activity Coefficients
02322     /*!
02323      *  See the Eq3/6 Manual for a thorough discussion.
02324      *
02325      */
02326     void s_updateScaling_pHScaling() const;
02327 
02328     //!  Apply the current phScale to a set of derivatives of the activity Coefficients
02329     //!  wrt temperature
02330     /*! 
02331      *  See the Eq3/6 Manual for a thorough discussion of the need
02332      */
02333     void s_updateScaling_pHScaling_dT() const;
02334 
02335     //!  Apply the current phScale to a set of 2nd derivatives of the activity Coefficients
02336     //!  wrt temperature
02337     /*! 
02338      *  See the Eq3/6 Manual for a thorough discussion of the need
02339      */
02340     void s_updateScaling_pHScaling_dT2() const;
02341 
02342     //!  Apply the current phScale to a set of derivatives of the activity Coefficients
02343     //!  wrt pressure
02344     /*! 
02345      *  See the Eq3/6 Manual for a thorough discussion of the need
02346      */
02347     void s_updateScaling_pHScaling_dP() const;
02348 
02349  
02350     //!  Calculate the Chlorine activity coefficient on the NBS scale
02351     /*!
02352      *  We assume here that the m_IionicMolality variable is up to date.
02353      */
02354     doublereal s_NBS_CLM_lnMolalityActCoeff() const;
02355 
02356     //!  Calculate the temperature derivative of the Chlorine activity coefficient
02357     //!  on the NBS scale
02358     /*! 
02359      *  We assume here that the m_IionicMolality variable is up to date.
02360      */
02361     doublereal s_NBS_CLM_dlnMolalityActCoeff_dT() const;
02362 
02363     //!  Calculate the second temperature derivative of the Chlorine activity coefficient
02364     //!  on the NBS scale
02365     /*! 
02366      *  We assume here that the m_IionicMolality variable is up to date.
02367      */
02368     doublereal s_NBS_CLM_d2lnMolalityActCoeff_dT2() const;
02369     
02370     //!  Calculate the pressure derivative of the Chlorine activity coefficient
02371     /*! 
02372      *  We assume here that the m_IionicMolality variable is up to date.
02373      */
02374     doublereal s_NBS_CLM_dlnMolalityActCoeff_dP() const;
02375 
02376     //@}
02377          
02378   private:
02379 
02380     /**
02381      * This is the form of the Pitzer parameterization
02382      * used in this model.
02383      * The options are described at the top of this document,
02384      * and in the general documentation.
02385      * The list is repeated here:
02386      *
02387      * PITZERFORM_BASE  = 0    (only one supported atm)
02388      *
02389      */
02390     int m_formPitzer;
02391 
02392     /**
02393      * This is the form of the temperature dependence of Pitzer
02394      * parameterization used in the model.
02395      *
02396      *       PITZER_TEMP_CONSTANT   0
02397      *       PITZER_TEMP_LINEAR     1
02398      *       PITZER_TEMP_COMPLEX1   2
02399      */
02400     int m_formPitzerTemp;
02401 
02402     /**
02403      * Format for the generalized concentration:
02404      *
02405      *  0 = unity
02406      *  1 = molar_volume
02407      *  2 = solvent_volume    (default)
02408      *
02409      * The generalized concentrations can have three different forms
02410      * depending on the value of the member attribute m_formGC, which
02411      * is supplied in the constructor.
02412      *                          <TABLE>
02413      *  <TR><TD> m_formGC </TD><TD> GeneralizedConc </TD><TD> StandardConc </TD></TR>
02414      *  <TR><TD> 0        </TD><TD> X_k             </TD><TD> 1.0          </TD></TR>
02415      *  <TR><TD> 1        </TD><TD> X_k / V_k       </TD><TD> 1.0 / V_k    </TD></TR>
02416      *  <TR><TD> 2        </TD><TD> X_k / V_N       </TD><TD> 1.0 / V_N    </TD></TR>
02417      *                         </TABLE>
02418      *
02419      * The value and form of the generalized concentration will affect
02420      * reaction rate constants involving species in this phase.
02421      *
02422      * (HKM Note: Using option #1 may lead to spurious results and
02423      *  has been included only with warnings. The reason is that it
02424      *  molar volumes of electrolytes may often be negative. The
02425      *  molar volume of H+ is defined to be zero too. Either options
02426      *  0 or 2 are the appropriate choice. Option 0 leads to 
02427      *  bulk reaction rate constants which have units of s-1.
02428      *  Option 2 leads to bulk reaction rate constants for 
02429      *  bimolecular rxns which have units of m-3 kmol-1 s-1.)
02430      */
02431     int m_formGC;
02432 
02433     //! Vector containing the electrolyte species type
02434     /*!
02435      * The possible types are:
02436      *  - solvent
02437      *  - Charged Species
02438      *  - weakAcidAssociated
02439      *  - strongAcidAssociated
02440      *  - polarNeutral
02441      *  - nonpolarNeutral  .
02442      */
02443     vector_int  m_electrolyteSpeciesType;
02444 
02445     /**
02446      * Species molar volumes \f$ m^3 kmol^-1 \f$
02447      *  -> m_speciesSize in Constituents.h
02448      */
02449     //array_fp m_speciesMolarVolume;
02450 
02451     /**
02452      *  a_k = Size of the ionic species in the DH formulation
02453      *        units = meters
02454      */
02455     array_fp m_Aionic;
02456 
02457     /**
02458      * Current value of the ionic strength on the molality scale
02459      * Associated Salts, if present in the mechanism,
02460      * don't contribute to the value of the ionic strength
02461      * in this version of the Ionic strength.
02462      */
02463     mutable double m_IionicMolality;
02464 
02465     /**
02466      * Maximum value of the ionic strength allowed in the
02467      * calculation of the activity coefficients.
02468      */
02469     double m_maxIionicStrength;
02470 
02471     /**
02472      * Reference Temperature for the Pitzer formulations.
02473      */
02474     double m_TempPitzerRef;
02475 
02476     /**
02477      * Stoichiometric ionic strength on the molality scale.
02478      * This differs from m_IionicMolality in the sense that
02479      * associated salts are treated as unassociated salts,
02480      * when calculating the Ionic strength by this method.
02481      */
02482     mutable double m_IionicMolalityStoich;
02483 
02484   public:
02485     /**
02486      * Form of the constant outside the Debye-Huckel term
02487      * called A. It's normally a function of temperature 
02488      * and pressure. However, it can be set from the
02489      * input file in order to aid in numerical comparisons.
02490      * Acceptable forms:
02491      *
02492      *       A_DEBYE_CONST  0
02493      *       A_DEBYE_WATER  1
02494      *
02495      * The A_DEBYE_WATER form may be used for water solvents
02496      * with needs to cover varying temperatures and pressures.
02497      * Note, the dielectric constant of water is a relatively
02498      * strong function of T, and its variability must be
02499      * accounted for,
02500      */
02501     int m_form_A_Debye;
02502 
02503   private:
02504     /**
02505      * A_Debye -> this expression appears on the top of the
02506      *            ln actCoeff term in the general Debye-Huckel
02507      *            expression
02508      *            It depends on temperature. And, therefore,
02509      *            most be recalculated whenever T or P changes.
02510      *            This variable is a local copy of the calculation.
02511      *            
02512      *            A_Debye = (F e B_Debye) / (8 Pi epsilon R T)
02513      *
02514      *                 where B_Debye = F / sqrt(epsilon R T/2) 
02515      *                                 (dw/1000)^(1/2)
02516      *
02517      *            A_Debye = (1/ (8 Pi)) (2 Na * dw/1000)^(1/2)
02518      *                       (e * e / (epsilon * kb * T))^(3/2) 
02519      *
02520      *            Units = sqrt(kg/gmol)
02521      *
02522      *            Nominal value = 1.172576 sqrt(kg/gmol)
02523      *                  based on:
02524      *                    epsilon/epsilon_0 = 78.54
02525      *                           (water at 25C)
02526      *                    epsilon_0 = 8.854187817E-12 C2 N-1 m-2
02527      *                    e = 1.60217653 E-19 C
02528      *                    F = 9.6485309E7 C kmol-1
02529      *                    R = 8.314472E3 kg m2 s-2 kmol-1 K-1
02530      *                    T = 298.15 K
02531      *                    B_Debye = 3.28640E9 sqrt(kg/gmol)/m
02532      *                    dw = C_0 * M_0 (density of water) (kg/m3)
02533      *                       = 1.0E3 at 25C
02534      */
02535     mutable double m_A_Debye;
02536 
02537     //!  Water standard state calculator
02538     /*!
02539      *  derived from the equation of state for water. 
02540      */
02541     PDSS *m_waterSS;
02542 
02543     //! density of standard-state water
02544     /*!
02545      * internal temporary variable
02546      */
02547     double m_densWaterSS;
02548 
02549     /**
02550      *  Pointer to the water property calculator
02551      */
02552     WaterProps *m_waterProps;
02553 
02554     /**
02555      * Vector containing the species reference exp(-G/RT) functions
02556      * at T = m_tlast
02557      */
02558     mutable vector_fp      m_expg0_RT;
02559 
02560     /**
02561      * Vector of potential energies for the species.
02562      */
02563     mutable vector_fp      m_pe;
02564 
02565     /**
02566      * Temporary array used in equilibrium calculations
02567      */
02568     mutable vector_fp      m_pp;
02569 
02570     /**
02571      * vector of size m_kk, used as a temporary holding area.
02572      */
02573     mutable vector_fp      m_tmpV;
02574 
02575     /**
02576      * Stoichiometric species charge -> This is for calculations
02577      * of the ionic strength which ignore ion-ion pairing into
02578      * neutral molecules. The Stoichiometric species charge is the
02579      * charge of one of the ion that would occur if the species broke
02580      * into two charged ion pairs.
02581      *  NaCl ->   m_speciesCharge_Stoich = -1;
02582      *  HSO4- -> H+ + SO42-              = -2
02583      *      -> The other charge is calculated.
02584      * For species that aren't ion pairs, its equal to the
02585      * m_speciesCharge[] value.
02586      */
02587     vector_fp  m_speciesCharge_Stoich;
02588 
02589     /**
02590      *  Array of 2D data used in the Pitzer/HMW formulation.
02591      *  Beta0_ij[i][j] is the value of the Beta0 coefficient
02592      *  for the ij salt. It will be nonzero iff i and j are
02593      *  both charged and have opposite sign. The array is also
02594      *  symmetric.
02595      *     counterIJ where counterIJ = m_counterIJ[i][j]
02596      *  is used to access this array.
02597      */
02598     mutable vector_fp  m_Beta0MX_ij;
02599 
02600     //! Derivative of Beta0_ij[i][j] wrt T
02601     /*!
02602      *  vector index is counterIJ
02603      */
02604     mutable vector_fp  m_Beta0MX_ij_L;
02605 
02606     //! Derivative of Beta0_ij[i][j] wrt TT
02607     /*!
02608      *  vector index is counterIJ
02609      */
02610     mutable vector_fp  m_Beta0MX_ij_LL;
02611 
02612     //! Derivative of Beta0_ij[i][j] wrt P
02613     /*!
02614      *  vector index is counterIJ
02615      */
02616     mutable vector_fp  m_Beta0MX_ij_P;
02617 
02618     //! Array of coefficients for Beta0, a variable in Pitzer's papers
02619     /*!
02620      *  column index is counterIJ
02621      *  m_Beta0MX_ij_coeff.ptrColumn(counterIJ) is a double* containing
02622      *  the vector of coefficients for the counterIJ interaction.
02623      */
02624     mutable Array2D    m_Beta0MX_ij_coeff;
02625 
02626     /*!
02627      *  Array of 2D data used in the Pitzer/HMW formulation.
02628      *  Beta1_ij[i][j] is the value of the Beta1 coefficient
02629      *  for the ij salt. It will be nonzero iff i and j are
02630      *  both charged and have opposite sign. The array is also
02631      *  symmetric. 
02632      *     counterIJ where counterIJ = m_counterIJ[i][j]
02633      *     is used to access this array.
02634      */
02635     mutable vector_fp m_Beta1MX_ij;
02636 
02637     //! Derivative of Beta1_ij[i][j] wrt T
02638     /*!
02639      *  vector index is counterIJ
02640      */
02641     mutable vector_fp m_Beta1MX_ij_L;
02642 
02643     //! Derivative of Beta1_ij[i][j] wrt TT
02644     /*!
02645      *  vector index is counterIJ
02646      */
02647     mutable vector_fp m_Beta1MX_ij_LL;
02648 
02649     //! Derivative of Beta1_ij[i][j] wrt P
02650     /*!
02651      *  vector index is counterIJ
02652      */
02653     mutable vector_fp m_Beta1MX_ij_P;
02654 
02655     //! Array of coefficients for Beta1, a variable in Pitzer's papers
02656     /*!
02657      * column index is counterIJ
02658      *  m_Beta1MX_ij_coeff.ptrColumn(counterIJ) is a double* containing
02659      *  the vector of coefficients for the counterIJ interaction.
02660      */
02661     mutable Array2D   m_Beta1MX_ij_coeff;
02662 
02663     /**
02664      *  Array of 2D data used in the Pitzer/HMW formulation.
02665      *  Beta2_ij[i][j] is the value of the Beta2 coefficient
02666      *  for the ij salt. It will be nonzero iff i and j are
02667      *  both charged and have opposite sign, and i and j
02668      *  both have charges of 2 or more. The array is also
02669      *  symmetric.
02670      *  counterIJ where counterIJ = m_counterIJ[i][j]
02671      *  is used to access this array.
02672      */
02673     mutable vector_fp m_Beta2MX_ij;
02674 
02675     //! Derivative of Beta2_ij[i][j] wrt T
02676     /*!
02677      *  vector index is counterIJ
02678      */
02679     mutable vector_fp m_Beta2MX_ij_L;
02680 
02681     //! Derivative of Beta2_ij[i][j] wrt TT
02682     /*!
02683      *  vector index is counterIJ
02684      */
02685     mutable vector_fp m_Beta2MX_ij_LL;
02686 
02687     //! Derivative of Beta2_ij[i][j] wrt P
02688     /*!
02689      *  vector index is counterIJ
02690      */
02691     mutable vector_fp m_Beta2MX_ij_P;
02692 
02693     //! Array of coefficients for Beta2, a variable in Pitzer's papers
02694     /*!
02695      * column index is counterIJ
02696      *  m_Beta2MX_ij_coeff.ptrColumn(counterIJ) is a double* containing
02697      *  the vector of coefficients for the counterIJ interaction.
02698      *  This was added for the YMP database version of the code since it
02699      *  contains temperature-dependent parameters for some 2-2 electrolytes.
02700      */
02701     mutable Array2D   m_Beta2MX_ij_coeff;
02702 
02703     /**
02704      *  Array of 2D data used in the Pitzer/HMW formulation.
02705      *  Alpha1MX_ij[i][j] is the value of the alpha1 coefficient
02706      *  for the ij interaction. It will be nonzero iff i and j are
02707      *  both charged and have opposite sign.
02708      *  It is symmetric wrt i, j.
02709      *  counterIJ where counterIJ = m_counterIJ[i][j]
02710      *  is used to access this array.
02711      */
02712     vector_fp m_Alpha1MX_ij;
02713 
02714    /**
02715      *  Array of 2D data used in the Pitzer/HMW formulation.
02716      *  Alpha2MX_ij[i][j] is the value of the alpha2 coefficient
02717      *  for the ij interaction. It will be nonzero iff i and j are
02718      *  both charged and have opposite sign, and i and j
02719      *  both have charges of 2 or more, usually.
02720      *  It is symmetric wrt i, j.
02721      *  counterIJ, where counterIJ = m_counterIJ[i][j],
02722      *  is used to access this array.
02723      */
02724     vector_fp m_Alpha2MX_ij;
02725 
02726     /**
02727      *  Array of 2D data used in the Pitzer/HMW formulation.
02728      *  CphiMX_ij[i][j] is the value of the Cphi coefficient
02729      *  for the ij interaction. It will be nonzero iff i and j are
02730      *  both charged and have opposite sign, and i and j
02731      *  both have charges of 2 or more. The array is also
02732      *  symmetric.
02733      *  counterIJ where counterIJ = m_counterIJ[i][j]
02734      *  is used to access this array.
02735      */
02736     mutable vector_fp m_CphiMX_ij;
02737 
02738     //! Derivative of Cphi_ij[i][j] wrt T
02739     /*!
02740      *  vector index is counterIJ
02741      */
02742     mutable vector_fp m_CphiMX_ij_L;
02743 
02744     //! Derivative of Cphi_ij[i][j] wrt TT
02745     /*!
02746      *  vector index is counterIJ
02747      */
02748     mutable vector_fp m_CphiMX_ij_LL;
02749 
02750     //! Derivative of Cphi_ij[i][j] wrt P
02751     /*!
02752      *  vector index is counterIJ
02753      */
02754     mutable vector_fp m_CphiMX_ij_P;
02755 
02756     //! Array of coefficients for CphiMX, a parameter in the activity
02757     //! coefficient formulation
02758     /*!
02759      *  Column index is counterIJ
02760      *  m_CphiMX_ij_coeff.ptrColumn(counterIJ) is a double* containing
02761      *  the vector of coefficients for the counterIJ interaction.
02762      */
02763     mutable Array2D m_CphiMX_ij_coeff;
02764 
02765     //! Array of 2D data for Theta_ij[i][j] in the Pitzer/HMW formulation.
02766     /*!
02767      *  Array of 2D data used in the Pitzer/HMW formulation.
02768      *  Theta_ij[i][j] is the value of the theta coefficient
02769      *  for the ij interaction. It will be nonzero for charged
02770      *  ions with the same sign. It is symmetric.
02771      *  counterIJ where counterIJ = m_counterIJ[i][j]
02772      *  is used to access this array.
02773      *
02774      *  HKM Recent Pitzer papers have used a functional form
02775      *      for Theta_ij, which depends on the ionic strength.
02776      */
02777     mutable vector_fp m_Theta_ij;
02778 
02779     //! Derivative of Theta_ij[i][j] wrt T
02780     /*!
02781      *  vector index is counterIJ
02782      */
02783     mutable vector_fp m_Theta_ij_L;
02784 
02785     //! Derivative of Theta_ij[i][j] wrt TT
02786     /*!
02787      *  vector index is counterIJ
02788      */
02789     mutable vector_fp m_Theta_ij_LL;
02790 
02791     //! Derivative of Theta_ij[i][j] wrt P
02792     /*!
02793      *  vector index is counterIJ
02794      */
02795     mutable vector_fp m_Theta_ij_P;
02796 
02797     //! Array of coefficients for Theta_ij[i][j] in the Pitzer/HMW formulation.
02798     /*!
02799      *  Theta_ij[i][j] is the value of the theta coefficient
02800      *  for the ij interaction. It will be nonzero for charged
02801      *  ions with the same sign. It is symmetric.
02802      *  Column index is counterIJ.
02803      *  counterIJ where counterIJ = m_counterIJ[i][j]
02804      *  is used to access this array.
02805      *
02806      *  m_Theta_ij_coeff.ptrColumn(counterIJ) is a double* containing
02807      *  the vector of coefficients for the counterIJ interaction.
02808      */
02809     Array2D m_Theta_ij_coeff;
02810 
02811     //! Array of 3D data used in the Pitzer/HMW formulation.
02812     /*!
02813      * Psi_ijk[n] is the value of the psi coefficient for the
02814      * ijk interaction where
02815      *
02816      *   n = k + j * m_kk + i * m_kk * m_kk;
02817      *
02818      * It is potentially nonzero everywhere. 
02819      * The first two coordinates are symmetric wrt cations,
02820      * and the last two coordinates are symmetric wrt anions.
02821      */
02822     mutable vector_fp m_Psi_ijk;
02823 
02824     //! Derivitive of Psi_ijk[n] wrt T
02825     /*! 
02826      *  see m_Psi_ijk for reference on the indexing into this variable. 
02827      */
02828     mutable vector_fp m_Psi_ijk_L;
02829 
02830     //! Derivitive of Psi_ijk[n] wrt TT
02831     /*! 
02832      *  see m_Psi_ijk for reference on the indexing into this variable. 
02833      */
02834     mutable vector_fp m_Psi_ijk_LL;
02835 
02836     //! Derivitive of Psi_ijk[n] wrt P
02837     /*! 
02838      *  see m_Psi_ijk for reference on the indexing into this variable. 
02839      */
02840     mutable vector_fp m_Psi_ijk_P;
02841 
02842     //! Array of coefficients for Psi_ijk[n] in the Pitzer/HMW formulation.
02843     /*!
02844      * Psi_ijk[n] is the value of the psi coefficient for the
02845      * ijk interaction where
02846      *
02847      *   n = k + j * m_kk + i * m_kk * m_kk;
02848      *
02849      * It is potentially nonzero everywhere. 
02850      * The first two coordinates are symmetric wrt cations,
02851      * and the last two coordinates are symmetric wrt anions.
02852      *
02853      *
02854      *  m_Psi_ijk_coeff.ptrColumn(n) is a double* containing
02855      *  the vector of coefficients for the n interaction.
02856      */
02857     Array2D m_Psi_ijk_coeff;
02858 
02859     //! Lambda coefficient for the ij interaction
02860     /*!
02861      * Array of 2D data used in the Pitzer/HMW formulation.
02862      * Lambda_nj[n][j] represents the lambda coefficient for the
02863      * ij interaction. This is a general interaction representing
02864      * neutral species. The neutral species occupy the first
02865      * index, i.e., n. The charged species occupy the j coordinate.
02866      * neutral, neutral interactions are also included here.
02867      */
02868     mutable Array2D   m_Lambda_nj;
02869 
02870     //! Derivative of  Lambda_nj[i][j] wrt T. see m_Lambda_ij
02871     mutable Array2D   m_Lambda_nj_L;
02872 
02873     //! Derivative of  Lambda_nj[i][j] wrt TT
02874     mutable Array2D   m_Lambda_nj_LL;
02875 
02876     //! Derivative of  Lambda_nj[i][j] wrt P
02877     mutable Array2D   m_Lambda_nj_P;
02878 
02879     //! Array of coefficients for Lambda_nj[i][j] in the Pitzer/HMW formulation.
02880     /*!
02881      *  Lambda_ij[i][j] is the value of the theta coefficient
02882      *  for the ij interaction.
02883      *  Array of 2D data used in the Pitzer/HMW formulation.
02884      *  Lambda_ij[i][j] represents the lambda coefficient for the
02885      *  ij interaction. This is a general interaction representing
02886      *  neutral species. The neutral species occupy the first
02887      *  index, i.e., i. The charged species occupy the j coordinate.
02888      *  Neutral, neutral interactions are also included here.
02889      *
02890      *      n = j + m_kk * i
02891      *
02892      *  m_Lambda_ij_coeff.ptrColumn(n) is a double* containing
02893      *  the vector of coefficients for the (i,j) interaction.
02894      */
02895     Array2D  m_Lambda_nj_coeff;
02896 
02897 
02898     //! Mu coefficient for the self-ternary neutral coefficient
02899     /*!
02900      * Array of 2D data used in the Pitzer/HMW formulation.
02901      * Mu_nnn[i] represents the Mu coefficient for the
02902      * nnn interaction. This is a general interaction representing
02903      * neutral species interacting with itself.
02904      */
02905     mutable vector_fp  m_Mu_nnn;
02906 
02907     //! Mu coefficient temperature derivative for the self-ternary neutral coefficient
02908     /*!
02909      * Array of 2D data used in the Pitzer/HMW formulation.
02910      * Mu_nnn_L[i] represents the Mu coefficient temperature derivative for the
02911      * nnn interaction. This is a general interaction representing
02912      * neutral species interacting with itself.
02913      */
02914     mutable vector_fp  m_Mu_nnn_L;
02915 
02916     //! Mu coefficient 2nd temperature derivative for the self-ternary neutral coefficient
02917     /*!
02918      * Array of 2D data used in the Pitzer/HMW formulation.
02919      * Mu_nnn_L[i] represents the Mu coefficient 2nd temperature derivative for the
02920      * nnn interaction. This is a general interaction representing
02921      * neutral species interacting with itself.
02922      */
02923     mutable vector_fp  m_Mu_nnn_LL;
02924 
02925     //! Mu coefficient pressure derivative for the self-ternary neutral coefficient
02926     /*!
02927      * Array of 2D data used in the Pitzer/HMW formulation.
02928      * Mu_nnn_L[i] represents the Mu coefficient pressure derivative for the
02929      * nnn interaction. This is a general interaction representing
02930      * neutral species interacting with itself.
02931      */
02932     mutable vector_fp  m_Mu_nnn_P;
02933 
02934     //! Array of coefficients form_Mu_nnn term
02935     /*!
02936      *
02937      */
02938     Array2D   m_Mu_nnn_coeff;
02939     
02940     
02941     //!  Logarithm of the activity coefficients on the molality
02942     //!  scale.
02943     /*!
02944      *       mutable because we change this if the composition
02945      *       or temperature or pressure changes.
02946      *
02947      *  index is the species index
02948      */
02949     mutable vector_fp m_lnActCoeffMolal_Scaled;
02950 
02951    //!  Logarithm of the activity coefficients on the molality
02952     //!  scale.
02953     /*!
02954      *       mutable because we change this if the composition
02955      *       or temperature or pressure changes.
02956      *
02957      *  index is the species index
02958      */
02959     mutable vector_fp m_lnActCoeffMolal_Unscaled;
02960 
02961     //!  Derivative of the Logarithm of the activity coefficients on the molality
02962     //!  scale wrt T
02963     /*!
02964      *  index is the species index
02965      */
02966     mutable vector_fp m_dlnActCoeffMolaldT_Scaled;
02967 
02968     //!  Derivative of the Logarithm of the activity coefficients on the molality
02969     //!  scale wrt T
02970     /*!
02971      *  index is the species index
02972      */
02973     mutable vector_fp m_dlnActCoeffMolaldT_Unscaled;
02974 
02975     //!  Derivative of the Logarithm of the activity coefficients on the molality
02976     //!  scale wrt TT
02977     /*!
02978      *  index is the species index
02979      */
02980     mutable vector_fp m_d2lnActCoeffMolaldT2_Scaled;
02981 
02982    //!  Derivative of the Logarithm of the activity coefficients on the molality
02983     //!  scale wrt TT
02984     /*!
02985      *  index is the species index
02986      */
02987     mutable vector_fp m_d2lnActCoeffMolaldT2_Unscaled;
02988 
02989     //!  Derivative of the Logarithm of the activity coefficients on the 
02990     //!  molality scale wrt P
02991     /*!
02992      *  index is the species index
02993      */
02994     mutable vector_fp m_dlnActCoeffMolaldP_Scaled;
02995 
02996     //!  Derivative of the Logarithm of the activity coefficients on the 
02997     //!  molality scale wrt P
02998     /*!
02999      *  index is the species index
03000      */
03001     mutable vector_fp m_dlnActCoeffMolaldP_Unscaled;
03002 
03003     /*
03004      * -------- Temporary Variables Used in the Activity Coeff Calc
03005      */
03006 
03007     //! Cropped and modified values of the molalities used in activity
03008     //! coefficient calculations
03009     mutable vector_fp m_molalitiesCropped;
03010 
03011     //! Boolean indicating whether the molalities are cropped or are modified
03012     mutable bool m_molalitiesAreCropped;
03013     
03014     //! a counter variable for keeping track of symmetric binary
03015     //! interactions amongst the solute species.
03016     /*!
03017      * n = m_kk*i + j
03018      * m_CounterIJ[n] = counterIJ
03019      */
03020     mutable array_int m_CounterIJ;
03021 
03022     /**
03023      *  This is elambda, MEC
03024      */
03025     mutable double elambda[17];
03026 
03027     /**
03028      *  This is elambda1, MEC
03029      */
03030     mutable double elambda1[17];
03031 
03032     /**
03033      *  Various temporary arrays used in the calculation of
03034      *  the Pitzer activity coefficents.
03035      *  The subscript, L, denotes the same quantity's derivative
03036      *  wrt temperature
03037      */
03038 
03039     //! This is the value of g(x) in Pitzer's papers
03040     /*!
03041      *  vector index is counterIJ
03042      */
03043     mutable vector_fp m_gfunc_IJ;
03044 
03045     //! This is the value of g2(x2) in Pitzer's papers
03046     /*!
03047      *  vector index is counterIJ
03048      */
03049     mutable vector_fp m_g2func_IJ;
03050 
03051     //! hfunc, was called gprime in Pitzer's paper. However,
03052     //! it's not the derivative of gfunc(x), so I renamed it.
03053     /*!
03054      *  vector index is counterIJ
03055      */
03056     mutable vector_fp m_hfunc_IJ;
03057 
03058     //! hfunc2, was called gprime in Pitzer's paper. However,
03059     //! it's not the derivative of gfunc(x), so I renamed it.
03060     /*!
03061      *  vector index is counterIJ
03062      */
03063     mutable vector_fp m_h2func_IJ;
03064 
03065     //! Intermediate variable called BMX in Pitzer's paper
03066     //! This is the basic cation - anion interaction
03067     /*!
03068      *  vector index is counterIJ
03069      */
03070     mutable vector_fp m_BMX_IJ;
03071 
03072     //! Derivative of BMX_IJ wrt T
03073     /*!
03074      *  vector index is counterIJ
03075      */
03076     mutable vector_fp m_BMX_IJ_L;
03077 
03078     //! Derivative of BMX_IJ wrt TT
03079     /*!
03080      *  vector index is counterIJ
03081      */
03082     mutable vector_fp m_BMX_IJ_LL;
03083 
03084     //! Derivative of BMX_IJ wrt P
03085     /*!
03086      *  vector index is counterIJ
03087      */
03088     mutable vector_fp m_BMX_IJ_P;
03089 
03090     //! Intermediate variable called BprimeMX in Pitzer's paper
03091     /*!
03092      *  vector index is counterIJ
03093      */
03094     mutable vector_fp m_BprimeMX_IJ;
03095 
03096     //! Derivative of BprimeMX wrt T 
03097     /*!
03098      *  vector index is counterIJ
03099      */
03100     mutable vector_fp m_BprimeMX_IJ_L;
03101 
03102     //! Derivative of BprimeMX wrt TT
03103     /*!
03104      *  vector index is counterIJ
03105      */
03106     mutable vector_fp m_BprimeMX_IJ_LL;
03107 
03108     //! Derivative of BprimeMX wrt P
03109     /*!
03110      *  vector index is counterIJ
03111      */
03112     mutable vector_fp m_BprimeMX_IJ_P;
03113 
03114     //! Intermediate variable called BphiMX in Pitzer's paper
03115     /*!
03116      *  vector index is counterIJ
03117      */
03118     mutable vector_fp m_BphiMX_IJ;
03119 
03120     //! Derivative of BphiMX_IJ wrt T
03121     /*!
03122      *  vector index is counterIJ
03123      */
03124     mutable vector_fp m_BphiMX_IJ_L;
03125 
03126     //! Derivative of BphiMX_IJ wrt TT
03127     /*!
03128      *  vector index is counterIJ
03129      */
03130     mutable vector_fp m_BphiMX_IJ_LL;
03131 
03132     //! Derivative of BphiMX_IJ wrt P
03133     /*!
03134      *  vector index is counterIJ
03135      */
03136     mutable vector_fp m_BphiMX_IJ_P;
03137 
03138     //! Intermediate variable called Phi in Pitzer's paper
03139     /*!
03140      *  vector index is counterIJ
03141      */
03142     mutable vector_fp m_Phi_IJ;
03143 
03144     //! Derivative of m_Phi_IJ wrt T
03145     /*!
03146      *  vector index is counterIJ
03147      */
03148     mutable vector_fp m_Phi_IJ_L;
03149 
03150     //! Derivative of m_Phi_IJ wrt TT
03151     /*!
03152      *  vector index is counterIJ
03153      */
03154     mutable vector_fp m_Phi_IJ_LL;
03155 
03156     //! Derivative of m_Phi_IJ wrt P
03157     /*!
03158      *  vector index is counterIJ
03159      */
03160     mutable vector_fp m_Phi_IJ_P;
03161 
03162     //! Intermediate variable called Phiprime in Pitzer's paper
03163     /*!
03164      *  vector index is counterIJ
03165      */
03166     mutable vector_fp m_Phiprime_IJ;
03167 
03168     //! Intermediate variable called PhiPhi in Pitzer's paper
03169     /*!
03170      *  vector index is counterIJ
03171      */
03172     mutable vector_fp m_PhiPhi_IJ;
03173 
03174     //! Derivative of m_PhiPhi_IJ wrt T
03175     /*!
03176      *  vector index is counterIJ
03177      */
03178     mutable vector_fp m_PhiPhi_IJ_L;
03179 
03180     //! Derivative of m_PhiPhi_IJ wrt TT
03181     /*!
03182      *  vector index is counterIJ
03183      */
03184     mutable vector_fp m_PhiPhi_IJ_LL;
03185 
03186     //! Derivative of m_PhiPhi_IJ wrt P
03187     /*!
03188      *  vector index is counterIJ
03189      */
03190     mutable vector_fp m_PhiPhi_IJ_P;
03191 
03192     //! Intermediate variable called CMX in Pitzer's paper
03193     /*!
03194      *  vector index is counterIJ
03195      */
03196     mutable vector_fp m_CMX_IJ;
03197 
03198     //! Derivative of m_CMX_IJ wrt T
03199     /*!
03200      *  vector index is counterIJ
03201      */
03202     mutable vector_fp m_CMX_IJ_L;
03203 
03204     //! Derivative of m_CMX_IJ wrt TT
03205     /*!
03206      *  vector index is counterIJ
03207      */
03208     mutable vector_fp m_CMX_IJ_LL;
03209 
03210     //! Derivative of m_CMX_IJ wrt P
03211     /*!
03212      *  vector index is counterIJ
03213      */
03214     mutable vector_fp m_CMX_IJ_P;
03215 
03216     //! Intermediate storage of the activity coefficient itself
03217     /*!
03218      *  vector index is the species index
03219      */
03220     mutable vector_fp m_gamma_tmp;
03221 
03222     //! Logarithm of the molal activity coefficients
03223     /*!
03224      *   Normally these are all one. However, stability schemes will change that
03225      */
03226     mutable vector_fp      IMS_lnActCoeffMolal_;
03227 
03228     //! IMS Cutoff type
03229     int IMS_typeCutoff_;
03230   
03231     //! value of the solute mole fraction that centers the cutoff polynomials
03232     //! for the cutoff =1 process;
03233     doublereal IMS_X_o_cutoff_;
03234 
03235     //! gamma_o value for the cutoff process at the zero solvent point
03236     doublereal IMS_gamma_o_min_;
03237 
03238     //! gamma_k minimun for the cutoff process at the zero solvent point
03239     doublereal IMS_gamma_k_min_;
03240 
03241     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03242     doublereal IMS_cCut_;
03243 
03244     //! Parameter in the polyExp cutoff treatment
03245     /*!
03246      *  This is the slope of the f function at the zero solvent point
03247      *  Default value is 0.6
03248      */
03249     doublereal IMS_slopefCut_;
03250     
03251     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03252     doublereal IMS_dfCut_;
03253 
03254     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03255     doublereal IMS_efCut_;
03256 
03257     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03258     doublereal IMS_afCut_;
03259 
03260     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03261     doublereal IMS_bfCut_;
03262 
03263     //! Parameter in the polyExp cutoff treatment
03264     /*!
03265      *  This is the slope of the g function at the zero solvent point
03266      *  Default value is 0.0
03267      */
03268     doublereal IMS_slopegCut_;
03269 
03270     //! Parameter in the polyExp cutoff treatment having to do with 
03271     //! rate of exp decay
03272     doublereal IMS_dgCut_;
03273 
03274     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03275     doublereal IMS_egCut_;
03276 
03277     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03278     doublereal IMS_agCut_;
03279 
03280     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
03281     doublereal IMS_bgCut_;
03282 
03283     //! value of the solvent mole fraction that centers the cutoff polynomials
03284     //! for the cutoff =1 process;
03285     doublereal MC_X_o_cutoff_;
03286 
03287     //! gamma_o value for the cutoff process at the zero solvent point
03288     doublereal MC_X_o_min_;
03289 
03290     //! Parameter in the Molality Exp cutoff treatment
03291     /*!
03292      *  This is the slope of the p function at the zero solvent point
03293      *  Default value is 0.0
03294      */
03295     doublereal MC_slopepCut_;
03296 
03297     //! Parameter in the Molality Exp cutoff treatment
03298     doublereal MC_dpCut_;
03299 
03300     //!  Parameter in the Molality Exp cutoff treatment
03301     doublereal MC_epCut_;
03302 
03303     //! Parameter in the Molality Exp cutoff treatment 
03304     doublereal MC_apCut_;
03305 
03306     //! Parameter in the Molality Exp cutoff treatment 
03307     doublereal MC_bpCut_;
03308 
03309     //! Parameter in the Molality Exp cutoff treatment
03310     doublereal MC_cpCut_;
03311 
03312     //! Parameter in the Molality Exp cutoff treatment
03313     doublereal CROP_ln_gamma_o_min;
03314 
03315     //! Parameter in the Molality Exp cutoff treatment
03316     doublereal CROP_ln_gamma_o_max;
03317 
03318     //! Parameter in the Molality Exp cutoff treatment
03319     doublereal CROP_ln_gamma_k_min;
03320 
03321     //! Parameter in the Molality Exp cutoff treatment
03322     doublereal CROP_ln_gamma_k_max;
03323 
03324     //! This is a boolean-type vector indicating whether
03325     //! a species's activity coefficient is in the cropped regime
03326     /*!
03327      *
03328      *  0 = Not in cropped regime
03329      *  1 = In a transition regime where it is altered but there
03330      *      still may be a temperature or pressure dependence
03331      *  2 = In a cropped regime where there is no temperature 
03332      *      or pressure dependence
03333      */
03334     mutable std::vector<int> CROP_speciesCropped_;
03335 
03336 
03337     //! Local error routine
03338     /*!
03339      * @param msg print out a message and error exit
03340      */
03341     doublereal err(std::string msg) const;
03342 
03343     //!  Initialize all of the species - dependent lengths in the object
03344     void initLengths();
03345 
03346     //! Apply the current phScale to a set of activity Coefficients or 
03347     //! activities
03348     /*!
03349      *  See the Eq3/6 Manual for a thorough discussion.
03350      *
03351      * @param acMolality input/Output vector containing the molality based 
03352      *                   activity coefficients. length: m_kk.
03353      */
03354     virtual void applyphScale(doublereal *acMolality) const;
03355 
03356   private:
03357     /*
03358      * This function will be called to update the internally storred
03359      * natural logarithm of the molality activity coefficients 
03360      */
03361     void s_update_lnMolalityActCoeff() const;
03362 
03363     //! This function calculates the temperature derivative of the
03364      //! natural logarithm of the molality activity coefficients.
03365      /*!
03366       * This is the private function. It does all of the direct work.
03367       */
03368     void s_update_dlnMolalityActCoeff_dT() const;
03369 
03370     /**
03371      * This function calculates the temperature second derivative
03372      * of the natural logarithm of the molality activity 
03373      * coefficients.
03374      */
03375     void s_update_d2lnMolalityActCoeff_dT2() const;
03376 
03377     /**
03378      * This function calculates the pressure derivative of the
03379      * natural logarithm of the molality activity coefficients.
03380      */
03381     void s_update_dlnMolalityActCoeff_dP() const;
03382 
03383     //! This function will be called to update the internally storred
03384     //! natural logarithm of the molality activity coefficients
03385     /*
03386      * Normally they are all one. However, sometimes they are not,
03387      * due to stability schemes
03388      *
03389      *    gamma_k_molar =  gamma_k_molal / Xmol_solvent
03390      *
03391      *    gamma_o_molar = gamma_o_molal
03392      */
03393     void s_updateIMS_lnMolalityActCoeff() const;
03394 
03395   private:
03396     /**
03397      * This function does the main pitzer coefficient 
03398      * calculation
03399      */
03400     void s_updatePitzer_lnMolalityActCoeff() const;
03401 
03402     //!  Calculates the temperature derivative of the
03403     //!  natural logarithm of the molality activity coefficients.
03404     /*!
03405      * Public function makes sure that all dependent data is
03406      * up to date, before calling a private function
03407      */
03408     void s_updatePitzer_dlnMolalityActCoeff_dT() const;
03409 
03410     /**
03411      * This function calculates the temperature second derivative
03412      * of the natural logarithm of the molality activity 
03413      * coefficients.
03414      */
03415     void s_updatePitzer_d2lnMolalityActCoeff_dT2() const;
03416 
03417     //!  Calculates the Pressure derivative of the
03418     //!  natural logarithm of the molality activity coefficients.
03419     /*!
03420      * Public function makes sure that all dependent data is
03421      * up to date, before calling a private function
03422      */
03423     void s_updatePitzer_dlnMolalityActCoeff_dP() const;
03424 
03425  
03426    
03427 
03428     //! Calculates the Pitzer coefficients' dependence on the temperature.
03429     /*!
03430      *  It will also calculate the temperature
03431      *  derivatives of the coefficients, as they are important
03432      *  in the calculation of the latent heats and the
03433      *  heat capacities of the mixtures.
03434      *
03435      *   @param doDerivs If >= 1, then the routine will calculate
03436      *                  the first derivative. If >= 2, the 
03437      *                  routine will calculate the first and second
03438      *                  temperature derivative.
03439      *                  default = 2
03440      */
03441     void s_updatePitzer_CoeffWRTemp(int doDerivs = 2) const;
03442 
03443   
03444     
03445     //! Calculate the lambda interactions.
03446     /*!
03447      *
03448      * Calculate E-lambda terms for charge combinations of like sign,
03449      *   using method of Pitzer (1975).
03450      *
03451      * @param is Ionic strength
03452      */
03453     void calc_lambdas(double is) const;
03454 
03455     /**
03456      *  Calculate etheta and etheta_prime
03457      *
03458      *  This interaction will be nonzero for species with the 
03459      *  same charge. this routine is not to be called for 
03460      *  neutral species; it core dumps or error exits.
03461      *
03462      * MEC implementation routine.
03463      *
03464      *  @param z1 charge of the first molecule
03465      *  @param z2 charge of the second molecule
03466      *  @param etheta return pointer containing etheta
03467      *  @param etheta_prime Return pointer containing etheta_prime.
03468      *
03469      *  This routine uses the internal variables, 
03470      *   elambda[] and elambda1[].
03471      *
03472      *  There is no prohibition against calling 
03473      *
03474      */
03475     void calc_thetas(int z1, int z2,
03476                      double *etheta, double *etheta_prime) const;
03477 
03478     //! Set up a counter variable for keeping track of symmetric binary
03479     //! interactactions amongst the solute species.
03480     /*!
03481      * The purpose of this is to squeeze the ij parameters into a
03482      * compressed single counter.
03483      *
03484      * n = m_kk*i + j 
03485      * m_Counter[n] = counter
03486      */
03487     void counterIJ_setup() const;
03488 
03489     //! Calculate the cropped molalities
03490     /*!
03491      * This is an internal routine that calculates values
03492      * of m_molalitiesCropped from m_molalities
03493      */
03494     void calcMolalitiesCropped() const;
03495 
03496     //! Process an XML node called "binarySaltParameters"
03497     /*!
03498      * This node contains all of the parameters necessary to describe
03499      * the Pitzer model for that particular binary salt.
03500      * This function reads the XML file and writes the coefficients
03501      * it finds to an internal data structures.
03502      *
03503      * @param BinSalt  reference to the XML_Node named binarySaltParameters
03504      *                 containing the
03505      *                 anion - cation interaction
03506      */
03507     void readXMLBinarySalt(XML_Node &BinSalt);
03508 
03509     //! Process an XML node called "thetaAnion"
03510     /*!
03511      * This node contains all of the parameters necessary to describe
03512      * the binary interactions between two anions.
03513      *
03514      * @param BinSalt  reference to the XML_Node named thetaAnion
03515      *                 containing the
03516      *                 anion - anion interaction
03517      */
03518     void readXMLThetaAnion(XML_Node &BinSalt);
03519 
03520     //! Process an XML node called "thetaCation"
03521     /*!
03522      * This node contains all of the parameters necessary to describe
03523      * the binary interactions between two cations.
03524      *
03525      * @param BinSalt  reference to the XML_Node named thetaCation
03526      *                 containing the
03527      *                 cation - cation interaction
03528      */
03529     void readXMLThetaCation(XML_Node &BinSalt);
03530 
03531     //! Process an XML node called "psiCommonAnion"
03532     /*!
03533      * This node contains all of the parameters necessary to describe
03534      * the ternary interactions between one anion and two cations.
03535      *
03536      * @param BinSalt  reference to the XML_Node named psiCommonAnion
03537      *                 containing the
03538      *                 anion - cation1 - cation2 interaction
03539      */
03540     void readXMLPsiCommonAnion(XML_Node &BinSalt);
03541 
03542     //! Process an XML node called "psiCommonCation"
03543     /*!
03544      * This node contains all of the parameters necessary to describe
03545      * the ternary interactions between one cation and two anions.
03546      *
03547      * @param BinSalt  reference to the XML_Node named psiCommonCation
03548      *                 containing the
03549      *                 cation - anion1 - anion2 interaction
03550      */
03551     void readXMLPsiCommonCation(XML_Node &BinSalt);
03552 
03553     //! Process an XML node called "lambdaNeutral"
03554     /*!
03555      * This node contains all of the parameters necessary to describe
03556      * the binary interactions between one neutral species and
03557      * any other species (neutral or otherwise) in the mechanism.
03558      *
03559      * @param BinSalt  reference to the XML_Node named lambdaNeutral
03560      *                 containing multiple
03561      *                 Neutral - species interactions
03562      */
03563     void readXMLLambdaNeutral(XML_Node &BinSalt);
03564 
03565     //! Process an XML node called "MunnnNeutral"
03566     /*!
03567      * This node contains all of the parameters necessary to describe
03568      * the self-ternary interactions for one neutral species.
03569      *
03570      * @param BinSalt  reference to the XML_Node named Munnn
03571      *                 containing the self-ternary interaction
03572      */
03573     void readXMLMunnnNeutral(XML_Node &BinSalt);
03574 
03575     //! Process an XML node called "zetaCation"
03576     /*!
03577      * This node contains all of the parameters necessary to describe
03578      * the ternary interactions between one neutral, one cation, and one anion.
03579      *
03580      * @param BinSalt  reference to the XML_Node named psiCommonCation
03581      *                 containing the
03582      *                 neutral - cation - anion interaction
03583      */
03584     void readXMLZetaCation(const XML_Node &BinSalt);
03585 
03586     //! Process an XML node called "croppingCoefficients"
03587     //! for the cropping coefficients values
03588     /*!
03589      * @param acNode Activity Coefficient XML Node
03590      */
03591     void readXMLCroppingCoefficients(const XML_Node &acNode);
03592 
03593     //! Precalculate the IMS Cutoff parameters for typeCutoff = 2
03594     void  calcIMSCutoffParams_();
03595 
03596     //! Calculate molality cut-off parameters
03597     void  calcMCCutoffParams_();
03598 
03599     //! Utility function to assign an integer value from a string
03600     //! for the ElectrolyteSpeciesType field.
03601     /*!
03602      *  @param estString string name of the electrolyte species type
03603      */
03604     static int interp_est(std::string estString);
03605 
03606   public:
03607     /*!
03608      * Turn on copious debug printing when this
03609      * is true and DEBUG_MODE is turned on.
03610      */
03611     mutable int m_debugCalc;
03612 
03613     //! Return int specifying the amount of debug printing
03614     /*!
03615      *  This will return 0 if DEBUG_MODE is not turned on
03616      */
03617     int debugPrinting();
03618   };
03619 
03620 }
03621         
03622 #endif
03623 
Generated by  doxygen 1.6.3