Merge pull request #9 from adelq/janaf-burcat

JANAF and Burcat fixes

Author: adelq

.travis.yml

@@ -3,6 +3,7 @@ python:
   - "2.7"
   - "3.4"
   - "3.5"
+  - "3.6"
 
 install:
   - sudo apt-get update

tests/__init__.py

@@ -11,5 +11,5 @@ def test_titanium():
 
 def test_rutile():
     db = Janafdb()
-    rutile = db.getphasedata(formula='O2Ti', name='Rutile', phase='cr')
+    rutile = db.getphasedata(formula='O2Ti', name='Rutile', phase='cr', cache=False)
     assert_allclose(rutile.cp([500, 550, 1800]), [67.203, 68.567, 78.283])

tests/test_janaf.py

@@ -1 +1 @@
-__version__ = '0.6.0'
+__version__ = '0.7.0'

thermochem/__init__.py

@@ -15,7 +15,7 @@
     from xml.etree.ElementTree import parse
 except ImportError:
     from elementtree import parse
-from numpy import empty, array, dot, log
+import numpy as np
 
 # Universal gas constant R
 R = 8.314472
@@ -62,11 +62,11 @@ def cpo(self, T):
         # I know perfectly that the most efficient way of evaluatin
         # polynomials is recursively but I want the implementation to
         # be as explicit as possible
-        Ta = array([1, T, T ** 2, T ** 3, T ** 4], 'd')
+        Ta = np.array([1, T, T ** 2, T ** 3, T ** 4], 'd')
         if T > 200 and T <= 1000:
-            return dot(self.Tmin_[:5], Ta) * R
+            return np.dot(self.Tmin_[:5], Ta) * R
         elif T > 1000 and T < 6000:
-            return dot(self._Tmax[:5], Ta) * R
+            return np.dot(self._Tmax[:5], Ta) * R
         else:
             raise ValueError("Temperature out of range")
 
@@ -87,11 +87,11 @@ def ho(self, T):
         """
         Computes the sensible enthalpy in J/mol
         """
-        Ta = array([1, T / 2, T ** 2 / 3, T ** 3 / 4, T ** 4 / 5, 1 / T], 'd')
+        Ta = np.array([1, T / 2, T ** 2 / 3, T ** 3 / 4, T ** 4 / 5, 1 / T], 'd')
         if T > 200 and T <= 1000:
-            return dot(self.Tmin_[:6], Ta) * R * T
+            return np.dot(self.Tmin_[:6], Ta) * R * T
         elif T > 1000 and T < 6000:
-            return dot(self._Tmax[:6], Ta) * R * T
+            return np.dot(self._Tmax[:6], Ta) * R * T
         else:
             raise ValueError("Temperature out of range")
 
@@ -105,11 +105,11 @@ def so(self, T):
         """
         Computes enthropy in J/mol K
         """
-        Ta = array([log(T), T, T ** 2 / 2, T ** 3 / 3, T ** 4 / 4, 0, 1], 'd')
+        Ta = np.array([np.log(T), T, T ** 2 / 2, T ** 3 / 3, T ** 4 / 4, 0, 1], 'd')
         if T > 200 and T <= 1000:
-            return dot(self.Tmin_, Ta) * R
+            return np.dot(self.Tmin_, Ta) * R
         elif T > 1000 and T < 6000:
-            return dot(self._Tmax, Ta) * R
+            return np.dot(self._Tmax, Ta) * R
         else:
             raise ValueError("Temperature out of range")
 
@@ -324,7 +324,6 @@ def __unicode__(self):
 
 
 class Elementdb(object):
