From 5de3c69245abcb24ed0e8a34b280974fbd64fb65 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Sat, 22 Oct 2005 00:18:16 +0000
Subject: [PATCH] General commit for a reworking of the Species reference state
 thermo calculation. This is a reclarification of the reference state thermo
 calculations for individual species, and an expansion to handle liquid phase
 thermo needs.

There is now a virtual base class for the calculation of reference state
thermo functions for  individual species. It is called
SpeciesThermoInterpType.

There is also a class which allows for a complete general
calculation of the reference state species thermo for a phase,
GeneralSpeciesThermo.

Some of this new functionality may be relegated to ifdef blocks
in the future to limit the amount of code for users who don't
need the new functionality.
---
 Cantera/src/ConstCpPoly.cpp           | 132 +++++++++
 Cantera/src/ConstCpPoly.h             |  66 +++++
 Cantera/src/GeneralSpeciesThermo.cpp  | 228 ++++++++++++++
 Cantera/src/GeneralSpeciesThermo.h    | 123 ++++++++
 Cantera/src/Makefile.in               |  11 +-
 Cantera/src/Mu0Poly.cpp               | 409 ++++++++++++++++++++++++++
 Cantera/src/Mu0Poly.h                 | 136 +++++++++
 Cantera/src/NasaPoly1.h               |  79 ++++-
 Cantera/src/NasaPoly2.h               | 217 ++++++++++++++
 Cantera/src/NasaThermo.h              |  19 +-
 Cantera/src/ShomatePoly.h             | 256 ++++++++++++++--
 Cantera/src/ShomateThermo.h           |  43 ++-
 Cantera/src/SimpleThermo.h            |   5 +-
 Cantera/src/SpeciesThermo.h           |  25 +-
 Cantera/src/SpeciesThermoFactory.cpp  |  54 +++-
 Cantera/src/SpeciesThermoInterpType.h |  49 +++
 Cantera/src/SpeciesThermoMgr.h        |  19 +-
 Cantera/src/speciesThermoTypes.h      |   8 +
 18 files changed, 1797 insertions(+), 82 deletions(-)
 create mode 100644 Cantera/src/ConstCpPoly.cpp
 create mode 100644 Cantera/src/ConstCpPoly.h
 create mode 100644 Cantera/src/GeneralSpeciesThermo.cpp
 create mode 100644 Cantera/src/GeneralSpeciesThermo.h
 create mode 100644 Cantera/src/Mu0Poly.cpp
 create mode 100644 Cantera/src/Mu0Poly.h
 create mode 100644 Cantera/src/NasaPoly2.h
 create mode 100644 Cantera/src/SpeciesThermoInterpType.h

diff --git a/Cantera/src/ConstCpPoly.cpp b/Cantera/src/ConstCpPoly.cpp
new file mode 100644
index 000000000..d0d998385
--- /dev/null
+++ b/Cantera/src/ConstCpPoly.cpp
@@ -0,0 +1,132 @@
+/**
+ *  @file ConstCpPoly.h
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+// Copyright 2001  California Institute of Technology
+
+
+
+#include "ConstCpPoly.h"
+
+namespace Cantera {
+
+    ConstCpPoly::ConstCpPoly() 
+	:  m_t0(0.0),
+	   m_cp0_R(0.0),
+	   m_h0_R(0.0),
+	   m_s0_R(0.0),
+	   m_logt0(0.0),
+	   m_lowT(0.0),
+	   m_highT(0.0),
+	   m_Pref(0.0),
+	   m_index(0) {
+    }
+    
+    ConstCpPoly::ConstCpPoly(int n, doublereal tlow, doublereal thigh, 
+			     doublereal pref,
+			     const doublereal* coeffs) :
+	m_lowT      (tlow),
+	m_highT     (thigh),
+	m_Pref      (pref),
+	m_index     (n) {
+	m_t0 = coeffs[0];
+	m_h0_R = coeffs[1]  / GasConstant;
+	m_s0_R = coeffs[2]  / GasConstant;
+	m_cp0_R = coeffs[3] / GasConstant;
+	m_logt0 = log(m_t0);
+    }
+
+    ConstCpPoly::ConstCpPoly(const ConstCpPoly& b) :
+	m_t0        (b.m_t0),
+	m_cp0_R     (b.m_cp0_R),
+	m_h0_R      (b.m_h0_R),
+	m_s0_R      (b.m_s0_R),
+	m_logt0     (b.m_logt0),
+	m_lowT      (b.m_lowT),
+	m_highT     (b.m_highT),
+	m_Pref      (b.m_Pref),
+	m_index     (b.m_index)
+    {
+    }    
+
+    ConstCpPoly& ConstCpPoly::operator=(const ConstCpPoly& b) {
+	if (&b != this) {
+	  m_t0    = b.m_t0;
+	  m_cp0_R = b.m_cp0_R;
+	  m_h0_R  = b.m_h0_R;
+	  m_s0_R  = b.m_s0_R;
+	  m_logt0 = b.m_logt0;
+	  m_lowT  = b.m_lowT;
+	  m_highT = b.m_highT;
+	  m_Pref  = b.m_Pref;
+	  m_index = b.m_index;
+	}
+	return *this;
+    }
+    
+    ConstCpPoly::~ConstCpPoly(){}
+
+    SpeciesThermoInterpType *
+    ConstCpPoly::duplMyselfAsSpeciesThermoInterpType() const {
+	ConstCpPoly* newCCP = new ConstCpPoly(*this);
+	return (SpeciesThermoInterpType*) newCCP;
+    }
+
+    doublereal ConstCpPoly::minTemp() const {
+	return m_lowT;
+    }
+    doublereal ConstCpPoly::maxTemp() const { 
+	return m_highT;
+    }
+    doublereal ConstCpPoly::refPressure() const { 
+	return m_Pref;
+    }
+      
+    void ConstCpPoly::updateProperties(const doublereal* tt, 
+				       doublereal* cp_R,
+				       doublereal* h_RT,
+				       doublereal* s_R) const {
+	double t = *tt;
+	doublereal logt = log(t);
+	doublereal rt = 1.0/t;
+	cp_R[m_index] = m_cp0_R;
+	h_RT[m_index] = rt*(m_h0_R + (t - m_t0) * m_cp0_R);
+	s_R[m_index]  = m_s0_R + m_cp0_R * (logt - m_logt0);
+    }
+
+    void ConstCpPoly::updatePropertiesTemp(const doublereal temp, 
+					   doublereal* cp_R,
+					   doublereal* h_RT,
+					   doublereal* s_R) const {
+	doublereal logt = log(temp);
+	doublereal rt = 1.0/temp;
+	cp_R[m_index] = m_cp0_R;
+	h_RT[m_index] = rt*(m_h0_R + (temp - m_t0) * m_cp0_R);
+	s_R[m_index]  = m_s0_R + m_cp0_R * (logt - m_logt0);
+    }
+
+    void ConstCpPoly::reportParameters(int &n, int &type,
+				       doublereal &tlow, doublereal &thigh,
+				       doublereal &pref,
+				       doublereal* const coeffs) const {
+	n = m_index;
+	type = CONSTANT_CP;
+	tlow = m_lowT;
+	thigh = m_highT;
+	pref = m_Pref;
+	coeffs[0] = m_t0;
+	coeffs[1] = m_h0_R * GasConstant;
+	coeffs[2] = m_s0_R * GasConstant;
+	coeffs[3] = m_cp0_R * GasConstant;
+    }
+
+}
+
+
+
+
+
diff --git a/Cantera/src/ConstCpPoly.h b/Cantera/src/ConstCpPoly.h
new file mode 100644
index 000000000..233901465
--- /dev/null
+++ b/Cantera/src/ConstCpPoly.h
@@ -0,0 +1,66 @@
+/**
+ *  @file ConstCpPoly.h
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+// Copyright 2001  California Institute of Technology
+
+
+#ifndef CT_CONSTCPPOLY_H
+#define CT_CONSTCPPOLY_H
+
+#include "SpeciesThermoInterpType.h"
+
+namespace Cantera {
+
+
+    class ConstCpPoly: public SpeciesThermoInterpType {
+
+    public:
+
+        ConstCpPoly();
+        ConstCpPoly(int n, doublereal tlow, doublereal thigh, 
+		    doublereal pref,
+		    const doublereal* coeffs);
+	ConstCpPoly(const ConstCpPoly&);
+	ConstCpPoly& operator=(const ConstCpPoly&);
+        virtual ~ConstCpPoly();
+	virtual SpeciesThermoInterpType *
+    	   duplMyselfAsSpeciesThermoInterpType() const;
+
+        doublereal minTemp() const;
+        doublereal maxTemp() const;
+        doublereal refPressure() const;
+        virtual int reportType() const { return CONSTANT_CP; }
+      
+        void updateProperties(const doublereal* tt, 
+			      doublereal* cp_R, doublereal* h_RT, 
+			      doublereal* s_R) const;
+
+        void updatePropertiesTemp(const doublereal temp, 
+				  doublereal* cp_R, doublereal* h_RT,
+				  doublereal* s_R) const;
+
+	void reportParameters(int &n, int &type,
+			      doublereal &tlow, doublereal &thigh,
+			      doublereal &pref,
+			      doublereal* const coeffs) const;
+    protected:
+	doublereal m_t0;
+        doublereal m_cp0_R;
+	doublereal m_h0_R;
+	doublereal m_s0_R;
+	doublereal m_logt0;
+        doublereal m_lowT, m_highT, m_Pref;
+        int m_index;
+    
+    private:
+
+    };
+
+}
+
+#endif
diff --git a/Cantera/src/GeneralSpeciesThermo.cpp b/Cantera/src/GeneralSpeciesThermo.cpp
new file mode 100644
index 000000000..15ba9da95
--- /dev/null
+++ b/Cantera/src/GeneralSpeciesThermo.cpp
@@ -0,0 +1,228 @@
+/**
+ *  @file GeneralSpeciesThermo.cpp
+ *
+ */
+                
+// Copyright 2001-2004  California Institute of Technology
+                
+#include "GeneralSpeciesThermo.h"
+#include "NasaPoly1.h"
+#include "NasaPoly2.h"
+#include "ShomatePoly.h"
+#include "ConstCpPoly.h"
+#include "Mu0Poly.h"
+#include "SpeciesThermoFactory.h"
+#include <iostream>
+using namespace std;
+                                                           
+namespace Cantera {
+
+
+    /*
+     * Constructors
+     */
+    GeneralSpeciesThermo::GeneralSpeciesThermo() :
+	SpeciesThermo(),
+	m_tlow_max(0.0), 
+	m_thigh_min(1.0E30),
+	m_p0(OneAtm),
+        m_kk(0) 
+    {
+	m_tlow_max = 0.0;
+	m_thigh_min = 1.0E30;
+    }
+ 
+    GeneralSpeciesThermo::
+    GeneralSpeciesThermo(const GeneralSpeciesThermo &b) :
+	m_tlow_max(b.m_tlow_max), 
+	m_thigh_min(b.m_thigh_min),
+	m_kk(b.m_kk) {
+	m_sp = b.m_sp;
+    }
+
+    const GeneralSpeciesThermo& 
+    GeneralSpeciesThermo::operator=(const GeneralSpeciesThermo &b) {
+	if (&b != this) {
+	  m_tlow_max = b.m_tlow_max;
+	  m_thigh_min = b.m_thigh_min;
+	  m_kk = b.m_kk;
+	  m_sp = b.m_sp;
+	}
+	return *this;
+    }
+
+    GeneralSpeciesThermo::~GeneralSpeciesThermo() {
+	for (int k = 0; k < m_kk; k++) {
+	  SpeciesThermoInterpType *sp = m_sp[k];
+	  if (sp) {
+	    delete (sp);
+	    m_sp[k] = 0;
+	  }
+	}
+    }
+
+
+    SpeciesThermo *
+    GeneralSpeciesThermo::duplMyselfAsSpeciesThermo() const {
+	GeneralSpeciesThermo *gsth = new GeneralSpeciesThermo(*this);
+	return (SpeciesThermo *) gsth;
+    }
+
+    
+    /**
+     * Install parameterization for a species.
+     * @param index    Species index
+     * @param type     ignored, since only NASA type is supported
+     * @param c        coefficients. These are
+     * - c[0]          midpoint temperature
+     * - c[1] - c[7]   coefficients for low T range
+     * - c[8] - c[14]  coefficients for high T range
+     */
+    void GeneralSpeciesThermo::install(string name,
+				       int index,
+				       int type, 
+				       const doublereal* c, 
+				       doublereal minTemp,
+				       doublereal maxTemp,
+				       doublereal refPressure) {
+	/*
+	 * Resize the arrays if necessary, filling the empty
+	 * slots with the zero pointer.
+	 */
+	if (index > m_kk - 1) {
+	  m_sp.resize(index+1, 0);
+          m_kk = index+1;
+	}
+
+	/*
+	 * Create the necessary object
+	 */
+	switch (type) {
+	case NASA1:
+	    m_sp[index] = new NasaPoly1(index, minTemp, maxTemp,
+					refPressure, c);
+	    break;
+	case SHOMATE1:
+	    m_sp[index] = new ShomatePoly(index, minTemp, maxTemp,
+					  refPressure, c);
+	    break;
+        case CONSTANT_CP:
+        case SIMPLE:
+	    m_sp[index] = new ConstCpPoly(index, minTemp, maxTemp,
+					  refPressure, c);
+	    break;
+	case MU0_INTERP:
+	    m_sp[index] = new Mu0Poly(index, minTemp, maxTemp,
+				      refPressure, c);
+	    break;
+        case SHOMATE2:
+	    m_sp[index] = new ShomatePoly2(index, minTemp, maxTemp,
+					  refPressure, c);
+	    break;
+	case NASA2:
+	    m_sp[index] = new NasaPoly2(index, minTemp, maxTemp,
+					refPressure, c);
+	    break;
+	default:
+	    throw UnknownSpeciesThermoModel(
+                "GeneralSpeciesThermo::install",
+		"unknown species type", int2str(type));
+            break;
+	}
+	m_tlow_max = max(minTemp, m_tlow_max);
+	m_thigh_min = min(maxTemp, m_thigh_min);
+    }
+
+    /**
+     *  Update the properties for all species;
+     */
+    void GeneralSpeciesThermo::
+    update_one(int k, doublereal t, doublereal* cp_R, 
+	       doublereal* h_RT, doublereal* s_R) const {
+	SpeciesThermoInterpType * sp_ptr = m_sp[k];
+	sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
+    }
+
+    
+    /**
+     *  Update the properties for all species;
+     */
+    void GeneralSpeciesThermo::
+    update(doublereal t, doublereal* cp_R, 
+	   doublereal* h_RT, doublereal* s_R) const {
+	vector<SpeciesThermoInterpType *>::const_iterator _begin, _end;
+	_begin  = m_sp.begin();
+	_end    = m_sp.end();
+	SpeciesThermoInterpType * sp_ptr;
+	for (; _begin != _end; ++_begin) {
+	  sp_ptr = *(_begin);
+	  sp_ptr->updatePropertiesTemp(t, cp_R, h_RT, s_R);
+	}
+    }
+
+    /**
+     * This utility function reports the type of parameterization
+     * used for the species, index.
+     */
+    int GeneralSpeciesThermo::reportType(int index) const {
+	SpeciesThermoInterpType *sp = m_sp[index];
+	return sp->reportType();
+    }
+    
+    /**
+     * This utility function reports back the type of 
+     * parameterization and all of the parameters for the 
+     * species, index.
+     *  For the NASA object, there are 15 coefficients.
+     */
+    void GeneralSpeciesThermo::
+    reportParams(int index, int &type, 
+		 doublereal * const c, 
+		 doublereal &minTemp, 
+		 doublereal &maxTemp, 
+		 doublereal &refPressure) {
+	SpeciesThermoInterpType *sp = m_sp[index];
+	int n;
+	sp->reportParameters(n, type, minTemp, maxTemp, 
+			     refPressure, c);      
+	if (n != index) {
+	  throw CanteraError("  ", "confused");
+	}
+    }
+
+    /**
+     * Return the lowest temperature at which the thermodynamic
+     * parameterization is valid.  If no argument is supplied, the
+     * value is the one for which all species parameterizations
+     * are valid. Otherwise, if an integer argument is given, the
+     * value applies only to the species with that index.
+     */
+    doublereal GeneralSpeciesThermo::minTemp(int k) const {
+	if (k < 0)
+	    return m_tlow_max;
+	else {
+	  SpeciesThermoInterpType *sp = m_sp[k];
+	  return sp->minTemp();	      
+	}
+    }
+
+    doublereal GeneralSpeciesThermo::maxTemp(int k) const {
+	if (k < 0) {
+	  return m_thigh_min;
+	} else {
+	  SpeciesThermoInterpType *sp = m_sp[k];
+	  return sp->maxTemp();
+	}
+    }
+
+    doublereal GeneralSpeciesThermo::refPressure(int k) const {
+	if (k < 0) {
+	  return m_p0;
+	} else {
+	  SpeciesThermoInterpType *sp = m_sp[k];
+	  return sp->refPressure();
+	}
+    }
+    
+
+}
diff --git a/Cantera/src/GeneralSpeciesThermo.h b/Cantera/src/GeneralSpeciesThermo.h
new file mode 100644
index 000000000..8ac2e6d35
--- /dev/null
+++ b/Cantera/src/GeneralSpeciesThermo.h
@@ -0,0 +1,123 @@
+/**
+ * @file GeneralSpeciesThermo.h
+ */
+
+/*
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+
+#ifndef CT_GENERALSPECIESTHERMO_H
+#define CT_GENERALSPECIESTHERMO_H
+#include <string>
+#include "ct_defs.h"
+#include "SpeciesThermoMgr.h"
+#include "NasaPoly1.h"
+#include "speciesThermoTypes.h"
+#include "polyfit.h"
+
+namespace Cantera {
+
+    /**
+     * A species thermodynamic property manager for a phase.
+     * This is a general manager that can handle a wide variety
+     * of species thermodynamic polynomials for individual species.
+     * It is slow, however.
+     *
+     *
+     */
+    class GeneralSpeciesThermo : public SpeciesThermo {
+    
+    public:
+
+	GeneralSpeciesThermo();
+	GeneralSpeciesThermo(const GeneralSpeciesThermo &);
+	const GeneralSpeciesThermo & operator=(const GeneralSpeciesThermo &);
+        virtual ~GeneralSpeciesThermo();
+	virtual SpeciesThermo *duplMyselfAsSpeciesThermo() const ;
+
+        /**
+         * Install parameterization for a species.
+         * @param index    Species index
+         * @param type     ignored, since only NASA type is supported
+         * @param c        coefficients. These are
+         * - c[0]          midpoint temperature
+         * - c[1] - c[7]   coefficients for low T range
+         * - c[8] - c[14]  coefficients for high T range
+         */
+        virtual void install(string name, int index, int type, 
+			     const doublereal* c, 
+			     doublereal minTemp, doublereal maxTemp,
+			     doublereal refPressure);
+        /** 
+         * update the properties for only one species.
+         */
+        virtual void update_one(int k, doublereal t, doublereal* cp_R, 
+				doublereal* h_RT,
+				doublereal* s_R) const;
+	
+        virtual void update(doublereal t, doublereal* cp_R, 
+			    doublereal* h_RT, doublereal* s_R) const;
+                
+        /**
+         * Return the lowest temperature at which the thermodynamic
+         * parameterization is valid.  If no argument is supplied, the
+         * value is the one for which all species parameterizations
+         * are valid. Otherwise, if an integer argument is given, the
+         * value applies only to the species with that index.
+         */
+        virtual doublereal minTemp(int k=-1) const;
+
+        virtual doublereal maxTemp(int k=-1) const;
+
+        virtual doublereal refPressure(int k = -1) const;
+
+	/**
+         * This utility function reports the type of parameterization
+         * used for the species, index.
+         */
+        virtual int reportType(int index) const;
+
+	/**
+	 * This utility function reports back the type of 
+	 * parameterization and all of the parameters for the 
+	 * species, index.
+	 *  For the NASA object, there are 15 coefficients.
+	 */
+	virtual void reportParams(int index, int &type, 
+				  doublereal * const c, 
+				  doublereal &minTemp, 
+				  doublereal &maxTemp,
+				  doublereal &refPressure);
+
+ protected:
+
+	/**
+	 * This is the main unknown in the object. It is 
+	 * a list of pointers to type SpeciesThermoInterpType.
+	 * Note, this object owns the objects, so they are deleted
+	 * in the destructor of this object.
+	 */
+	vector<SpeciesThermoInterpType *>  m_sp;
+        doublereal                         m_tlow_max;
+        doublereal                         m_thigh_min;
+        doublereal                         m_p0;
+
+	/**
+	 * Internal variable indicating the length of the 
+	 * number of species in the phase.
+	 */
+	int m_kk;
+
+    private:
+
+  
+
+    };
+
+}
+
+#endif
+
diff --git a/Cantera/src/Makefile.in b/Cantera/src/Makefile.in
index 7e49be887..d230f0bec 100755
--- a/Cantera/src/Makefile.in
+++ b/Cantera/src/Makefile.in
@@ -39,13 +39,16 @@ BASE       = $(BASE_OBJ)
 
 # thermodynamic properties
 THERMO_OBJ = ThermoPhase.o IdealGasPhase.o ConstDensityThermo.o \
-             SpeciesThermoFactory.o \
+             SpeciesThermoFactory.o ConstCpPoly.o \
+             Mu0Poly.o GeneralSpeciesThermo.o \
              ThermoFactory.o @phase_object_files@
 
 THERMO_H   = ThermoPhase.h IdealGasPhase.h ConstDensityThermo.h \