-
     """
     Class that reads the Alexander Burcat's thermochemical database
     for combustion.
@@ -416,41 +415,32 @@ def getelementdata(self, formula):
         """
         Returns an element instance given the name of the element.
         """
-        Tmin_ = empty((7), 'd')
-        _Tmax = empty((7), 'd')
+        Tmin_ = np.zeros(7)
+        _Tmax = np.zeros(7)
         comp = []
         for specie in self.db:
-            try:
-                for element in specie:
-                    try:
-                        if element.tag == "phase":
-                            if formula == element.find("formula").text:
-                                phase = element
-                                coefficients = phase.find("coefficients")
-                                low = coefficients.find("range_Tmin_to_1000")
-                                for (i, c) in zip(range(7), low):
-                                    Tmin_[i] = float(c.text)
-
-                                high = coefficients.find("range_1000_to_Tmax")
-                                for (i, c) in zip(range(7), high):
-                                    _Tmax[i] = float(c.text)
-
-                                elements = phase.find("elements")
-                                elements = elements.getchildren()
-                                for elem in elements:
-                                    it = elem.items()
-                                    # First is name of element, second is number
-                                    # of atoms
-                                    comp.append((it[0][1], int(it[1][1])))
-
-                                mm = float(phase.find("molecular_weight").text) / 1000
-                                hfr = float(coefficients.find("hf298_div_r").text)
-
-                                return Element(formula, Tmin_, _Tmax, mm, hfr, comp)
-                    except:
-                        pass
-            except:
-                pass
+            for element in specie:
+                if element.tag == "phase":
+                    if formula == element.find("formula").text:
+                        phase = element
+                        coefficients = phase.find("coefficients")
+                        low = coefficients.find("range_Tmin_to_1000")
+                        for (i, c) in zip(range(7), low):
+                            Tmin_[i] = float(c.text)
+
+                        high = coefficients.find("range_1000_to_Tmax")
+                        for (i, c) in zip(range(7), high):
+                            _Tmax[i] = float(c.text)
+
+                        elements = phase.find("elements").getchildren()
+                        for elem in elements:
+                            elem_data = elem.attrib
+                            comp.append((elem_data['name'], int(elem_data['num_of_atoms'])))
+
+                        mm = float(phase.find("molecular_weight").text) / 1000
+                        hfr = float(coefficients.find("hf298_div_r").text)
+
+                        return Element(formula, Tmin_, _Tmax, mm, hfr, comp)
 
     def getmixturedata(self, components):
         """

thermochem/burcat.py

@@ -11,6 +11,7 @@
 from __future__ import print_function
 
 import os
+import sys
 import numpy as np
 import pandas as pd
 from scipy.interpolate import interp1d
@@ -71,8 +72,9 @@ def __init__(self, rawdata_text):
             StringIO(self.rawdata_text),
             skiprows=2,
             header=None,
-            delimiter=r'[\t\s]*',
-            engine='python')
+            delimiter=r'[\t\s]+',
+            engine='python'
+        )
         data.columns = ['T', 'Cp', 'S', '[G-H(Tr)]/T', 'H-H(Tr)', 'Delta_fH',
                         'Delta_fG', 'log(Kf)']
         self.rawdata = data
@@ -86,7 +88,7 @@ def __init__(self, rawdata_text):
             if np.issubdtype(data.dtypes[c], np.floating):
                 continue
             # Change INFINITE to inf
-            data.loc[data[c] == 'INFINITE', c]
+            data.loc[data[c] == 'INFINITE', c] = np.inf
             # Anything else becomes a nan.
             # Convert to floats.
             data[c] = pd.to_numeric(data[c], errors='coerce')
@@ -162,11 +164,11 @@ def search(self, searchstr):
 
         return self.db[formulasearch | namesearch]
 
-    def getphasedata(self, formula=None, name=None, phase=None, nocache=False):
+    def getphasedata(self, formula=None, name=None, phase=None, cache=True):
         """
         Returns an element instance given the name of the element.
         formula, name and phase match the respective fields in the JANAF index.
-        nocache = True means that we will always get the data from the web.
+        cache = False means that we will always get the data from the web.
 
         >>> db = Janafdb()
         >>> db.getphasedata(formula='O2Ti', phase='cr')
@@ -184,9 +186,8 @@ def getphasedata(self, formula=None, name=None, phase=None, nocache=False):
         """
 
         # Check that the phase type requested is valid.
-        if phase not in Janafdb.VALIDPHASETYPES:
-            raise ValueError("Valid phase types are " + str(
-                Janafdb.VALIDPHASETYPES) + ".")
+        if phase not in self.VALIDPHASETYPES:
+            raise ValueError("Valid phase types are %s." % self.VALIDPHASETYPES)
 
         # We can search on either an exact formula, partial text match in the
         # name, and exact phase type.
@@ -212,8 +213,10 @@ def getphasedata(self, formula=None, name=None, phase=None, nocache=False):
 
         # At this point we have one record.  Check if we have that file cached.
         cachedfilename = os.path.join(
-            self.JANAF_cachedir, PhaseRecord['filename'].values[0] + '.txt')
-        if os.path.exists(cachedfilename) and not nocache:
+            self.JANAF_cachedir,
+            "%s.txt" % PhaseRecord['filename'].values[0]
+        )
+        if cache and os.path.exists(cachedfilename):
             # Yes it was cached, so let's read it into memory.
             with open(cachedfilename, 'r') as f:
                 textdata = f.read()
@@ -222,10 +225,12 @@ def getphasedata(self, formula=None, name=None, phase=None, nocache=False):
             response = urllib2.urlopen(Janafdb.JANAF_URL %
                                        PhaseRecord['filename'].values[0])
             textdata = response.read()
+            if sys.version_info[0] > 2:
+                textdata = textdata.decode()
 
             # And cache the data so we aren't making unnecessary trips to the
             # web.
-            if not nocache:
+            if cache:
                 with open(cachedfilename, 'w') as f:
                     f.write(textdata)