-             SpeciesThermoFactory.h \
-             ThermoFactory.h NasaPoly1.h NasaThermo.h ShomateThermo.h \
-             ShomatePoly.h SimpleThermo.h SpeciesThermoMgr.h \
+             SpeciesThermoFactory.h ThermoFactory.h \
+             NasaPoly1.h NasaPoly2.h NasaThermo.h \
+             ShomateThermo.h ShomatePoly.h ConstCpPoly.h\
+             SimpleThermo.h SpeciesThermoMgr.h \
+             GeneralSpeciesThermo.h Mu0Poly.h \
              speciesThermoTypes.h SpeciesThermo.h SurfPhase.h \
              EdgePhase.h polyfit.h  Func1.h \
              FuncEval.h StoichManager.h @phase_header_files@
diff --git a/Cantera/src/Mu0Poly.cpp b/Cantera/src/Mu0Poly.cpp
new file mode 100644
index 000000000..ff590f952
--- /dev/null
+++ b/Cantera/src/Mu0Poly.cpp
@@ -0,0 +1,409 @@
+/**
+ *  @file Mu0Poly.h
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+
+#include "Mu0Poly.h"
+#include "ctexceptions.h"
+#include "speciesThermoTypes.h"
+#include "SpeciesThermo.h"
+#include "xml.h"
+#include "ctml.h"
+
+using namespace ctml;
+
+namespace Cantera {
+
+    /**
+     *   The Mu0Poly class implements a linear interpolation
+     *   of the standard state chemical potential of one
+     *   species at a single reference pressure.
+     *   The chemical potential is input as a series of (T, mu0)
+     *   values. The first temperature is assumed to be equal
+     *   to 298.15 K; however, this may be relaxed in the future. 
+     *   This information, and an assumption of a constant
+     *   heat capacity within each interval is enough to 
+     *   calculate all thermodynamic functions.
+     *   
+     *  The basic equation for going from point 1 to point 2
+     *  are as follows for T, T1 <= T <= T2
+     *
+     *      mu1 = H1 - T1 * S1
+     * 
+     *      mu2 - mu1 = Cp1(T2 - T1) - Cp1(ln(T2/T1)) - S1(T2 - T1)
+     *
+     *      S2 = S1 + Cp1(ln(T2/T1))
+     * 
+     *      H2 = H1 + Cp1(T2 - T1)
+     *
+     *   In the future, a better assumption about the heat
+     *   capacity may be employed, so that it can be continuous.
+     *
+     *   Notes about temperature interpolation for T < T1 and T > Tn
+     *     These are achieved by assuming a constant heat capacity
+     *     equal to the value in the closest temperature interval.
+     *     No error is thrown.
+     */
+    Mu0Poly::Mu0Poly() : m_numIntervals(0),
+			 m_H298(0.0),
+			 m_lowT(0.0),
+			 m_highT(0.0),
+			 m_Pref(0.0),
+			 m_index(0) {
+    }
+
+    /**
+     * Mu0Poly():
+     *
+     * In the constructor, we calculate and store the
+     * piecewise linear approximation to the thermodynamic
+     * functions.
+     *
+     *  coeffs[0] = number of points (integer)
+     *         1  = H298(J/kmol)
+     *         2  = T1  (Kelvin)
+     *         3  = mu1 (J/kmol) 
+     *         4  = T2  (Kelvin)
+     *         5  = mu2 (J/kmol) 
+     *         6  = T3  (Kelvin)
+     *         7  = mu3 (J/kmol)
+     *         ........
+     */
+    Mu0Poly::Mu0Poly(int n, doublereal tlow, doublereal thigh, 
+		     doublereal pref,
+		     const doublereal* coeffs) :
+	m_numIntervals(0),
+	m_H298(0.0),
+	m_lowT      (tlow),
+	m_highT     (thigh),
+	m_Pref      (pref),
+	m_index     (n) {
+
+	int i, iindex;
+	double T1, T2;
+	int nPoints = (int) coeffs[0];
+	if (nPoints < 2) {
+	  throw CanteraError("Mu0Poly", 
+			     "nPoints must be >= 2");
+	}
+	m_numIntervals = nPoints - 1;
+	m_H298       = coeffs[1] / GasConstant;
+        int iT298 = 0;
+	/*
+	 * Resize according to the number of points
+	 */
+	m_t0_int.resize(nPoints);
+	m_h0_R_int.resize(nPoints);
+	m_s0_R_int.resize(nPoints);
+	m_cp0_R_int.resize(nPoints);
+	m_mu0_R_int.resize(nPoints);
+	/*
+	 * Calculate the T298 interval and make sure that
+	 * the temperatures are strictly monotonic.
+	 * Also distribute the data into the internal arrays.
+	 */
+        bool ifound = false;
+        for (i = 0, iindex = 2; i < nPoints; i++) {
+	  T1 = coeffs[iindex];
+	  m_t0_int[i] = T1;
+	  m_mu0_R_int[i] =  coeffs[iindex+1] / GasConstant;
+	  if (T1 == 298.15) {
+	    iT298 = i;
+	    ifound = true;
+	  }
+	  if (i < nPoints - 1) {
+	    T2 = coeffs[iindex+2];
+	    if (T2 <= T1) {
+	      throw CanteraError("Mu0Poly", 
+				 "Temperatures are not monotonic increasing");
+	    }
+	  }
+	  iindex += 2;
+	}
+	if (!ifound) {
+	  throw CanteraError("Mu0Poly",
+			     "One temperature has to be 298.15");
+	}
+
+	/*
+	 * Starting from the interval with T298, we go up
+	 */
+	doublereal mu2, s1, s2, h1, h2, cpi, deltaMu, deltaT;
+	T1 = m_t0_int[iT298];
+	doublereal mu1 = m_mu0_R_int[iT298];
+	m_h0_R_int[iT298] = m_H298;
+	m_s0_R_int[iT298] = - (mu1 - m_h0_R_int[iT298]) / T1;
+	for (i = iT298; i < m_numIntervals; i++) {
+	  T1      = m_t0_int[i];
+	  s1      = m_s0_R_int[i];
+	  h1      = m_h0_R_int[i];
+	  mu1     = m_mu0_R_int[i];
+	  T2      = m_t0_int[i+1];
+	  mu2     = m_mu0_R_int[i+1];
+	  deltaMu = mu2 - mu1;
+	  deltaT  = T2 - T1;
+	  cpi = (deltaMu - T1 * s1 + T2 * s1) / (deltaT - T2 * log(T2/T1));
+	  h2 = h1 + cpi * deltaT;
+	  s2 = s1 + cpi * log(T2/T1);
+	  m_cp0_R_int[i]   = cpi;
+	  m_h0_R_int[i+1]  = h2;
+	  m_s0_R_int[i+1]  = s2;
+	  m_cp0_R_int[i+1] = cpi;
+	}
+
+	/*
+	 * Starting from the interval with T298, we go down
+	 */
+	if (iT298 > 0) {
+	  T2 = m_t0_int[iT298];
+	  mu2 = m_mu0_R_int[iT298];
+	  m_h0_R_int[iT298] = m_H298;
+	  m_s0_R_int[iT298] = - (mu2 - m_h0_R_int[iT298]) / T2;
+	  for (i = iT298 - 1; i >= 0; i--) {
+	    T1      = m_t0_int[i];
+	    mu1     = m_mu0_R_int[i];
+	    T2      = m_t0_int[i+1];
+	    mu2     = m_mu0_R_int[i+1];
+	    s2      = m_s0_R_int[i+1];
+	    h2      = m_h0_R_int[i+1];
+	    deltaMu = mu2 - mu1;
+	    deltaT  = T2 - T1;
+	    cpi = (deltaMu - T1 * s2 + T2 * s2) / (deltaT - T1 * log(T2/T1)); 
+	    h1 = h2 - cpi * deltaT;
+	    s1 = s2 - cpi * log(T2/T1);
+	    m_cp0_R_int[i]   = cpi;
+	    m_h0_R_int[i]    = h1;
+	    m_s0_R_int[i]    = s1;
+	    if (i == (m_numIntervals-1)) {
+	      m_cp0_R_int[i+1] = cpi;
+	    }
+	  }
+	}
+#ifdef DEBUG_HKM_NOT
+	printf("    Temp     mu0(J/kmol)   cp0(J/kmol/K)   "
+	       " h0(J/kmol)   s0(J/kmol/K) \n");
+        for (i = 0; i < nPoints; i++) {
+	  printf("%12.3g %12.5g %12.5g %12.5g %12.5g\n",
+		 m_t0_int[i],  m_mu0_R_int[i] * GasConstant,
+		 m_cp0_R_int[i]* GasConstant, 
+		 m_h0_R_int[i]* GasConstant,
+		 m_s0_R_int[i]* GasConstant);
+	  fflush(stdout);
+	}
+#endif
+    }
+
+
+    Mu0Poly::Mu0Poly(const Mu0Poly &b)
+	: m_numIntervals (b.m_numIntervals),
+	  m_H298         (b.m_H298),
+	  m_t0_int       (b.m_t0_int),
+          m_mu0_R_int    (b.m_mu0_R_int),
+	  m_h0_R_int     (b.m_h0_R_int),
+	  m_s0_R_int     (b.m_s0_R_int),
+	  m_cp0_R_int    (b.m_cp0_R_int),
+	  m_lowT         (b.m_lowT),
+	  m_highT        (b.m_highT),
+	  m_Pref         (b.m_Pref),
+	  m_index        (b.m_index) {
+    }
+
+    Mu0Poly& Mu0Poly::operator=(const Mu0Poly& b) {
+	if (&b != this) {
+	  m_numIntervals = b.m_numIntervals;
+	  m_H298         = b.m_H298;
+	  m_t0_int       = b.m_t0_int;
+          m_mu0_R_int    = b.m_mu0_R_int;
+	  m_h0_R_int     = b.m_h0_R_int;
+	  m_s0_R_int     = b.m_s0_R_int;
+	  m_cp0_R_int    = b.m_cp0_R_int;
+	  m_lowT         = b.m_lowT;
+	  m_highT        = b.m_highT;
+	  m_Pref         = b.m_Pref;
+	  m_index        = b.m_index;
+	}
+	return *this;
+    }
+    
+    /**
+     * Destructor:
+     */
+    Mu0Poly::~Mu0Poly(){
+    }
+    
+    SpeciesThermoInterpType *
+    Mu0Poly::duplMyselfAsSpeciesThermoInterpType() const {
+	Mu0Poly* mp = new Mu0Poly(*this);
+	return (SpeciesThermoInterpType *) mp;
+    }
+	
+    doublereal Mu0Poly::minTemp() const { return m_lowT;}
+    doublereal  Mu0Poly::maxTemp() const { return m_highT;}
+    doublereal Mu0Poly::refPressure() const { return m_Pref; }
+    
+    /**
+     *  updateProperties is the main workhorse program. 
+     *  Given a temperature (*tt), it calculates the thermodynamic
+     *  functions H/RT, S_R, and cp_R, and returns the answer.
+     *
+     *  Note, it returns an answer by inserting the values into the
+     *  index position, m_index in vectors of  H/RT, S_R, and cp_R.
+     *
+     *
+     *   Input
+     *  -------
+     *       *tt = Temperature (Kelvin)
+     *      
+     */
+    void  Mu0Poly::
+    updateProperties(const doublereal* tt,  doublereal* cp_R,
+		     doublereal* h_RT, doublereal* s_R) const {
+	int j = m_numIntervals;
+	double T = *tt;
+	for (int i = 0; i < m_numIntervals; i++) {
+	  double T2 =  m_t0_int[i+1];
+	  if (T <=T2) {
+	    j = i;
+	    break;
+	  }
+	}
+	double T1 = m_t0_int[j];
+	double cp_Rj = m_cp0_R_int[j];
+
+	doublereal rt = 1.0/T;
+	cp_R[m_index] = cp_Rj;
+	h_RT[m_index] = rt*(m_h0_R_int[j] + (T - T1) * cp_Rj);
+	s_R[m_index]  = m_s0_R_int[j] + cp_Rj * (log(T/T1));
+    }
+
+    void  Mu0Poly::
+    updatePropertiesTemp(const doublereal T,
+			 doublereal* cp_R,
+			 doublereal* h_RT,
+			 doublereal* s_R) const {
+	updateProperties(&T, cp_R, h_RT, s_R);
+    }
+
+    /**
+     * report all of the parameters that make up this 
+     * interpolation.
+     */
+    void Mu0Poly::reportParameters(int &n, int &type,
+				   doublereal &tlow, doublereal &thigh,
+				   doublereal &pref,
+				   doublereal* const coeffs) const {
+	n = m_index;
+	type = MU0_INTERP;
+	tlow = m_lowT;
+	thigh = m_highT;
+	pref = m_Pref;
+	coeffs[0] = m_numIntervals+1;
+	coeffs[1] = m_H298 * GasConstant;
+	int j = 2;
+	for (int i = 0; i < m_numIntervals+1; i++) {
+	  coeffs[j] = m_t0_int[i];
+	  coeffs[j+1] = m_mu0_R_int[i] * GasConstant;
+	  j += 2;
+	}
+    }
+
+
+    /** 
+     * Install a Mu0 polynomial thermodynamic reference state property
+     * parameterization for species k into a SpeciesThermo instance,
+     * getting the information from an XML database.
+     */
+    void installMu0ThermoFromXML(string speciesName,
+				 SpeciesThermo& sp, int k, 
+				 const XML_Node* Mu0Node_ptr) {
+
+        doublereal tmin, tmax;
+	bool dimensionlessMu0Values = false;
+        const XML_Node& Mu0Node = *Mu0Node_ptr;
+  
+        tmin = fpValue(Mu0Node["Tmin"]);
+        tmax = fpValue(Mu0Node["Tmax"]);
+	doublereal pref = fpValue(Mu0Node["Pref"]);
+
+	doublereal h298 = 0.0;
+	if (Mu0Node.hasChild("H298")) {
+	  h298 = getFloat(Mu0Node, "H298", "actEnergy");
+	}
+
+	int numPoints = 1;
+	if (Mu0Node.hasChild("numPoints")) {
+	  numPoints = getInteger(Mu0Node, "numPoints");
+	}
+
+	vector_fp cValues(numPoints);
+	const XML_Node *valNode_ptr = 
+	    getByTitle(const_cast<XML_Node&>(Mu0Node), "Mu0Values");
+	if (!valNode_ptr) {
+	  throw CanteraError("installMu0ThermoFromXML",
+			     "missing required while processing "
+			     + speciesName);
+	}
+	getFloatArray(*valNode_ptr, cValues, true, "actEnergy");
+        /*
+	 * Check to see whether the Mu0's were input in a dimensionless
+	 * form. If they were, then the assumed temperature needs to be
+	 * adjusted from the assumed T = 273.15
+	 */
+	string uuu = (*valNode_ptr)["units"];
+	if (uuu == "Dimensionless") {
+	   dimensionlessMu0Values = true;
+	}
+	int ns = cValues.size();
+	if (ns != numPoints) {
+	  throw CanteraError("installMu0ThermoFromXML",
+			     "numPoints inconsistent while processing "
+			     + speciesName);
+	}
+	
+	vector_fp cTemperatures(numPoints);
+	const XML_Node *tempNode_ptr = 
+	    getByTitle(const_cast<XML_Node&>(Mu0Node), "Mu0Temperatures");
+	if (!tempNode_ptr) {
+	  throw CanteraError("installMu0ThermoFromXML",
+			     "missing required while processing + "
+			     + speciesName);
+	}
+	getFloatArray(*tempNode_ptr, cTemperatures, false);
+        ns = cTemperatures.size();
+	if (ns != numPoints) {
+	  throw CanteraError("installMu0ThermoFromXML",
+			     "numPoints inconsistent while processing "
+			     + speciesName);
+	}
+	
+	/*
+	 * Fix up dimensionless Mu0 values if input
+	 */
+        if (dimensionlessMu0Values) {
+	  for (int i = 0; i < numPoints; i++) {
+	    cValues[i] *= cTemperatures[i] / 273.15;
+	  }
+	}
+
+
+	vector_fp c(2 + 2 * numPoints);
+     
+	c[0] = numPoints;
+	c[1] = h298;
+	for (int i = 0; i < numPoints; i++) {
+	  c[2+i*2]   = cTemperatures[i];
+	  c[2+i*2+1] = cValues[i];
+	}
+      
+        sp.install(speciesName, k, MU0_INTERP, c.begin(), tmin, tmax, pref);
+    }
+}
+
+
+
+
+
diff --git a/Cantera/src/Mu0Poly.h b/Cantera/src/Mu0Poly.h
new file mode 100644
index 000000000..ddb3f09d1
--- /dev/null
+++ b/Cantera/src/Mu0Poly.h
@@ -0,0 +1,136 @@
+/**
+ *  @file Mu0Poly.h
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+
+
+#ifndef CT_MU0POLY_H
+#define CT_MU0POLY_H
+
+#include "SpeciesThermoInterpType.h"
+
+namespace Cantera {
+    class SpeciesThermo;
+    class XML_Node;
+    /**
+     *   The Mu0Poly class implements a linear interpolation
+     *   of the standard state chemical potential of one
+     *   species at a single reference pressure.
+     *   The chemical potential is input as a series of (T, mu0)
+     *   values. The first temperature is assumed to be equal
+     *   to 298.15 K; however, this may be relaxed in the future. 
+     *   This information, and an assumption of a constant
+     *   heat capacity within each interval is enough to 
+     *   calculate all thermodynamic functions.
+     *   
+     *   The basic equation for going from point 1 to point 2
+     *  are as follows for T, T1 <= T <= T2
+     *
+     *      mu1 = H1 - T1 * S1
+     * 
+     *      mu2 - mu1 = Cp1(T2 - T1) - Cp1(ln(T2/T1)) - S1(T2 - T1)
+     *
+     *      S2 = S1 + Cp1(ln(T2/T1))
+     * 
+     *      H2 = H1 + Cp1(T2 - T1)
+     *
+     *   In the future, a better assumption about the heat
+     *   capacity may be employed, so that it can be continuous.
+     *
+     *   Notes about temperature interpolation for T < T1 and T > Tn
+     *     These are achieved by assuming a constant heat capacity
+     *     equal to the value in the closest temperature interval.
+     *     No error is thrown.
+     */
+    class Mu0Poly: public SpeciesThermoInterpType {
+
+    public:
+
+        Mu0Poly();
+
+        Mu0Poly(int n, doublereal tlow, doublereal thigh, 
+		doublereal pref, const doublereal* coeffs);
+	Mu0Poly(const Mu0Poly &);
+	Mu0Poly& operator=(const Mu0Poly&);
+        virtual ~Mu0Poly();
+	SpeciesThermoInterpType *
+	duplMyselfAsSpeciesThermoInterpType() const;
+	
+        doublereal minTemp() const;
+        doublereal maxTemp() const;
+        doublereal refPressure() const;
+        virtual int reportType() const { return MU0_INTERP; }
+
+	/**
+	 * Update all of the properties, using the polynomial
+	 * tPoly[]
+	 *
+	 * tPoly[0] = temp (Kelvin)
+	 */
+        void updateProperties(const doublereal* tPoly, 
+			      doublereal* cp_R, doublereal* h_RT, 
+			      doublereal* s_R) const ;
+
+	void updatePropertiesTemp(const doublereal temp, 
+				  doublereal* cp_R,
+				  doublereal* h_RT, 
+				  doublereal* s_R) const ;
+
+	/**
+	 * report all of the parameters that make up this 
+	 * interpolation.
+	 */
+	void reportParameters(int &n, int &type,
+			      doublereal &tlow, doublereal &thigh,
+			      doublereal &pref,
+			      doublereal* const coeffs) const;
+    protected:
+	/**
+	 * Number of intervals in the interpolating linear
+	 * approximation. Number of points is one more than the
+	 * number of intervals.
+	 */
+	int m_numIntervals;
+	/**
+	 * Value of the enthalpy at T = 298.15.
+	 *  This value is tied to the Heat of formation of 
+	 *  the species at 298.15.
+	 */
+	doublereal m_H298;
+	/**
+	 * Points at which the standard state chemical potential
+	 * are given. 
+	 */
+	vector_fp m_t0_int;
+
+	/*
+	 * Mu0's are primary input data. They aren't strictly
+	 * needed, but are kept here for convenience. 
+	 */
+	vector_fp m_mu0_R_int;
+	vector_fp m_h0_R_int;
+	vector_fp m_s0_R_int;
+	vector_fp m_cp0_R_int;
+        doublereal m_lowT, m_highT, m_Pref;
+        int m_index;
+    
+    private:
+
+    };
+
+    void installMu0ThermoFromXML(string speciesName,
+				 SpeciesThermo& sp, int k, 
+				 const XML_Node* Mu0Node_ptr);
+}
+
+#endif
+
+
+
+
+
+
diff --git a/Cantera/src/NasaPoly1.h b/Cantera/src/NasaPoly1.h
index 1b3921244..0fc4983f3 100755
--- a/Cantera/src/NasaPoly1.h
+++ b/Cantera/src/NasaPoly1.h
@@ -13,6 +13,8 @@
 #ifndef CT_NASAPOLY1_H
 #define CT_NASAPOLY1_H
 
+#include "SpeciesThermoInterpType.h"
+
 namespace Cantera {
 
     /**
@@ -41,7 +43,7 @@ namespace Cantera {
      * This class is designed specifically for use by class NasaThermo.
      * @ingroup spthermo
      */
-    class NasaPoly1 {
+    class NasaPoly1 : public SpeciesThermoInterpType {
 
     public:
 
@@ -59,20 +61,59 @@ namespace Cantera {
             copy(coeffs, coeffs + 7, m_coeff.begin());
         }
 
-        ~NasaPoly1(){}
+	NasaPoly1(const NasaPoly1& b) :
+	    m_lowT      (b.m_lowT),
+            m_highT     (b.m_highT),
+            m_Pref      (b.m_Pref),
+            m_index     (b.m_index),
+            m_coeff     (array_fp(7)) {
+            copy(b.m_coeff.begin(),
+		 b.m_coeff.begin() + 7,
+		 m_coeff.begin());
+	}
+
+	NasaPoly1& operator=(const NasaPoly1& b) {
+	    if (&b != this) {
+	      m_lowT   = b.m_lowT;
+	      m_highT  = b.m_highT;
+	      m_Pref   = b.m_Pref;
+	      m_index  = b.m_index;
+	      copy(b.m_coeff.begin(),
+		   b.m_coeff.begin() + 7,
+		   m_coeff.begin());
+	    }
+	    return *this;
+	}
+
+        virtual ~NasaPoly1(){}
+
+	virtual SpeciesThermoInterpType *
+	duplMyselfAsSpeciesThermoInterpType() const {
+	    NasaPoly1* np = new NasaPoly1(*this);
+	    return (SpeciesThermoInterpType *) np;
+	}
       
         doublereal minTemp() const     { return m_lowT;}
         doublereal maxTemp() const     { return m_highT;}
         doublereal refPressure() const { return m_Pref; }
+        virtual int reportType() const { return NASA1; }
       
         /**
          * Update the properties for this species. This method is called 
          * with a pointer to an array containing the functions of
-	 *  temperature needed by this 
+	 * temperature needed by this 
          * parameterization, and three pointers to arrays where the
-	 *  computed property values
+	 * computed property values
          * should be written. This method updates only one value in 
 	 * each array.  
+	 *
+	 * Temperature Polynomial:
+	 *  tt[0] = t;
+	 *  tt[1] = t*t;
+	 *  tt[2] = m_t[1]*t;
+	 *  tt[3] = m_t[2]*t;
+	 *  tt[4] = 1.0/t;
+	 *  tt[5] = log(t);
          */
         void updateProperties(const doublereal* tt, 
             doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
@@ -97,10 +138,32 @@ namespace Cantera {
             s_R[m_index] = s;
         }
 
-	void reportParameters(int &n, doublereal &tlow, doublereal &thigh,
+	/**
+	 *   updatePropertiesTemp():
+	 *     This formulation creates its own temperature 
+	 *     polynomial. Then, it calls updateProperties();
+	 *
+	 *     (note: this is slow, but it is general)
+	 */
+	void updatePropertiesTemp(const doublereal temp, 
+				  doublereal* cp_R, doublereal* h_RT, 
+				  doublereal* s_R) const {
+	    double tPoly[6];
+	    tPoly[0]  = temp;
+	    tPoly[1]  = temp * temp;
+	    tPoly[2]  = tPoly[1] * temp;
+	    tPoly[3]  = tPoly[2] * temp;
+	    tPoly[4]  = 1.0 / temp;
+	    tPoly[5]  = log(temp);
+	    updateProperties(tPoly, cp_R, h_RT, s_R);
+	}
+
+	void reportParameters(int &n, int &type,
+			      doublereal &tlow, doublereal &thigh,
 			      doublereal &pref,
 			      doublereal* const coeffs) const {
 	    n = m_index;
+	    type = NASA1;
 	    tlow = m_lowT;
 	    thigh = m_highT;
 	    pref = m_Pref;
@@ -122,11 +185,5 @@ namespace Cantera {
     };
 
 }
-
 #endif
 
-
-
-
-
-
diff --git a/Cantera/src/NasaPoly2.h b/Cantera/src/NasaPoly2.h
new file mode 100644
index 000000000..3acec7970
--- /dev/null
+++ b/Cantera/src/NasaPoly2.h
@@ -0,0 +1,217 @@
+/**
+ *  @file NasaPoly1.h
+ */
+
+/* $Author$
+ * $Revision$
+ * $Date$
+ */
+
+// Copyright 2001  California Institute of Technology
+
+
+#ifndef CT_NASAPOLY2_H
+#define CT_NASAPOLY2_H
+
+#include "SpeciesThermoInterpType.h"
+
+namespace Cantera {
+
+    /**
+     *
+     * 
+     * The NASA polynomial parameterization for one temperature range.
+     * This parameterization expresses the heat capacity as a
+     * fourth-order polynomial. Note that this is the form used in the
+     * 1971 NASA equilibrium program and by the Chemkin software
+     * package, but differs from the form used in the more recent NASA
+     * equilibrium program.
+     *
+     * Seven coefficients \f$(a_0,\dots,a_6)\f$ are used to represent
+     * \f$ c_p^0(T)\f$, \f$ h^0(T)\f$, and \f$ s^0(T) \f$ as 
+     * polynomials in \f$ T \f$ :  
+     * \f[
+     * \frac{c_p(T)}{R} = a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4
+     * \f]
+     * \f[
+     * \frac{h^0(T)}{RT} = a_0 + \frac{a_1}{2} T + \frac{a_2}{3} T^2 
+     * + \frac{a_3}{4} T^3 + \frac{a_4}{5} T^4  + \frac{a_5}{T}.
+     * \f]
+     * \f[
+     * \frac{s^0(T)}{R} = a_0\ln T + a_1 T + \frac{a_2}{2} T^2 
+     + \frac{a_3}{3} T^3 + \frac{a_4}{4} T^4  + a_6.
+     * \f]
+     * 
+     * This class is designed specifically for use by class 
+     * GeneralSpeciesThermo.
+     * @ingroup spthermo
+     */
+    class NasaPoly2 : public SpeciesThermoInterpType {
+
+    public:
+
+        NasaPoly2() 
+            : m_lowT(0.0),
+	      m_midT(0.0), 
+	      m_highT (0.0),
+              m_Pref(0.0),
+	      mnp_low(0),
+	      mnp_high(0), 
+	      m_index(0),
+	      m_coeff(array_fp(15)) {
+	}
+
+        NasaPoly2(int n, doublereal tlow, doublereal thigh, doublereal pref,
+            const doublereal* coeffs) :
+            m_lowT(tlow),
+            m_highT(thigh),
+            m_Pref(pref),
+	    mnp_low(0),
+	    mnp_high(0),
+            m_index(n),
+            m_coeff(array_fp(15)) {
+
+            copy(coeffs, coeffs + 15, m_coeff.begin());
+	    m_midT = coeffs[0];
+	    mnp_low  = new NasaPoly1(m_index, m_lowT, m_midT,
+				     m_Pref, m_coeff.begin()+1);
+	    mnp_high = new NasaPoly1(m_index, m_midT, m_highT,
+				     m_Pref, m_coeff.begin()+8);
+        }
+
+	NasaPoly2(const NasaPoly2& b) :
+	    m_lowT(b.m_lowT),
+	    m_midT(b.m_midT),
+            m_highT(b.m_highT),
+            m_Pref(b.m_Pref),
+	    mnp_low(0),
+	    mnp_high(0),
+            m_index(b.m_index),
+            m_coeff(array_fp(15)) {
+
+            copy(b.m_coeff.begin(),
+		 b.m_coeff.begin() + 15,
+		 m_coeff.begin());
+	    mnp_low  = new NasaPoly1(m_index, m_lowT, m_midT,
+				     m_Pref, m_coeff.begin()+1);
+	    mnp_high = new NasaPoly1(m_index, m_midT, m_highT,
+				     m_Pref, m_coeff.begin()+8);
+	}
+
+	NasaPoly2& operator=(const NasaPoly2& b) {
+	    if (&b != this) {
+	      m_lowT   = b.m_lowT;
+	      m_midT   = b.m_midT;
+	      m_highT  = b.m_highT;
+	      m_Pref   = b.m_Pref;
+	      m_index  = b.m_index;
+	      copy(b.m_coeff.begin(),
+		   b.m_coeff.begin() + 15,
+		   m_coeff.begin());
+	      if (mnp_low) delete mnp_low;
+	      if (mnp_high) delete mnp_high;
+	      mnp_low  = new NasaPoly1(m_index, m_lowT, m_midT, 
+				       m_Pref, m_coeff.begin()+1);
+	      mnp_high = new NasaPoly1(m_index, m_midT, m_highT, 
+				       m_Pref, m_coeff.begin()+8);
+	    }
+	    return *this;
+	}
+
+        virtual ~NasaPoly2(){
+	    delete mnp_low;
+	    delete mnp_high;
+	}
+
+	virtual SpeciesThermoInterpType *
+	duplMyselfAsSpeciesThermoInterpType() const {
+	    NasaPoly2* np = new NasaPoly2(*this);
+	    return (SpeciesThermoInterpType *) np;
+	}
+      
+        doublereal minTemp() const     { return m_lowT;}
+        doublereal maxTemp() const     { return m_highT;}
+        doublereal refPressure() const { return m_Pref; }
+        virtual int reportType() const { return NASA2; }
+      
+        /**
+         * Update the properties for this species. This method is called 
+         * with a pointer to an array containing the functions of
+	 * temperature needed by this 
+         * parameterization, and three pointers to arrays where the
+	 * computed property values
+         * should be written. This method updates only one value in 
+	 * each array.  
+	 *
+	 * Temperature Polynomial:
+	 *  tt[0] = t;
+	 *  tt[1] = t*t;
+	 *  tt[2] = m_t[1]*t;
+	 *  tt[3] = m_t[2]*t;
+	 *  tt[4] = 1.0/t;
+	 *  tt[5] = log(t);
+         */
+        void updateProperties(const doublereal* tt, 
+            doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
+          
+       	    double T = tt[0];
+	    if (T <= m_midT) {
+	      mnp_low->updateProperties(tt, cp_R, h_RT, s_R);
+	    } else {
+	      mnp_high->updateProperties(tt, cp_R, h_RT, s_R);
+	    }
+        }
+
+	/**
+	 *   updatePropertiesTemp():
+	 *     This formulation creates its own temperature 
+	 *     polynomial. Then, it calls updateProperties();
+	 *
+	 *     (note: this is slow, but it is general)
+	 */
+	void updatePropertiesTemp(const doublereal temp, 
+				  doublereal* cp_R,
+				  doublereal* h_RT, 
+				  doublereal* s_R) const {
+	    if (temp <= m_midT) {
+	      mnp_low->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
+	    } else {
+	      mnp_high->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
+	    }
+	}
+
+	void reportParameters(int &n, int &type,
+			      doublereal &tlow, doublereal &thigh,
+			      doublereal &pref,
+			      doublereal* const coeffs) const {
+	    n = m_index;
+	    type = NASA2;
+	    tlow = m_lowT;
+	    thigh = m_highT;
+	    pref = m_Pref;
+	    for (int i = 0; i < 15; i++) {
+	      coeffs[i] = m_coeff[i];
+	    }
+	}
+
+    protected:
+        
+        doublereal m_lowT;     // lowest valid temperature
+	doublereal m_midT;
+        doublereal m_highT;    // highest valid temperature
+        doublereal m_Pref;     // standard-state pressure
+	NasaPoly1 *mnp_low;
+	NasaPoly1 *mnp_high;
+        int m_index;           // species index
+        array_fp m_coeff;      // array of polynomial coefficients
+    
+    private:
+
+    };
+
+}
+#endif
+
+
+
+
diff --git a/Cantera/src/NasaThermo.h b/Cantera/src/NasaThermo.h
index 8416eb1c6..31aba5a18 100755
--- a/Cantera/src/NasaThermo.h
+++ b/Cantera/src/NasaThermo.h
@@ -65,7 +65,8 @@ namespace Cantera {
          * - c[1] - c[7]   coefficients for low T range
          * - c[8] - c[14]  coefficients for high T range
          */
-        virtual void install(string name, int index, int type, const doublereal* c, 
+        virtual void install(string name, int index, int type, 
+			     const doublereal* c, 
 			     doublereal minTemp, doublereal maxTemp,
 			     doublereal refPressure) { 
 
@@ -184,7 +185,9 @@ namespace Cantera {
                 return m_thigh[k];
         }
 
-        virtual doublereal refPressure() const {return m_p0;}
+        virtual doublereal refPressure(int k = -1) const {
+            return m_p0;
+        }
 
 	/**
          * This utility function reports the type of parameterization
@@ -211,21 +214,27 @@ namespace Cantera {
 	      const vector<NasaPoly1> &mhg = m_high[grp-1];
 	      const NasaPoly1 *lowPoly  = &(mlg[pos]);
 	      const NasaPoly1 *highPoly = &(mhg[pos]);
-
+	      int itype = NASA;
 	      doublereal tmid = lowPoly->maxTemp();
 	      c[0] = tmid;
 	      int n;
 	      double ttemp;
-	      lowPoly->reportParameters(n, minTemp, ttemp, refPressure,
+	      lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
 					c + 1);
 	      if (n != index) {
 		throw CanteraError("  ", "confused");
 	      }
-	      highPoly->reportParameters(n, ttemp, maxTemp, refPressure,
+	      if (itype != NASA1) {
+		throw CanteraError("  ", "confused");
+	      }
+	      highPoly->reportParameters(n, itype, ttemp, maxTemp, refPressure,
 					c + 8);
 	      if (n != index) {
 		throw CanteraError("  ", "confused");
 	      }
+	      if (itype != NASA1) {
+		throw CanteraError("  ", "confused");
+	      }
 	    } else {
 	      throw CanteraError(" ", "confused");
 	    }
diff --git a/Cantera/src/ShomatePoly.h b/Cantera/src/ShomatePoly.h
index 05bbb9a11..3f6c754d3 100755
--- a/Cantera/src/ShomatePoly.h
+++ b/Cantera/src/ShomatePoly.h
@@ -12,6 +12,8 @@
 #ifndef CT_SHOMATEPOLY1_H
 #define CT_SHOMATEPOLY1_H
 
+#include "SpeciesThermoInterpType.h"
+
 namespace Cantera {
 
     /**
@@ -32,7 +34,7 @@ namespace Cantera {
      * \f]
      */
 
-    class ShomatePoly {
+    class ShomatePoly : public SpeciesThermoInterpType {
 
     public:
 
@@ -50,23 +52,59 @@ namespace Cantera {
             copy(coeffs, coeffs + 7, m_coeff.begin());
         }
 
-        ~ShomatePoly(){}
+	ShomatePoly(const ShomatePoly& b) :
+	    m_lowT      (b.m_lowT),
+            m_highT     (b.m_highT),
+            m_Pref      (b.m_Pref),
+            m_coeff     (array_fp(7)),
+	    m_index     (b.m_index) {
+            copy(b.m_coeff.begin(),
+		 b.m_coeff.begin() + 7,
+		 m_coeff.begin());
+	}
+
+	ShomatePoly& operator=(const ShomatePoly& b) {
+	    if (&b != this) {
+	      m_lowT   = b.m_lowT;
+	      m_highT  = b.m_highT;
+	      m_Pref   = b.m_Pref;
+	      m_index  = b.m_index;
+	      copy(b.m_coeff.begin(),
+		   b.m_coeff.begin() + 7,
+		   m_coeff.begin());
+	    }
+	    return *this;
+	}
+        virtual ~ShomatePoly(){}
+
+	virtual SpeciesThermoInterpType *
+	duplMyselfAsSpeciesThermoInterpType() const {
+	    ShomatePoly* sp = new ShomatePoly(*this);
+	    return (SpeciesThermoInterpType *) sp;
+	}
       
         doublereal minTemp() const { return m_lowT;}
         doublereal maxTemp() const { return m_highT;}
         doublereal refPressure() const { return m_Pref; }
+        virtual int reportType() const { return SHOMATE; }
       
         /**
-         * t is T/1000.
-         *
-         * tt[0]      t
-         * tt[1]      t*t
-         * tt[2]      t*t*t
-         * tt[3]      t^4
-         * tt[4]      ln t
+	 *  This formulation calculates the thermo functions
+	 *  given the native formulation of the temperature
+	 *  polynomial 
+	 *
+         *  tt is T/1000. 
+	 *  m_t[0] = tt;
+	 *  m_t[1] = tt*tt;
+	 *  m_t[2] = m_t[1]*tt;
+	 *  m_t[3] = 1.0/m_t[1];
+	 *  m_t[4] = log(tt);
+	 *  m_t[5] = 1.0/GasConstant;
+	 *  m_t[6] = 1.0/(GasConstant * T);
          */
         void updateProperties(const doublereal* tt, 
-            doublereal* cp_R, doublereal* h_RT, doublereal* s_R) const {
+			      doublereal* cp_R, doublereal* h_RT, 
+			      doublereal* s_R) const {
 
             doublereal A      = m_coeff[0]; 
             doublereal Bt     = m_coeff[1]*tt[0];    
@@ -80,18 +118,49 @@ namespace Cantera {
             cp = A + Bt + Ct2 + Dt3 + Etm2;
             h = tt[0]*(A + 0.5*Bt + OneThird*Ct2 + 0.25*Dt3 - Etm2) + F;
             s = A*tt[4] + Bt + 0.5*Ct2 + OneThird*Dt3 - 0.5*Etm2 + G; 
-            h *= 1.e6;
 
+            /*
+             *  Shomate polynomials parameterizes assuming units of
+             *  J/(gmol*K) for cp_r and s_R and kJ/(gmol) for h.
+             *  However, Cantera assumes default MKS units of 
+             *  J/(kmol*K). This requires us to multiply cp and s
+             *  by 1.e3 and h by 1.e6, before we then nondimensionlize
+             *  the results by dividing by (GasConstant * T),
+             *  where GasConstant has units of J/(kmol * K).
+             */
             cp_R[m_index] = 1.e3 * cp * tt[5];
-            h_RT[m_index] = h * tt[6];
-            s_R[m_index] = 1.e3 * s * tt[5];
+            h_RT[m_index] = 1.e6 * h  * tt[6];
+            s_R[m_index]  = 1.e3 * s  * tt[5];
         }
 
+	/**
+	 *   updatePropertiesTemp():
+	 *     This formulation creates its own temperature 
+	 *     polynomial. Then, it calls updateProperties();
+	 *     -> general, but slow.
+	 */
+	void updatePropertiesTemp(const doublereal temp, 
+				  doublereal* cp_R, doublereal* h_RT, 
+				  doublereal* s_R) const {
+	    double tPoly[7];
+	    doublereal tt = 1.e-3*temp;
+	    tPoly[0] = tt;
+	    tPoly[1] = tt * tt;
+	    tPoly[2] = tPoly[1] * tt;
+	    tPoly[3] = 1.0/tPoly[1];
+	    tPoly[4] = log(tt);
+	    tPoly[5] = 1.0/GasConstant;
+	    tPoly[6] = 1.0/(GasConstant * temp);
+	    updateProperties(tPoly, cp_R, h_RT, s_R);
+	}
 
-	void reportParameters(int &n, doublereal &tlow, doublereal &thigh,
+
+	void reportParameters(int &n, int &type,
+			      doublereal &tlow, doublereal &thigh,
 			      doublereal &pref,
 			      doublereal* const coeffs) const {
 	    n = m_index;
+	    type = SHOMATE;
 	    tlow = m_lowT;
 	    thigh = m_highT;
 	    pref = m_Pref;
@@ -109,12 +178,159 @@ namespace Cantera {
 
     };
 
+
+
+    class ShomatePoly2 : public SpeciesThermoInterpType {
+    public:
+
+	ShomatePoly2() 
+            : m_lowT(0.0),
+	      m_midT(0.0), 
+	      m_highT (0.0),
+	      m_Pref(0.0),
+	      msp_low(0),
+	      msp_high(0),
+              m_index(0) {
+	    m_coeff.resize(15);
+	}
+
+	ShomatePoly2(int n, doublereal tlow, doublereal thigh, doublereal pref,
+		    const doublereal* coeffs) :
+            m_lowT      (tlow),
+	    m_midT(0.0),
+            m_highT     (thigh),
+	    m_Pref      (pref),
+	    msp_low(0),
+	    msp_high(0),
+            m_index     (n)  {
+            m_coeff.resize(15);
+            copy(coeffs, coeffs + 15, m_coeff.begin());
+	    m_midT = coeffs[0];
+	    msp_low  = new ShomatePoly(n, tlow, m_midT, pref, coeffs+1);
+	    msp_high = new ShomatePoly(n, m_midT, thigh, pref, coeffs+8);
+        }
+
+	ShomatePoly2(const ShomatePoly2& b) :
+	    m_lowT      (b.m_lowT),
+	    m_midT      (b.m_midT),
+            m_highT     (b.m_highT),
+            m_Pref      (b.m_Pref),
+	    msp_low(0),
+	    msp_high(0),
+            m_coeff     (array_fp(15)),
+	    m_index     (b.m_index) {
+            copy(b.m_coeff.begin(),
+		 b.m_coeff.begin() + 15,
+		 m_coeff.begin());
+	    msp_low  = new ShomatePoly(m_index, m_lowT, m_midT, 
+				       m_Pref, m_coeff.begin()+1);
+	    msp_high = new ShomatePoly(m_index, m_midT, m_highT, 
+				       m_Pref, m_coeff.begin()+8);
+	}
+
+	ShomatePoly2& operator=(const ShomatePoly2& b) {
+	    if (&b != this) {
+	      m_lowT   = b.m_lowT;
+	      m_midT   = b.m_midT;
+	      m_highT  = b.m_highT;
+	      m_Pref   = b.m_Pref;
+	      m_index  = b.m_index;
+	      copy(b.m_coeff.begin(),
+		   b.m_coeff.begin() + 15,
+		   m_coeff.begin());
+	      if (msp_low) delete msp_low;
+	      if (msp_high) delete msp_high;
+	      msp_low  = new ShomatePoly(m_index, m_lowT, m_midT, 
+					 m_Pref, m_coeff.begin()+1);
+	      msp_high = new ShomatePoly(m_index, m_midT, m_highT, 
+					 m_Pref, m_coeff.begin()+8);
+	    }
+	    return *this;
+	}
+
+	virtual ~ShomatePoly2(){
+	    delete msp_low;
+	    delete msp_high;
+	}
+	
+	virtual SpeciesThermoInterpType *
+	duplMyselfAsSpeciesThermoInterpType() const {
+	    ShomatePoly2* sp = new ShomatePoly2(*this);
+	    return (SpeciesThermoInterpType *) sp;
+	}
+
+	doublereal minTemp() const { return m_lowT;}
+        doublereal maxTemp() const { return m_highT;}
+        doublereal refPressure() const { return m_Pref; }
+        virtual int reportType() const { return SHOMATE2; }
+
+        /**
+	 *  This formulation calculates the thermo functions
+	 *  given the native formulation of the temperature
+	 *  polynomial 
+	 *
+         *  tt is T/1000. 
+	 *  m_t[0] = tt;
+	 *  m_t[1] = tt*tt;
+	 *  m_t[2] = m_t[1]*tt;
+	 *  m_t[3] = 1.0/m_t[1];
+	 *  m_t[4] = log(tt);
+	 *  m_t[5] = 1.0/GasConstant;
+	 *  m_t[6] = 1.0/(GasConstant * T);
+         */
+	void updateProperties(const doublereal* tt, 
+			      doublereal* cp_R, doublereal* h_RT, 
+			      doublereal* s_R) const {
+	    double T = 1000 * tt[0];
+	    if (T <= m_midT) {
+	      msp_low->updateProperties(tt, cp_R, h_RT, s_R);
+	    } else {
+	      msp_high->updateProperties(tt, cp_R, h_RT, s_R);
+	    }
+	    
+	}
+
+	/**
+	 *   updatePropertiesTemp():
+	 *     This formulation creates its own temperature 
+	 *     polynomial. Then, it calls updateProperties();
+	 *     -> general, but slow.
+	 */
+	void updatePropertiesTemp(const doublereal temp, 
+				  doublereal* cp_R,
+				  doublereal* h_RT, 
+				  doublereal* s_R) const {
+	    if (temp <= m_midT) {
+	      msp_low->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
+	    } else {
+	      msp_high->updatePropertiesTemp(temp, cp_R, h_RT, s_R);
+	    }
+	}
+
+	void reportParameters(int &n, int &type,
+			      doublereal &tlow, doublereal &thigh,
+			      doublereal &pref,
+			      doublereal* const coeffs) const {
+	    n = m_index;
+	    type = SHOMATE2;
+	    tlow = m_lowT;
+	    thigh = m_highT;
+	    pref = m_Pref;
+	    for (int i = 0; i < 15; i++) {
+	      coeffs[i] = m_coeff[i];
+	    }
+	}
+
+    protected:
+        
+        doublereal m_lowT, m_midT;
+	doublereal m_highT;
+	doublereal m_Pref;
+	ShomatePoly *msp_low;
+	ShomatePoly *msp_high;
+        array_fp m_coeff;
+        int m_index;
+    };
 }
 
 #endif
-
-
-
-
-
-
diff --git a/Cantera/src/ShomateThermo.h b/Cantera/src/ShomateThermo.h
index b3d0054d6..6249e98f8 100755
--- a/Cantera/src/ShomateThermo.h
+++ b/Cantera/src/ShomateThermo.h
@@ -10,12 +10,20 @@
  *
  *  \f[ S = A*ln(t) + B*t + C*t^2/2 + D*t^3/3 - E/(2*t^2) + G \f]
  *
- *    - Cp = heat capacity (J/mol*K)
- *    - H = standard enthalpy (kJ/mol)
- *    - \f$ \Delta_f H_298.15 \f$ = enthalpy of formation at 298.15 K (kJ/mol)
- *    - S = standard entropy (J/mol*K)
+ *    - Cp = heat capacity (J/gmol*K)
+ *    - H = standard enthalpy (kJ/gmol)
+ *    - \f$ \Delta_f H_298.15 \f$ = enthalpy of formation at 298.15 K (kJ/gmol)
+ *    - S = standard entropy (J/gmol*K)
  *    - t = temperature (K) / 1000.
  *
+ *  Note, the polynomial data (i.e., A, ... , G) is entered in dimensional form.
+ *  This is in contrast to the NASA database polynomials which are entered in 
+ *  nondimensional form (i.e., NASA parameterizes C_p/R, while Shomate
+ *  parameterizes C_p assuming units of J/gmol*K - and kJ/gmol*K for H).
+ *  Note, also that the H - H_298.15 equation has units of kJ/gmol, because of
+ *  the implicit integration of (t = T 1000), which provides a 
+ *  multiplier of 1000 to the Enthalpy equation.
+ *
  */
 
 #ifndef CT_SHOMATETHERMO_H
@@ -54,8 +62,10 @@ namespace Cantera {
          * in the same units as used in the NIST Chemistry WebBook.
          *
          */
-        virtual void install(string name, int index, int type, const doublereal* c,
-            doublereal minTemp, doublereal maxTemp, doublereal refPressure) {
+        virtual void install(string name, int index, int type, 
+			     const doublereal* c,
+			     doublereal minTemp, doublereal maxTemp, 
+			     doublereal refPressure) {
             int imid = int(c[0]);       // midpoint temp converted to integer
             int igrp = m_index[imid];   // has this value been seen before?
             if (igrp == 0) {            // if not, prepare new group
@@ -139,7 +149,7 @@ namespace Cantera {
                     _end    = m_low[i].end();
                 }
                 for (; _begin != _end; ++_begin) {
-                    _begin->updateProperties(m_t.begin(), cp_R, h_RT, s_R);
+                  _begin->updateProperties(m_t.begin(), cp_R, h_RT, s_R);
                 }
             }
         }
@@ -158,7 +168,9 @@ namespace Cantera {
                 return m_thigh[k];
         }
 
-        virtual doublereal refPressure() const {return m_p0;}
+        virtual doublereal refPressure(int k=-1) const {
+            return m_p0;
+        }
 
 	virtual int reportType(int index) const { return SHOMATE; }
 
@@ -174,9 +186,10 @@ namespace Cantera {
 				  doublereal &maxTemp, 
 				  doublereal &refPressure) {
 	    type = reportType(index);
-	    if (type == NASA) {
+	    if (type == SHOMATE) {
 	      int grp = m_group_map[index];
 	      int pos = m_posInGroup_map[index];
+	      int itype = SHOMATE;
 	      const vector<ShomatePoly> &mlg = m_low[grp-1];
 	      const vector<ShomatePoly> &mhg = m_high[grp-1];
 	      const ShomatePoly *lowPoly  = &(mlg[pos]);
@@ -185,16 +198,22 @@ namespace Cantera {
 	      c[0] = tmid;
 	      int n;
 	      double ttemp;
-	      lowPoly->reportParameters(n, minTemp, ttemp, refPressure,
+	      lowPoly->reportParameters(n, itype, minTemp, ttemp, refPressure,
 					c + 1);
 	      if (n != index) {
 		throw CanteraError("  ", "confused");
 	      }
-	      highPoly->reportParameters(n, ttemp, maxTemp, refPressure,
-					c + 8);
+	      if (itype != SHOMATE && itype != SHOMATE1) {
+		throw CanteraError("  ", "confused");
+	      }
+	      highPoly->reportParameters(n, itype,  ttemp, maxTemp,
+					 refPressure, c + 8);
 	      if (n != index) {
 		throw CanteraError("  ", "confused");
 	      }
+	      if (itype != SHOMATE && itype != SHOMATE1) {
+		throw CanteraError("  ", "confused");
+	      }
 	    } else {
 	      throw CanteraError(" ", "confused");
 	    }
diff --git a/Cantera/src/SimpleThermo.h b/Cantera/src/SimpleThermo.h
index 2d534d922..aba2bbc4e 100644
--- a/Cantera/src/SimpleThermo.h
+++ b/Cantera/src/SimpleThermo.h
@@ -27,7 +27,8 @@ namespace Cantera {
 
         virtual ~SimpleThermo() {}
 
-        virtual void install(string name, int index, int type, const doublereal* c,
+        virtual void install(string name, int index, int type, 
+			     const doublereal* c,
             doublereal minTemp, doublereal maxTemp, doublereal refPressure) {
             m_logt0.push_back(log(c[0]));
             m_t0.push_back(c[0]);
@@ -80,7 +81,7 @@ namespace Cantera {
                 return m_thigh[k];
         }
 
-        virtual doublereal refPressure() const {return m_p0;}
+        virtual doublereal refPressure(int k=-1) const {return m_p0;}
 
 	virtual int reportType(int index) const { return SIMPLE; }
 
diff --git a/Cantera/src/SpeciesThermo.h b/Cantera/src/SpeciesThermo.h
index a93d5bf7a..56c972bb2 100755
--- a/Cantera/src/SpeciesThermo.h
+++ b/Cantera/src/SpeciesThermo.h
@@ -64,14 +64,18 @@ namespace Cantera {
          * parameterization. 
          * @see speciesThermoTypes.h 
          */
-        virtual void install(string name, int index, int type, const doublereal* c, 
-            doublereal minTemp, doublereal maxTemp, doublereal refPressure)=0;
+        virtual void install(string name, int index, int type, 
+			     const doublereal* c, 
+			     doublereal minTemp,
+			     doublereal maxTemp,
+			     doublereal refPressure)=0;
 
         /**
          * Compute the standard-state properties for all species.
          * Given temperature T in K, this method updates the values of
          * the non-dimensional heat capacity at constant pressure,
-         * enthalpy, and entropy.
+         * enthalpy, and entropy, at the reference pressure Pref
+         * of each of the standard states.
          */
         virtual void update(doublereal T, 
             doublereal* cp_R, 
@@ -79,7 +83,7 @@ namespace Cantera {
             doublereal* s_R) const=0;
 
         /**
-         * Like update, but only updates the species k.
+         * Like update(), but only updates the single species k.
          */
         virtual void update_one(int k, doublereal T, 
             doublereal* cp_R, 
@@ -107,10 +111,17 @@ namespace Cantera {
         virtual doublereal maxTemp(int k=-1) const =0;
 
         /**
-         * The standard-state pressure. All parameterizations must be
-         * for the same standard-state pressure.
+         * The reference-state pressure for species k.
+         *
+         * returns the reference state pressure in Pascals for
+         * species k. If k is left out of the argument list,
+         * it returns the reference state pressure for the first
+         * species.
+         * Note that some SpeciesThermo implementations, such
+         * as those for ideal gases, require that all species
+         * in the same phase have the same reference state pressures.
          */
-        virtual doublereal refPressure() const =0;
+        virtual doublereal refPressure(int k=-1) const =0;
 
 	/**
          * This utility function reports the type of parameterization
diff --git a/Cantera/src/SpeciesThermoFactory.cpp b/Cantera/src/SpeciesThermoFactory.cpp
index 20ff7732a..eb5c59982 100755
--- a/Cantera/src/SpeciesThermoFactory.cpp
+++ b/Cantera/src/SpeciesThermoFactory.cpp
@@ -16,6 +16,8 @@
 #include "ShomateThermo.h"
 //#include "PolyThermoMgr.h"
 #include "SimpleThermo.h"
+#include "GeneralSpeciesThermo.h"
+#include "Mu0Poly.h"
 
 #include "SpeciesThermoMgr.h"
 #include "speciesThermoTypes.h"
@@ -31,7 +33,8 @@ namespace Cantera {
 
 
     static void getSpeciesThermoTypes(XML_Node* node, 
-        int& has_nasa, int& has_shomate, int& has_simple) {
+				      int& has_nasa, int& has_shomate, int& has_simple,
+				      int &has_other) {
         const XML_Node& sparray = *node;
         vector<XML_Node*> sp;
         sparray.getChildren("species",sp);
@@ -44,10 +47,11 @@ namespace Cantera {
             if (th.hasChild("Shomate")) has_shomate = 1;
             if (th.hasChild("const_cp")) has_simple = 1;
             if (th.hasChild("poly")) {
-                if (th.child("poly")["order"] == "1") has_simple = 1;
-                else throw CanteraError("newSpeciesThermo",
-                    "poly with order > 1 not yet supported");
+	      if (th.child("poly")["order"] == "1") has_simple = 1;
+	      else throw CanteraError("newSpeciesThermo",
+				      "poly with order > 1 not yet supported");
             }
+	    if (th.hasChild("Mu0")) has_other = 1;
 	  } else {
 	    throw UnknownSpeciesThermoModel("getSpeciesThermoTypes:",
 					    spNode->attrib("name"), "missing");
@@ -61,8 +65,15 @@ namespace Cantera {
      * specified in a CTML phase specification.
      */
     SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(XML_Node* node) {
-        int inasa = 0, ishomate = 0, isimple = 0;
-        getSpeciesThermoTypes(node, inasa, ishomate, isimple);
+        int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
+	try {
+	  getSpeciesThermoTypes(node, inasa, ishomate, isimple, iother);
+	} catch (UnknownSpeciesThermoModel) {
+	  iother = 1;
+	}
+	if (iother) {
+	  return new GeneralSpeciesThermo();
+	}
         return newSpeciesThermo(NASA*inasa
             + SHOMATE*ishomate + SIMPLE*isimple);
     }
@@ -70,10 +81,17 @@ namespace Cantera {
     SpeciesThermo* SpeciesThermoFactory::
     newSpeciesThermo(vector<XML_Node*> nodes) {
         int n = static_cast<int>(nodes.size());
-        int inasa = 0, ishomate = 0, isimple = 0;
+        int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
         for (int j = 0; j < n; j++) {
-            getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple);
+	  try {
+            getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother);
+	  } catch (UnknownSpeciesThermoModel) {
+	    iother = 1;
+	  }
         }
+	if (iother) {
+	  return new GeneralSpeciesThermo();
+	}
         return newSpeciesThermo(NASA*inasa
             + SHOMATE*ishomate + SIMPLE*isimple);
     }
@@ -82,14 +100,18 @@ namespace Cantera {
     SpeciesThermo* SpeciesThermoFactory::
     newSpeciesThermoOpt(vector<XML_Node*> nodes) {
         int n = static_cast<int>(nodes.size());
-        int inasa = 0, ishomate = 0, isimple = 0;
+        int inasa = 0, ishomate = 0, isimple = 0, iother = 0;
         for (int j = 0; j < n; j++) {
 	  try {
-            getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple);
+            getSpeciesThermoTypes(nodes[j], inasa, ishomate, isimple, iother);
 	  } catch (UnknownSpeciesThermoModel) {
-
+	    iother = 1;
+	    popError();
 	  }
         }
+	if (iother) {
+	  return new GeneralSpeciesThermo();
+	}
         return newSpeciesThermo(NASA*inasa
             + SHOMATE*ishomate + SIMPLE*isimple);
     }
@@ -98,7 +120,6 @@ namespace Cantera {
     SpeciesThermo* SpeciesThermoFactory::newSpeciesThermo(int type) {
         
         switch (type) {
-
         case NASA:
             return new NasaThermo;
         case SHOMATE:
@@ -308,8 +329,10 @@ namespace Cantera {
             getFloatArray(f0.child("floatArray"), c0, false);
             if (dualRange)
                 getFloatArray(f1ptr->child("floatArray"), c1, false);
-            else
-                c1.resize(7,0.0);
+            else {
+	      c1.resize(7,0.0);
+	      copy(c0.begin(), c0.begin()+7, c1.begin());
+	    }
         }
         else if (fabs(tmax1 - tmin0) < 0.01) {
             tmin = tmin1;
@@ -385,6 +408,9 @@ namespace Cantera {
             else if (f->name() == "NASA") {
                 installNasaThermoFromXML(s["name"], spthermo, k, f, 0);
             }
+	    else if (f->name() == "Mu0") {
+	      installMu0ThermoFromXML(s["name"], spthermo, k, f);
+	    }
             else {
                 throw UnknownSpeciesThermoModel("installSpecies", 
 						s["name"], f->name());
diff --git a/Cantera/src/SpeciesThermoInterpType.h b/Cantera/src/SpeciesThermoInterpType.h
new file mode 100644
index 000000000..375caca65
--- /dev/null
+++ b/Cantera/src/SpeciesThermoInterpType.h
@@ -0,0 +1,49 @@
+/**
+ *  @file SpeciesThermoInterpType.h
+ *
+ * $Author$
+ * $Revision$
+ * $Date$
+ */
+
+// Copyright 2001  California Institute of Technology
+
+#include "speciesThermoTypes.h"
+#ifndef CT_SPECIESTHERMOINTERPTYPE_H
+#define CT_SPECIESTHERMOINTERPTYPE_H
+
+namespace Cantera {
+
+    
+    class SpeciesThermoInterpType {
+
+    public:
+    
+        virtual SpeciesThermoInterpType * 
+          duplMyselfAsSpeciesThermoInterpType() const = 0;
+       
+        virtual doublereal minTemp() const = 0;
+        virtual doublereal maxTemp() const = 0;
+	virtual doublereal refPressure() const = 0;
+        virtual int reportType() const = 0;
+
+        virtual void updateProperties(const doublereal* tempPoly, 
+				      doublereal* cp_R,
+				      doublereal* h_RT,
+				      doublereal* s_R) const = 0;
+
+        virtual void updatePropertiesTemp(const doublereal temp, 
+				      doublereal* cp_R,
+				      doublereal* h_RT,
+				      doublereal* s_R) const = 0;
+
+	virtual void reportParameters(int &n, int &type,
+				      doublereal &tlow, doublereal &thigh,
+				      doublereal &pref,
+				      doublereal* const coeffs) const = 0;
+    };
+
+}
+
+#endif
+
diff --git a/Cantera/src/SpeciesThermoMgr.h b/Cantera/src/SpeciesThermoMgr.h
index 7861443a7..8e3d4f0f6 100755
--- a/Cantera/src/SpeciesThermoMgr.h
+++ b/Cantera/src/SpeciesThermoMgr.h
@@ -104,14 +104,19 @@ namespace Cantera {
 	SpeciesThermoDuo() {}
         virtual ~SpeciesThermoDuo(){}
                 
-	virtual void install(string name, int sp, int type, const doublereal* c,
-            doublereal minTemp, doublereal maxTemp, doublereal refPressure) {
+	virtual void install(string name, int sp, int type,
+			     const doublereal* c,
+			     doublereal minTemp,
+			     doublereal maxTemp, 
+			     doublereal refPressure) {
             m_p0 = refPressure;
             if (type == m_thermo1.ID) {
-                m_thermo1.install(name, sp, 0, c, minTemp, maxTemp, refPressure);
+                m_thermo1.install(name, sp, 0, c, minTemp, maxTemp,
+				  refPressure);
 		speciesToType[sp] = m_thermo1.ID;
             } else if (type == m_thermo2.ID) {
-                m_thermo2.install(name, sp, 0, c, minTemp, maxTemp, refPressure);
+                m_thermo2.install(name, sp, 0, c, minTemp, maxTemp,
+				  refPressure);
 		speciesToType[sp] = m_thermo2.ID;
             } else {
                 throw UnknownSpeciesThermo("SpeciesThermoDuo:install",type);
@@ -136,7 +141,7 @@ namespace Cantera {
             return (tm1 < tm2 ? tm1 : tm2);
         }                        
 
-        virtual doublereal refPressure() const {
+        virtual doublereal refPressure(int k = -1) const {
             return m_p0;
         }
 
@@ -158,7 +163,7 @@ namespace Cantera {
 	    if (ctype == m_thermo1.ID) {
 	      m_thermo1.reportParams(index, type, c, minTemp, maxTemp, 
 				     refPressure);
-	    } else if (ctype == m_thermo1.ID) {
+	    } else if (ctype == m_thermo2.ID) {
 	      m_thermo2.reportParams(index, type, c, minTemp, maxTemp, 
 				     refPressure);
 	    } else {
@@ -225,7 +230,7 @@ namespace Cantera {
                 return m_thermo[k].maxTemp();
         }                        
 
-        virtual doublereal refPressure() const {
+        virtual doublereal refPressure(int k = -1) const {
             return m_pref;
         }
 
diff --git a/Cantera/src/speciesThermoTypes.h b/Cantera/src/speciesThermoTypes.h
index 47166cbbc..95a04d498 100755
--- a/Cantera/src/speciesThermoTypes.h
+++ b/Cantera/src/speciesThermoTypes.h
@@ -20,15 +20,23 @@
 
 // NASA Polynomials
 #define NASA 4
+#define NASA2 4
 
 // Shomate Polynomials used in NIST database
 #define SHOMATE 8
+#define SHOMATE2 8
 
 // Tiger Polynomials
 #define TIGER 16
 
 #define SIMPLE 32
 
+#define MU0_INTERP 64
+
+#define SHOMATE1 128
+
+#define NASA1  256
+
 #include "ct_defs.h"
 
 namespace Cantera {