#!/usr/bin/env python
# encoding: utf-8

################################################################################
#
#   Copyright (c) 2009-2011 by the RMG Team (rmg_dev@mit.edu)
#
#   Permission is hereby granted, free of charge, to any person obtaining a
#   copy of this software and associated documentation files (the 'Software'),
#   to deal in the Software without restriction, including without limitation
#   the rights to use, copy, modify, merge, publish, distribute, sublicense,
#   and/or sell copies of the Software, and to permit persons to whom the
#   Software is furnished to do so, subject to the following conditions:
#
#   The above copyright notice and this permission notice shall be included in
#   all copies or substantial portions of the Software.
#
#   THE SOFTWARE IS PROVIDED 'AS IS', WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
#   IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
#   FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
#   AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
#   LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
#   FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
#   DEALINGS IN THE SOFTWARE.
#
################################################################################

"""
This module contains functions for converting Chemkin-format input files to
Cantera input files (CTI).
"""

import logging
import types
import os.path
import numpy as np
import re

QUANTITY_UNITS = {'MOL': 'mol',
                  'MOLE': 'mol',
                  'MOLES': 'mol',
                  'MOLEC': 'molec',
                  'MOLECULES': 'molec'}

ENERGY_UNITS = {'CAL/': 'cal/mol',
                'CAL/MOL': 'cal/mol',
                'CAL/MOLE': 'cal/mol',
                'EVOL': 'eV',
                'EVOLTS': 'eV',
                'JOUL': 'J/mol',
                'JOULES/MOL': 'J/mol',
                'JOULES/MOLE': 'J/mol',
                'KCAL': 'kcal/mol',
                'KCAL/MOL': 'kcal/mol',
                'KCAL/MOLE': 'kcal/mol',
                'KELV': 'K',
                'KELVIN': 'K',
                'KELVINS': 'K',
                'KJOU': 'kJ/mol',
                'KJOULES/MOL': 'kJ/mol',
                'KJOULES/MOLE': 'kJ/mol'}


def compatible_quantities(quantity_basis, units):
    if quantity_basis == 'mol':
        return 'molec' not in units
    elif quantity_basis == 'molec':
        return 'molec' in units or 'mol' not in units
    else:
        raise Exception('Unknown quantity basis: "{0}"'.format(quantity_basis))


class InputParseError(Exception):
    """
    An exception class for exceptional behavior involving Chemkin-format
    mechanism files. Pass a string describing the circumstances that caused
    the exceptional behavior.
    """
    pass


class Species(object):
    def __init__(self, label):
        self.label = label
        self.thermo = None
        self.transport = None
        self.note = None
        self.composition = None

    def __str__(self):
        return self.label

    def to_cti(self, indent=0):
        lines = []
        atoms = ' '.join('{0}:{1}'.format(*a)
                         for a in self.composition.iteritems())

        prefix = ' '*(indent+8)

        lines.append('species(name={0!r},'.format(self.label))
        lines.append(prefix + 'atoms={0!r},'.format(atoms))
        if self.thermo:
            lines.append(prefix +
                         'thermo={0},'.format(self.thermo.to_cti(15+indent)))
        if self.transport:
            lines.append(prefix +
                         'transport={0},'.format(self.transport.to_cti(14+indent)))
        if self.note:
            lines.append(prefix + 'note={0!r},'.format(self.note))

        lines[-1] = lines[-1][:-1] + ')'
        lines.append('')

        return '\n'.join(lines)


class ThermoModel(object):
    """
    A base class for thermodynamics models, containing several attributes
    common to all models:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `Tmin`          ``float``           The minimum temperature at which the model is valid, or ``None`` if unknown or undefined
    `Tmax`          ``float``           The maximum temperature at which the model is valid, or ``None`` if unknown or undefined
    `comment`       ``str``             Information about the model (e.g. its source)
    =============== =================== ========================================

    """

    def __init__(self, Tmin=None, Tmax=None, comment=''):
        if Tmin is not None:
            self.Tmin = Tmin
        else:
            self.Tmin = None
        if Tmax is not None:
            self.Tmax = Tmax
        else:
            self.Tmax = None
        self.comment = comment


class NASA(ThermoModel):
    """
    A single NASA polynomial for thermodynamic data. The `coeffs` attribute
    stores the seven or nine polynomial coefficients
    :math:`\\mathbf{a} = \\left[a_{-2}\\ a_{-1}\\ a_0\\ a_1\\ a_2\\ a_3\\ a_4\\ a_5\\ a_6 \\right]`
    from which the relevant thermodynamic parameters are evaluated via the
    expressions

    .. math:: \\frac{C_\\mathrm{p}(T)}{R} = a_{-2} T^{-2} + a_{-1} T^{-1} + a_0 + a_1 T + a_2 T^2 + a_3 T^3 + a_4 T^4

    .. math:: \\frac{H(T)}{RT} = - a_{-2} T^{-2} + a_{-1} T^{-1} \\ln T + a_0 + \\frac{1}{2} a_1 T + \\frac{1}{3} a_2 T^2 + \\frac{1}{4} a_3 T^3 + \\frac{1}{5} a_4 T^4 + \\frac{a_5}{T}

    .. math:: \\frac{S(T)}{R} = -\\frac{1}{2} a_{-2} T^{-2} - a_{-1} T^{-1} + a_0 \\ln T + a_1 T + \\frac{1}{2} a_2 T^2 + \\frac{1}{3} a_3 T^3 + \\frac{1}{4} a_4 T^4 + a_6

    For the 7 coefficient form, the first two coefficients are taken to be zero.
    """

    def __init__(self, coeffs, **kwargs):
        ThermoModel.__init__(self, **kwargs)
        if len(coeffs) not in (7,9):
            raise InputParseError('Invalid number of NASA polynomial coefficients; '
                                  'should be 7 or 9.')
        self.coeffs = coeffs

    def to_cti(self, indent=0):
        prefix = ' '*indent
        vals = ['{0: 15.8E}'.format(i) for i in self.coeffs]
        if len(self.coeffs) == 7:
            lines = ['NASA([{0:.2f}, {1:.2f}],'.format(self.Tmin[0], self.Tmax[0]),
                     prefix+'     [{0}, {1}, {2},'.format(*vals[0:3]),
                     prefix+'      {0}, {1}, {2},'.format(*vals[3:6]),
                     prefix+'      {0}]),'.format(vals[6])]
        else:
            lines = ['NASA9([{0:.2f}, {1:.2f}],'.format(self.Tmin[0], self.Tmax[0]),
                     prefix+'      [{0}, {1}, {2},'.format(*vals[0:3]),
                     prefix+'       {0}, {1}, {2},'.format(*vals[3:6]),
                     prefix+'       {0}, {1}, {2}]),'.format(*vals[6:9])]

        return '\n'.join(lines)


class MultiNASA(ThermoModel):
    """
    A set of thermodynamic parameters given by NASA polynomials. This class
    stores a list of :class:`NASA` objects in the `polynomials`
    attribute. When evaluating a thermodynamic quantity, a polynomial that
    contains the desired temperature within its valid range will be used.
    """

    def __init__(self, polynomials=None, **kwargs):
        ThermoModel.__init__(self, **kwargs)
        self.polynomials = polynomials or []

    def to_cti(self, indent=0):
        prefix = ' '*indent
        lines = []
        for i,p in enumerate(self.polynomials):
            if i == 0:
                lines.append('({0}'.format(p.to_cti(indent+1)))
            elif i != len(self.polynomials)-1:
                lines.append(prefix + ' {0}'.format(p.to_cti(indent+1)))
            else:
                lines.append(prefix + ' {0})'.format(p.to_cti(indent+1)[:-1]))

        return '\n'.join(lines)


class Reaction(object):
    """
    A chemical reaction. The attributes are:

    =================== =========================== ============================
    Attribute           Type                        Description
    =================== =========================== ============================
    `index`             :class:`int`                A unique nonnegative integer index
    `reactants`         :class:`list`               The reactant species (as :class:`Species` objects)
    `products`          :class:`list`               The product species (as :class:`Species` objects)
    `kinetics`          :class:`KineticsModel`      The kinetics model to use for the reaction
    `reversible`        ``bool``                    ``True`` if the reaction is reversible, ``False`` if not
    `duplicate`         ``bool``                    ``True`` if the reaction is known to be a duplicate, ``False`` if not
    `fwdOrders`         ``dict``                    Reaction order (value) for each specified species (key)
    =================== =========================== ============================

    """

    def __init__(self, index=-1, reactants=None, products=None, kinetics=None,
                 reversible=True, duplicate=False, fwdOrders=None):
        self.index = index
        self.reactants = reactants  # list of (stoichiometry, species) tuples
        self.products = products  # list of (stoichiometry, specis) tuples
        self.kinetics = kinetics
        self.reversible = reversible
        self.duplicate = duplicate
        self.fwdOrders = fwdOrders if fwdOrders is not None else {}
        self.thirdBody = None

    def _coeff_string(self, coeffs):
        L = []
        for stoichiometry,species in coeffs:
            if stoichiometry != 1:
                L.append('{0} {1}'.format(stoichiometry, species))
            else:
                L.append(str(species))
        expression = ' + '.join(L)
        if self.thirdBody:
            expression += ' (+ {0})'.format(self.thirdBody)

        return expression

    @property
    def reactantString(self):
        return self._coeff_string(self.reactants)

    @property
    def productString(self):
        return self._coeff_string(self.products)

    def __str__(self):
        """
        Return a string representation of the reaction, in the form 'A + B <=> C + D'.
        """
        arrow = ' <=> ' if self.reversible else ' -> '
        return arrow.join([self.reactantString, self.productString])

    def to_cti(self, indent=0):
        arrow = ' <=> ' if self.reversible else ' => '

        kinstr = self.kinetics.to_cti(self.reactantString, arrow,
                                      self.productString, indent)

        k_indent = ' ' * (kinstr.find('(') + 1)

        if self.duplicate:
            kinstr = kinstr[:-1] + ",\n{0}options='duplicate')".format(k_indent)

        if self.fwdOrders:
            order = ' '.join('{0}:{1}'.format(k,v)
                             for (k,v) in self.fwdOrders.items())
            kinstr = kinstr[:-1] + ",\n{0}order='{1}')".format(k_indent, order)

        return kinstr


class KineticsModel(object):
    """
    A base class for kinetics models, containing several attributes common to
    all models:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `Tmin`          :class:`Quantity`   The minimum absolute temperature in K at which the model is valid
    `Tmax`          :class:`Quantity`   The maximum absolute temperature in K at which the model is valid
    `Pmin`          :class:`Quantity`   The minimum absolute pressure in Pa at which the model is valid
    `Pmax`          :class:`Quantity`   The maximum absolute pressure in Pa at which the model is valid
    `comment`       :class:`str`        A string containing information about the model (e.g. its source)
    =============== =================== ========================================

    """

    def __init__(self, Tmin=None, Tmax=None, Pmin=None, Pmax=None, comment='',
                 parser=None):
        self.Tmin = Tmin
        self.Tmax = Tmax
        self.Pmin = Pmin
        self.Pmax = Pmax
        self.comment = comment
        self.parser = parser

    def isPressureDependent(self):
        """
        Return ``True`` if the kinetics are pressure-dependent or ``False`` if
        they are pressure-independent. This method must be overloaded in the
        derived class.
        """
        raise InputParseError('Unexpected call to KineticsModel.isPressureDependent();'
                              ' you should be using a class derived from KineticsModel.')

    def to_cti(self, reactantstr, arrow, productstr):
        raise InputParseError('to_cti is not implemented for objects of class {0}'.format(self.__class__.__name__))

    def efficiencyString(self):
        if hasattr(self, 'efficiencies'):
            return ' '.join('{0}:{1}'.format(mol, eff)
                            for mol,eff in self.efficiencies.iteritems())
        else:
            return ''


class KineticsData(KineticsModel):
    """
    A kinetics model based around a set of discrete (high-pressure limit)
    rate coefficients at various temperatures. The attributes are:

    =========== =================== ============================================
    Attribute   Type                Description
    =========== =================== ============================================
    `Tdata`     :class:`Quantity`   The temperatures at which the heat capacity data is provided
    `kdata`     :class:`Quantity`   The rate coefficients in SI units at each temperature in `Tdata`
    =========== =================== ============================================

    """

    def __init__(self, Tdata=None, kdata=None, **kwargs):
        KineticsModel.__init__(self, **kwargs)
        self.Tdata = Tdata
        self.kdata = kdata

    def isPressureDependent(self):
        """
        Returns ``False`` since KineticsData kinetics are not
        pressure-dependent.
        """
        return False


class Arrhenius(KineticsModel):
    """
    Represent a set of modified Arrhenius kinetics. The kinetic expression has
    the form

    .. math:: k(T) = A \\left( \\frac{T}{T_0} \\right)^n \\exp \\left( - \\frac{E_\\mathrm{a}}{RT} \\right)

    where :math:`A`, :math:`n`, :math:`E_\\mathrm{a}`, and :math:`T_0` are the
    parameters to be set, :math:`T` is absolute temperature, and :math:`R` is
    the gas law constant. The attributes are:

    =============== =================== ========================================
    Attribute       Type                Description
    =============== =================== ========================================
    `A`             :class:`Quantity`   The preexponential factor in s^-1, m^3/mol*s, etc.
    `T0`            :class:`Quantity`   The reference temperature in K
    `n`             :class:`Quantity`   The temperature exponent
    `Ea`            :class:`Quantity`   The activation energy in J/mol
    =============== =================== ========================================

    """

    def __init__(self, A=0.0, n=0.0, Ea=0.0, T0=1.0, **kwargs):
        KineticsModel.__init__(self, **kwargs)
        self.A = A
        self.T0 = T0
        self.n = n
        self.Ea = Ea

    def isPressureDependent(self):
        """
        Returns ``False`` since Arrhenius kinetics are not pressure-dependent.
        """
        return False

    def rateStr(self):
        if compatible_quantities(self.parser.quantity_units, self.A[1]):
            A = '{0:e}'.format(self.A[0])
        else:
            A = "({0:e}, '{1}')".format(*self.A)

        if self.Ea[1] == self.parser.energy_units:
            Ea = str(self.Ea[0])
        else:
            Ea = "({0}, '{1}')".format(*self.Ea)

        return '[{0}, {1}, {2}]'.format(A, self.n, Ea)

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstring = reactantstr + arrow + productstr
        return 'reaction({0!r}, {1})'.format(rxnstring, self.rateStr())


class PDepArrhenius(KineticsModel):
    """
    A kinetic model of a phenomenological rate coefficient k(T, P) using the
    expression

    .. math:: k(T,P) = A(P) T^{n(P)} \\exp \\left[ \\frac{-E_\\mathrm{a}(P)}{RT} \\right]

    where the modified Arrhenius parameters are stored at a variety of pressures
    and interpolated between on a logarithmic scale. The attributes are:

    =============== ================== ============================================
    Attribute       Type               Description
    =============== ================== ============================================
    `pressures`     :class:`list`      The list of pressures in Pa
    `arrhenius`     :class:`list`      The list of :class:`Arrhenius` objects at each pressure
    `highPlimit`    :class:`Arrhenius` The high (infinite) pressure limiting :class:`Arrhenius` expression
    =============== ================== ============================================

    Note that `highPlimit` is not used in evaluating k(T,P).
    """

    def __init__(self, pressures=None, arrhenius=None, highPlimit=None, **kwargs):
        KineticsModel.__init__(self, **kwargs)
        self.pressures = pressures
        self.arrhenius = arrhenius or []
        self.highPlimit = highPlimit or None

    def isPressureDependent(self):
        """
        Returns ``True`` since PDepArrhenius kinetics are pressure-dependent.
        """
        return True

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstring = reactantstr + arrow + productstr
        lines = ['pdep_arrhenius({0!r},'.format(rxnstring)]
        prefix = ' '*(indent+15)
        template = '[({0}, {1!r}), {2.A[0]:e}, {2.n}, {2.Ea[0]}],'
        for pressure,arrhenius in zip(self.pressures[0], self.arrhenius):
            lines.append(prefix + template.format(pressure,
                                                  self.pressures[1],
                                                  arrhenius))
        lines[-1] = lines[-1][:-1] + ')'
        return '\n'.join(lines)


class Chebyshev(KineticsModel):
    """
    A kinetic model of a phenomenological rate coefficient k(T, P) using the
    expression

    .. math:: \\log k(T,P) = \\sum_{t=1}^{N_T} \\sum_{p=1}^{N_P} \\alpha_{tp} \\phi_t(\\tilde{T}) \\phi_p(\\tilde{P})

    where :math:`\\alpha_{tp}` is a constant, :math:`\\phi_n(x)` is the
    Chebyshev polynomial of degree :math:`n` evaluated at :math:`x`, and

    .. math:: \\tilde{T} \\equiv \\frac{2T^{-1} - T_\\mathrm{min}^{-1} - T_\\mathrm{max}^{-1}}{T_\\mathrm{max}^{-1} - T_\\mathrm{min}^{-1}}

    .. math:: \\tilde{P} \\equiv \\frac{2 \\log P - \\log P_\\mathrm{min} - \\log P_\\mathrm{max}}{\\log P_\\mathrm{max} - \\log P_\\mathrm{min}}

    are reduced temperature and reduced pressures designed to map the ranges
    :math:`(T_\\mathrm{min}, T_\\mathrm{max})` and
    :math:`(P_\\mathrm{min}, P_\\mathrm{max})` to :math:`(-1, 1)`.
    The attributes are:

    =============== =============== ============================================
    Attribute       Type            Description
    =============== =============== ============================================
    `coeffs`        :class:`list`   Matrix of Chebyshev coefficients
    `kunits`        ``str``         The units of the generated k(T, P) values
    `degreeT`       :class:`int`    The number of terms in the inverse temperature direction
    `degreeP`       :class:`int`    The number of terms in the log pressure direction
    =============== =============== ============================================

    """

    def __init__(self, coeffs=None, kunits='', **kwargs):
        KineticsModel.__init__(self, **kwargs)
        if coeffs is not None:
            self.coeffs = np.array(coeffs, np.float64)
            self.degreeT = self.coeffs.shape[0]
            self.degreeP = self.coeffs.shape[1]
        else:
            self.coeffs = None
            self.degreeT = 0
            self.degreeP = 0
        self.kunits = kunits

    def isPressureDependent(self):
        """
        Returns ``True`` since Chebyshev polynomial kinetics are
        pressure-dependent.
        """
        return True

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstr = reactantstr + arrow + productstr
        prefix = ' '*(indent+19)
        lines = ['chebyshev_reaction({0!r},'.format(rxnstr),
                 prefix + 'Tmin={0.Tmin}, Tmax={0.Tmax},'.format(self),
                 prefix + 'Pmin={0.Pmin}, Pmax={0.Pmax},'.format(self)]
        for i in range(self.degreeT):
            coeffline = ', '.join('{0: 12.5e}'.format(self.coeffs[i,j]) for j in range(self.degreeP))
            if i == 0:
                lines.append(prefix + 'coeffs=[[{0}],'.format(coeffline))
            else:
                lines.append(prefix + '        [{0}],'.format(coeffline))

        lines[-1] = lines[-1][:-1] + '])'
        return '\n'.join(lines)


class ThirdBody(KineticsModel):
    """
    A kinetic model of a phenomenological rate coefficient k(T, P) using the
    expression

    .. math:: k(T,P) = k(T) [\\ce{M}]

    where :math:`k(T)` is an Arrhenius expression and
    :math:`[\\ce{M}] \\approx P/RT` is the concentration of the third body
    (i.e. the bath gas). A collision efficiency can be used to further correct
    the value of :math:`k(T,P)`.

    The attributes are:

    =============== ======================= ====================================
    Attribute       Type                    Description
    =============== ======================= ====================================
    `arrheniusHigh` :class:`Arrhenius`      The Arrhenius kinetics
    `efficiencies`  ``dict``                A mapping of species to collider efficiencies
    =============== ======================= ====================================

    """

    def __init__(self, arrheniusHigh=None, efficiencies=None, **kwargs):
        KineticsModel.__init__(self, **kwargs)
        self.arrheniusHigh = arrheniusHigh
        self.efficiencies = {}
        if efficiencies is not None:
            for mol, eff in efficiencies.iteritems():
                self.efficiencies[mol] = eff

    def isPressureDependent(self):
        """
        Returns ``True`` since third-body kinetics are pressure-dependent.
        """
        return True

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstr = reactantstr + ' + M' + arrow + productstr + ' + M'
        prefix = ' '*(indent + 20)
        lines = ['three_body_reaction({0!r}, {1},'.format(rxnstr, self.arrheniusHigh.rateStr())]
        if self.efficiencies:
            lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString()))

        lines[-1] = lines[-1][:-1] + ')'
        return '\n'.join(lines)


class Lindemann(ThirdBody):
    """
    A kinetic model of a phenomenological rate coefficient k(T, P) using the
    expression

    .. math:: k(T,P) = k_\\infty(T) \\left[ \\frac{P_\\mathrm{r}}{1 + P_\\mathrm{r}} \\right] F

    where

    .. math::

        P_\\mathrm{r} &= \\frac{k_0(T)}{k_\\infty(T)} [\\ce{M}]

        k_0(T) &= A_0 T^{n_0} \\exp \\left( - \\frac{E_0}{RT} \\right)

        k_\\infty(T) &= A_\\infty T^{n_\\infty} \\exp \\left( - \\frac{E_\\infty}{RT} \\right)

    and :math:`[\\ce{M}] \\approx P/RT` is the concentration of the
    bath gas. The Arrhenius expressions :math:`k_0(T)` and :math:`k_\\infty(T)`
    represent the low-pressure and high-pressure limit kinetics, respectively.
    The former is necessarily one reaction order higher than the latter. For
    the Lindemann model, :math:`F = 1`. A collision efficiency can be used to
    further correct the value of :math:`k(T,P)`.

    The attributes are:

    =============== ======================= ====================================
    Attribute       Type                    Description
    =============== ======================= ====================================
    `arrheniusLow`  :class:`Arrhenius`      The Arrhenius kinetics at the low-pressure limit
    `arrheniusHigh` :class:`Arrhenius`      The Arrhenius kinetics at the high-pressure limit
    `efficiencies`  ``dict``                A mapping of species to collider efficiencies
    =============== ======================= ====================================

    """

    def __init__(self, arrheniusLow=None, **kwargs):
        ThirdBody.__init__(self, **kwargs)
        self.arrheniusLow = arrheniusLow

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstr = reactantstr + arrow + productstr
        prefix = ' '*(indent + 17)
        lines = ['falloff_reaction({0!r},'.format(rxnstr)]
        lines.append(prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()))
        lines.append(prefix + 'kf0={0},'.format(self.arrheniusLow.rateStr()))
        if self.efficiencies:
            lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString()))

        lines[-1] = lines[-1][:-1] + ')'
        return '\n'.join(lines)


class Troe(Lindemann):
    """
    A kinetic model of a phenomenological rate coefficient k(T, P) using the
    expression

    .. math:: k(T,P) = k_\\infty(T) \\left[ \\frac{P_\\mathrm{r}}{1 + P_\\mathrm{r}} \\right] F

    where

    .. math::

        P_\\mathrm{r} &= \\frac{k_0(T)}{k_\\infty(T)} [\\ce{M}]

        k_0(T) &= A_0 T^{n_0} \\exp \\left( - \\frac{E_0}{RT} \\right)

        k_\\infty(T) &= A_\\infty T^{n_\\infty} \\exp \\left( - \\frac{E_\\infty}{RT} \\right)

    and :math:`[\\ce{M}] \\approx P/RT` is the concentration of the
    bath gas. The Arrhenius expressions :math:`k_0(T)` and :math:`k_\\infty(T)`
    represent the low-pressure and high-pressure limit kinetics, respectively.
    The former is necessarily one reaction order higher than the latter. A
    collision efficiency can be used to further correct the value of
    :math:`k(T,P)`.

    For the Troe model the parameter :math:`F` is computed via

    .. math::

        \\log F &= \\left\\{1 + \\left[ \\frac{\\log P_\\mathrm{r} + c}{n - d (\\log P_\\mathrm{r} + c)} \\right]^2 \\right\\}^{-1} \\log F_\\mathrm{cent}

        c &= -0.4 - 0.67 \\log F_\\mathrm{cent}

        n &= 0.75 - 1.27 \\log F_\\mathrm{cent}

        d &= 0.14

        F_\\mathrm{cent} &= (1 - \\alpha) \\exp \\left( -T/T_3 \\right) + \\alpha \\exp \\left( -T/T_1 \\right) + \\exp \\left( -T_2/T \\right)

    The attributes are:

    =============== ======================= ====================================
    Attribute       Type                    Description
    =============== ======================= ====================================
    `arrheniusLow`  :class:`Arrhenius`      The Arrhenius kinetics at the low-pressure limit
    `arrheniusHigh` :class:`Arrhenius`      The Arrhenius kinetics at the high-pressure limit
    `efficiencies`  ``dict``                A mapping of species to collider efficiencies
    `alpha`         :class:`Quantity`       The :math:`\\alpha` parameter
    `T1`            :class:`Quantity`       The :math:`T_1` parameter
    `T2`            :class:`Quantity`       The :math:`T_2` parameter
    `T3`            :class:`Quantity`       The :math:`T_3` parameter
    =============== ======================= ====================================

    """

    def __init__(self, alpha=0.0, T3=0.0, T1=0.0, T2=None, **kwargs):
        Lindemann.__init__(self, **kwargs)
        self.alpha = alpha
        self.T3 = T3
        self.T1 = T1
        self.T2 = T2

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstr = reactantstr + arrow + productstr
        prefix = ' '*17
        lines = ['falloff_reaction({0!r},'.format(rxnstr),
                 prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()),
                 prefix + 'kf0={0},'.format(self.arrheniusLow.rateStr())]

        if self.T2:
            troeArgs = 'A={0.alpha[0]}, T3={0.T3[0]}, T1={0.T1[0]}, T2={0.T2[0]}'.format(self)
        else:
            troeArgs = 'A={0.alpha[0]}, T3={0.T3[0]}, T1={0.T1[0]}'.format(self)
        lines.append(prefix + 'falloff=Troe({0}),'.format(troeArgs))

        if self.efficiencies:
            lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString()))

        # replace trailing comma
        lines[-1] = lines[-1][:-1] + ')'
        return '\n'.join(lines)


class Sri(Lindemann):
    """
    A kinetic model of a phenomenological rate coefficient k(T, P) using the
    "SRI" formulation of the blending function :math:`F` using either 3 or
    5 parameters. See :ref:`sec-sri-falloff`.

    The attributes are:

    =============== ======================= ====================================
    Attribute       Type                    Description
    =============== ======================= ====================================
    `arrheniusLow`  :class:`Arrhenius`      The Arrhenius kinetics at the low-pressure limit
    `arrheniusHigh` :class:`Arrhenius`      The Arrhenius kinetics at the high-pressure limit
    `efficiencies`  ``dict``                A mapping of species to collider efficiencies
    `A`            ``float``                The :math:`a` parameter
    `B`            ``float``                The :math:`b` parameter
    `C`            ``float``                The :math:`c` parameter
    `D`            ``float``                The :math:`d` parameter
    `E`            ``float``                The :math:`e` parameter
    =============== ======================= ====================================
    """

    def __init__(self, A=0.0, B=0.0, C=0.0, D=1.0, E=0.0, **kwargs):
        Lindemann.__init__(self, **kwargs)
        self.A = A
        self.B = B
        self.C = C
        self.D = D
        self.E = E

    def to_cti(self, reactantstr, arrow, productstr, indent=0):
        rxnstr = reactantstr + arrow + productstr
        prefix = ' '*17
        lines = ['falloff_reaction({0!r},'.format(rxnstr),
                 prefix + 'kf={0},'.format(self.arrheniusHigh.rateStr()),
                 prefix + 'kf0={0},'.format(self.arrheniusLow.rateStr())]

        if self.D == 1.0 and self.E == 0.0:
            sriArgs = 'A={0.A}, B={0.B}, C={0.C}'.format(self)
        else:
            sriArgs = 'A={0.A}, B={0.B}, C={0.C}, D={0.D}, E={0.E}'.format(self)
        lines.append(prefix + 'falloff=SRI({0}),'.format(sriArgs))

        if self.efficiencies:
            lines.append(prefix + 'efficiencies={0!r},'.format(self.efficiencyString()))

        # replace trailing comma
        lines[-1] = lines[-1][:-1] + ')'
        return '\n'.join(lines)


class TransportData(object):
    geometryFlags = ['atom', 'linear', 'nonlinear']

    def __init__(self, label, geometry, wellDepth, collisionDiameter,
                 dipoleMoment, polarizability, zRot, comment=None):

        if not isinstance(label, types.StringTypes):
            raise InputParseError("Bad label for transport data: " + repr(label))
        if int(geometry) not in (0,1,2):
            raise InputParseError("Bad geometry flag '{0}' for species '{1}'".format(geometry, label))

        self.label = label
        self.geometry = self.geometryFlags[int(geometry)]
        self.wellDepth = float(wellDepth)
        self.collisionDiameter = float(collisionDiameter)
        self.dipoleMoment = float(dipoleMoment)
        self.polarizability = float(polarizability)
        self.zRot = float(zRot)
        self.comment = comment or ''  # @todo: include this in the CTI

    def __repr__(self):
        return ('TransportData({label!r}, {geometry!r}, {wellDepth!r}, '
                '{collisionDiameter!r}, {dipoleMoment!r}, {polarizability!r}, '
                '{zRot!r}, {comment!r})').format(**self.__dict__)

    def to_cti(self, indent=0):
        prefix = ' '*(indent+18)
        lines = ['gas_transport(geom={0!r},'.format(self.geometry),
                 prefix+'diam={0},'.format(self.collisionDiameter),
                 prefix+'well_depth={0},'.format(self.wellDepth)]
        if self.dipoleMoment:
            lines.append(prefix+'dipole={0},'.format(self.dipoleMoment))
        if self.polarizability:
            lines.append(prefix+'polar={0},'.format(self.polarizability))
        if self.zRot:
            lines.append(prefix+'rot_relax={0},'.format(self.zRot))

        lines[-1] = lines[-1][:-1] + ')'
        return '\n'.join(lines)


def fortFloat(s):
    """
    Convert a string representation of a floating point value to a float,
    allowing for some of the peculiarities of allowable Fortran representations.
    """
    s = s.strip()
    s = s.replace('D', 'E').replace('d', 'e')
    s = s.replace('E ', 'E+').replace('e ', 'e+')
    return float(s)


class Parser(object):
    def __init__(self):
        self.processed_units = False
        self.energy_units = 'cal/mol'
        self.quantity_units = 'mol'
        self.warning_as_error = True

        self.elements = []
        self.speciesList = []
        self.speciesDict = {}
        self.reactions = []

    def warn(self, message):
        if self.warning_as_error:
            raise InputParseError(message)
        else:
            logging.warning(message)

    def parseComposition(self, elements, nElements, width):
        """
        Parse the elemental composition from a 7 or 9 coefficient NASA polynomial
        entry.
        """
        composition = {}
        for i in range(nElements):
            symbol = elements[width*i:width*i+2].strip()
            count = elements[width*i+2:width*i+width].strip()
            if not symbol:
                continue
            try:
                count = int(float(count))
                if count:
                    composition[symbol.capitalize()] = count
            except ValueError:
                pass
        return composition

    def getRateConstantUnits(self, length_dims, length_units, quantity_dims,
                             quantity_units, time_dims=1, time_units='s'):

        units = ''
        if length_dims:
            units += length_units
        if length_dims > 1:
            units += str(length_dims)
        if quantity_dims:
            units += '/' + quantity_units
        if quantity_dims > 1:
            units += str(quantity_dims)
        if time_dims:
            units += '/' + time_units
        if time_dims > 1:
            units += str(time_dims)
        if units.startswith('/'):
            units = '1' + units
        return units

    def readThermoEntry(self, entry, TintDefault):
        """
        Read a thermodynamics `entry` for one species in a Chemkin-format file
        (consisting of two 7-coefficient NASA polynomials). Returns the label of
        the species, the thermodynamics model as a :class:`MultiNASA` object, the
        elemental composition of the species, and the comment/note associated with
        the thermo entry.
        """
        lines = entry.splitlines()
        identifier = lines[0][0:24].split()
        species = identifier[0].strip()

        if len(identifier) > 1:
            note = ''.join(identifier[1:]).strip()
        else:
            note = ''

        # Extract the NASA polynomial coefficients
        # Remember that the high-T polynomial comes first!
        try:
            Tmin = fortFloat(lines[0][45:55])
            Tmax = fortFloat(lines[0][55:65])
            try:
                Tint = fortFloat(lines[0][65:75])
            except ValueError:
                Tint = TintDefault

            coeffs_high = [fortFloat(lines[i][j:k])
                           for i,j,k in [(1,0,15), (1,15,30), (1,30,45), (1,45,60),
                                         (1,60,75), (2,0,15), (2,15,30)]]
            coeffs_low = [fortFloat(lines[i][j:k])
                           for i,j,k in [(2,30,45), (2,45,60), (2,60,75), (3,0,15),
                                         (3,15,30), (3,30,45), (3,45,60)]]

        except (IndexError, ValueError):
            raise InputParseError('Error while reading thermo entry for species {0}'.format(species))

        composition = self.parseComposition(lines[0][24:44], 4, 5)

        # Non-standard extended elemental composition data may be located beyond
        # column 80 on the first line of the thermo entry
        if len(lines[0]) > 80:
            elements = lines[0][80:]
            composition2 = self.parseComposition(elements, len(elements)/10, 10)
            composition.update(composition2)

        # Construct and return the thermodynamics model
        thermo = MultiNASA(
            polynomials=[
                NASA(Tmin=(Tmin,"K"), Tmax=(Tint,"K"), coeffs=coeffs_low),
                NASA(Tmin=(Tint,"K"), Tmax=(Tmax,"K"), coeffs=coeffs_high)
            ],
            Tmin=(Tmin,"K"),
            Tmax=(Tmax,"K"),
        )

        return species, thermo, composition, note

    def readNasa9Entry(self, entry):
        """
        Read a thermodynamics `entry` for one species given as one or more
        9-coefficient NASA polynomials, written in the format described in
        Appendix A of NASA Reference Publication 1311 (McBride and Gordon, 1996).
        Returns the label of the species, the thermodynamics model as a
        :class:`MultiNASA` object, the elemental composition of the species, and
        the comment/note associated with the thermo entry.
        """
        tokens = entry[0].split()
        species = tokens[0]
        note = ' '.join(tokens[1:])
        N = int(entry[1][:2])
        note2 = entry[1][3:9].strip()
        if note and note2:
            note = '{0} [{1}]'.format(note, note2)
        elif note2:
            note = note2

        composition = self.parseComposition(entry[1][10:50], 5, 8)

        polys = []
        totalTmin = 1e100
        totalTmax = -1e100
        try:
            for i in range(N):
                A,B,C = entry[2+3*i:2+3*(i+1)]
                Tmin = fortFloat(A[1:11])
                Tmax = fortFloat(A[11:21])
                coeffs = [fortFloat(B[0:16]), fortFloat(B[16:32]),
                          fortFloat(B[32:48]), fortFloat(B[48:64]),
                          fortFloat(B[64:80]), fortFloat(C[0:16]),
                          fortFloat(C[16:32]), fortFloat(C[48:64]),
                          fortFloat(C[64:80])]
                polys.append(NASA(Tmin=(Tmin,"K"), Tmax=(Tmax,"K"), coeffs=coeffs))
                totalTmin = min(Tmin, totalTmin)
                totalTmax = max(Tmax, totalTmax)
        except (IndexError, ValueError):
            raise InputParseError('Error while reading thermo entry for species {0}'.format(species))

        thermo = MultiNASA(polynomials=polys,
                           Tmin=(totalTmin,"K"),
                           Tmax=(totalTmax,"K"))

        return species, thermo, composition, note

    def readKineticsEntry(self, entry):
        """
        Read a kinetics `entry` for a single reaction as loaded from a
        Chemkin-format file. Returns a :class:`Reaction` object with the
        reaction and its associated kinetics.
        """

        # Handle non-default units which apply to this entry
        energy_units = self.energy_units
        quantity_units = self.quantity_units
        if 'units' in entry.lower():
            for units in sorted(QUANTITY_UNITS, key=lambda k: -len(k)):
                m = re.search(r'units *\/ *%s *\/' % re.escape(units),
                              entry, re.IGNORECASE)
                if m:
                    entry = re.sub(r'units *\/ *%s *\/' % re.escape(units), '',
                                   entry, flags=re.IGNORECASE)
                    quantity_units = QUANTITY_UNITS[units]
                    break

            for units in sorted(ENERGY_UNITS, key=lambda k: -len(k)):
                m = re.search(r'units *\/ *%s *\/' % re.escape(units),
                              entry, re.IGNORECASE)
                if m:
                    entry = re.sub(r'units *\/ *%s *\/' % re.escape(units), '',
                                   entry, flags=re.IGNORECASE)
                    energy_units = ENERGY_UNITS[units]
                    break

        lines = entry.strip().splitlines()

        # The first line contains the reaction equation and a set of modified Arrhenius parameters
        tokens = lines[0].split()
        A = float(tokens[-3])
        n = float(tokens[-2])
        Ea = float(tokens[-1])
        reaction = ''.join(tokens[:-3])
        revReaction = None

        # Split the reaction equation into reactants and products
        if '<=>' in reaction:
            reversible = True
            reactants, products = reaction.split('<=>')
        elif '=>' in reaction:
            reversible = False
            reactants, products = reaction.split('=>')
        elif '=' in reaction:
            reversible = True
            reactants, products = reaction.split('=')
        else:
            raise InputParseError("Failed to find reactant/product delimiter in reaction string.")

        # Create a new Reaction object for this reaction
        reaction = Reaction(reactants=[], products=[], reversible=reversible)

        def parseExpression(expression, dest):
            falloff3b = None
            thirdBody = False  # simple third body reaction (non-falloff)

            # Look for third-body species for falloff reactions
            if re.search(r'\(\+[Mm]\)', expression):
                falloff3b = 'M'
                expression = re.sub(r'(\(\+[Mm]\))', '', expression)
            elif re.search(r'\(\+.*\)', expression):
                # See if it matches a known species
                for species in self.speciesDict:
                    if re.search(r'\(\+%s\)' % re.escape(species), expression):
                        falloff3b = species
                        expression = re.sub(r'(\(\+%s\))' % re.escape(species),
                                            '', expression)
                        break

            for term in expression.split('+'):
                term = term.strip()
                if not term[0].isalpha():
                    # This allows for for non-unity stoichiometric coefficients, e.g.
                    # 2A=B+C or .85A+.15B=>C
                    j = [i for i,c in enumerate(term) if c.isalpha()][0]
                    if term[:j].isdigit():
                        stoichiometry = int(term[:j])
                    else:
                        stoichiometry = float(term[:j])
                    species = term[j:]
                else:
                    species = term
                    stoichiometry = 1

                if species == 'M' or species == 'm':
                    thirdBody = True
                elif species not in self.speciesDict:
                    raise InputParseError('Unexpected species "{0}" in reaction expression "{1}".'.format(species, expression))
                else:
                    dest.append((stoichiometry, self.speciesDict[species]))

            return falloff3b, thirdBody

        falloff_3b_r, thirdBody = parseExpression(reactants, reaction.reactants)
        falloff_3b_p, thirdBody = parseExpression(products, reaction.products)

        if falloff_3b_r != falloff_3b_p:
            raise InputParseError('Third bodies do not match: "{0}" and "{1}" in'
                ' reaction entry:\n\n{2}'.format(falloff_3b_r, falloff_3b_p, entry))

        reaction.thirdBody = falloff_3b_r

        # Determine the appropriate units for k(T) and k(T,P) based on the number of reactants
        # This assumes elementary kinetics for all reactions
        rStoich = sum(r[0] for r in reaction.reactants) + (1 if thirdBody else 0)
        if rStoich > 3 or rStoich < 1:
            raise InputParseError('Invalid number of reactant species ({0}) for reaction {1}.'.format(rStoich, reaction))

        length_dim = 3 * (rStoich - 1)
        quantity_dim = rStoich - 1
        kunits = self.getRateConstantUnits(length_dim, 'cm',
                                           quantity_dim, quantity_units)
        klow_units = self.getRateConstantUnits(length_dim + 3, 'cm',
                                               quantity_dim + 1, quantity_units)

        # The rest of the first line contains the high-P limit Arrhenius parameters (if available)
        #tokens = lines[0][52:].split()
        tokens = lines[0].split()[1:]
        arrheniusHigh = Arrhenius(
            A=(A,kunits),
            n=n,
            Ea=(Ea, energy_units),
            T0=(1,"K"),
            parser=self
        )

        if len(lines) == 1:
            # If there's only one line then we know to use the high-P limit kinetics as-is
            reaction.kinetics = arrheniusHigh
        else:
            # There's more kinetics information to be read
            arrheniusLow = None
            troe = None
            sri = None
            chebyshev = None
            pdepArrhenius = None
            efficiencies = {}
            chebyshevCoeffs = []
            revReaction = None

            # Note that the subsequent lines could be in any order
            for line in lines[1:]:
                tokens = line.split('/')
                if 'dup' in line.lower():
                    # Duplicate reaction
                    reaction.duplicate = True

                elif 'low' in line.lower():
                    # Low-pressure-limit Arrhenius parameters
                    tokens = tokens[1].split()
                    arrheniusLow = Arrhenius(
                        A=(float(tokens[0].strip()),klow_units),
                        n=float(tokens[1].strip()),
                        Ea=(float(tokens[2].strip()),energy_units),
                        T0=(1,"K"),
                        parser=self
                    )

                elif 'rev' in line.lower():
                    reaction.reversible = False

                    # Create a reaction proceeding in the opposite direction
                    revReaction = Reaction(reactants=reaction.products,
                                           products=reaction.reactants,
                                           reversible=False)
                    tokens = tokens[1].split()
                    revReaction.kinetics = Arrhenius(
                        A=(float(tokens[0].strip()),klow_units),
                        n=float(tokens[1].strip()),
                        Ea=(float(tokens[2].strip()),energy_units),
                        T0=(1,"K"),
                        parser=self
                    )

                elif 'ford' in line.lower():
                    tokens = tokens[1].split()
                    reaction.fwdOrders[tokens[0].strip()] = tokens[1].strip()

                elif 'troe' in line.lower():
                    # Troe falloff parameters
                    tokens = tokens[1].split()
                    alpha = float(tokens[0].strip())
                    T3 = float(tokens[1].strip())
                    T1 = float(tokens[2].strip())
                    try:
                        T2 = float(tokens[3].strip())
                    except (IndexError, ValueError):
                        T2 = None

                    troe = Troe(
                        alpha=(alpha,''),
                        T3=(T3,"K"),
                        T1=(T1,"K"),
                        T2=(T2,"K") if T2 is not None else None,
                        parser=self
                    )
                elif 'sri' in line.lower():
                    # SRI falloff parameters
                    tokens = tokens[1].split()
                    A = float(tokens[0].strip())
                    B = float(tokens[1].strip())
                    C = float(tokens[2].strip())
                    try:
                        D = float(tokens[3].strip())
                        E = float(tokens[4].strip())
                    except (IndexError, ValueError):
                        D = None
                        E = None

                    if D is None or E is None:
                        sri = Sri(A=A, B=B, C=C, parser=self)
                    else:
                        sri = Sri(A=A, B=B, C=C, D=D, E=E, parser=self)

                elif 'cheb' in line.lower():
                    # Chebyshev parameters
                    if chebyshev is None:
                        chebyshev = Chebyshev()
                    tokens = [t.strip() for t in tokens]
                    if 'TCHEB' in line:
                        index = tokens.index('TCHEB')
                        tokens2 = tokens[index+1].split()
                        chebyshev.Tmin = float(tokens2[0].strip())
                        chebyshev.Tmax = float(tokens2[1].strip())
                    if 'PCHEB' in line:
                        index = tokens.index('PCHEB')
                        tokens2 = tokens[index+1].split()
                        chebyshev.Pmin = (float(tokens2[0].strip()), 'atm')
                        chebyshev.Pmax = (float(tokens2[1].strip()), 'atm')
                    if 'TCHEB' in line or 'PCHEB' in line:
                        pass
                    elif chebyshev.degreeT == 0 or chebyshev.degreeP == 0:
                        tokens2 = tokens[1].split()
                        chebyshev.degreeT = int(float(tokens2[0].strip()))
                        chebyshev.degreeP = int(float(tokens2[1].strip()))
                        chebyshev.coeffs = np.zeros((chebyshev.degreeT,chebyshev.degreeP), np.float64)
                    else:
                        tokens2 = tokens[1].split()
                        chebyshevCoeffs.extend([float(t.strip()) for t in tokens2])

                elif 'plog' in line.lower():
                    # Pressure-dependent Arrhenius parameters
                    if pdepArrhenius is None:
                        pdepArrhenius = []
                    tokens = tokens[1].split()
                    pdepArrhenius.append([float(tokens[0].strip()), Arrhenius(
                        A=(float(tokens[1].strip()),kunits),
                        n=float(tokens[2].strip()),
                        Ea=(float(tokens[3].strip()),energy_units),
                        T0=(1,"K"),
                        parser=self
                    )])
                else:
                    # Assume a list of collider efficiencies
                    for collider, efficiency in zip(tokens[0::2], tokens[1::2]):
                        efficiencies[collider.strip()] = float(efficiency.strip())

            if revReaction:
                revReaction.duplicate = reaction.duplicate

            # Decide which kinetics to keep and store them on the reaction object
            # Only one of these should be true at a time!
            if chebyshev is not None:
                if chebyshev.Tmin is None or chebyshev.Tmax is None:
                    raise InputParseError('Missing TCHEB line for reaction {0}'.format(reaction))
                if chebyshev.Pmin is None or chebyshev.Pmax is None:
                    raise InputParseError('Missing PCHEB line for reaction {0}'.format(reaction))
                index = 0
                for t in range(chebyshev.degreeT):
                    for p in range(chebyshev.degreeP):
                        chebyshev.coeffs[t,p] = chebyshevCoeffs[index]
                        index += 1
                reaction.kinetics = chebyshev
            elif pdepArrhenius is not None:
                reaction.kinetics = PDepArrhenius(
                    pressures=([P for P, arrh in pdepArrhenius],"atm"),
                    arrhenius=[arrh for P, arrh in pdepArrhenius],
                    parser=self
                )
            elif troe is not None:
                troe.arrheniusHigh = arrheniusHigh
                troe.arrheniusLow = arrheniusLow
                troe.efficiencies = efficiencies
                reaction.kinetics = troe
            elif sri is not None:
                sri.arrheniusHigh = arrheniusHigh
                sri.arrheniusLow = arrheniusLow
                sri.efficiencies = efficiencies
                reaction.kinetics = sri
            elif arrheniusLow is not None:
                reaction.kinetics = Lindemann(arrheniusHigh=arrheniusHigh,
                                              arrheniusLow=arrheniusLow,
                                              parser=self)
                reaction.kinetics.efficiencies = efficiencies
            elif thirdBody:
                reaction.kinetics = ThirdBody(arrheniusHigh=arrheniusHigh,
                                              parser=self)
                reaction.kinetics.efficiencies = efficiencies
            else:
                reaction.kinetics = arrheniusHigh

        return reaction, revReaction

    def loadChemkinFile(self, path):
        """
        Load a Chemkin-format input file to `path` on disk.
        """

        transportLines = []

        def removeCommentFromLine(line):
            if '!' in line:
                index = line.index('!')
                comment = line[index+1:-1]
                line = line[0:index] + '\n'
                return line, comment
            else:
                comment = ''
                return line, comment

        with open(path, 'r') as f:
            line = f.readline()
            while line != '':
                line = removeCommentFromLine(line)[0]
                line = line.strip()
                tokens = line.split()

                if 'ELEMENTS' in line:
                    index = tokens.index('ELEMENTS')
                    tokens = tokens[index+1:]
                    while 'END' not in tokens:
                        line = f.readline()
                        line = removeCommentFromLine(line)[0]
                        line = line.strip()
                        tokens.extend(line.split())

                    for token in tokens:
                        if token == 'END':
                            break
                        self.elements.append(token.capitalize())

                elif 'SPECIES' in line:
                    # List of species identifiers
                    index = tokens.index('SPECIES')
                    tokens = tokens[index+1:]
                    while 'END' not in tokens:
                        line = f.readline()
                        line = removeCommentFromLine(line)[0]
                        line = line.strip()
                        tokens.extend(line.split())

                    for token in tokens:
                        if token == 'END':
                            break
                        if token in self.speciesDict:
                            species = self.speciesDict[token]
                        else:
                            species = Species(label=token)
                            self.speciesDict[token] = species
                        self.speciesList.append(species)

                elif 'THERM' in line.upper() and 'NASA9' in line:
                    entryPosition = 0
                    entryLength = None
                    entry = []
                    while not line.startswith('END'):
                        line = f.readline()
                        line = removeCommentFromLine(line)[0]
                        if not line:
                            continue

                        if entryLength is None:
                            entryLength = 0
                            # special case if (redundant) temperature ranges are
                            # given as the first line
                            try:
                                s = line.split()
                                float(s[0]), float(s[1]), float(s[2])
                                continue
                            except (IndexError, ValueError):
                                pass

                        if entryPosition == 0:
                            entry.append(line)
                        elif entryPosition == 1:
                            entryLength = 2 + 3 * int(line.split()[0])
                            entry.append(line)
                        elif entryPosition < entryLength:
                            entry.append(line)

                        if entryPosition == entryLength-1:
                            label, thermo, comp, note = self.readNasa9Entry(entry)
                            try:
                                self.speciesDict[label].thermo = thermo
                                self.speciesDict[label].composition = comp
                                self.speciesDict[label].note = note
                            except KeyError:
                                logging.info('Skipping unexpected species "{0}" while reading thermodynamics entry.'.format(label))

                            entryPosition = -1
                            entry = []

                        entryPosition += 1

                elif 'THERM' in line:
                    # List of thermodynamics (hopefully one per species!)
                    line = f.readline()
                    TintDefault = float(line.split()[1])
                    thermo = ''
                    while line != '' and 'END' not in line:
                        line = removeCommentFromLine(line)[0]
                        if len(line) >= 80 and line[79] in ['1', '2', '3', '4']:
                            thermo += line
                            if line[79] == '4':
                                label, thermo, comp, note = self.readThermoEntry(thermo, TintDefault)
                                try:
                                    self.speciesDict[label].thermo = thermo
                                    self.speciesDict[label].composition = comp
                                    self.speciesDict[label].note = note
                                except KeyError:
                                    logging.info('Skipping unexpected species "{0}" while reading thermodynamics entry.'.format(label))
                                thermo = ''
                        line = f.readline()

                elif 'REACTIONS' in line:
                    # Reactions section
                    energyUnits = 'CAL/MOL'
                    moleculeUnits = 'MOLES'
                    try:
                        energyUnits = tokens[1]
                        moleculeUnits = tokens[2]
                    except IndexError:
                        pass

                    if not self.processed_units:
                        self.processed_units = True
                        self.energy_units = ENERGY_UNITS[energyUnits]
                        self.quantity_units = QUANTITY_UNITS[moleculeUnits]
                    else:
                        if (self.energy_units != ENERGY_UNITS[energyUnits] or
                            self.quantity_units != QUANTITY_UNITS[moleculeUnits]):
                            raise InputParseError("Multiple REACTIONS sections with "
                                                  "different units are not supported.")

                    kineticsList = []
                    commentsList = []
                    kinetics = ''
                    comments = ''

                    line = f.readline()
                    while line != '' and 'END' not in line:

                        lineStartsWithComment = line.startswith('!')
                        line, comment = removeCommentFromLine(line)
                        line = line.strip()
                        comment = comment.strip()

                        if '=' in line and not lineStartsWithComment:
                            # Finish previous record
                            kineticsList.append(kinetics)
                            commentsList.append(comments)
                            kinetics = ''
                            comments = ''

                        if line:
                            kinetics += line + '\n'
                        if comment:
                            comments += comment + '\n'

                        line = f.readline()

                    # Don't forget the last reaction!
                    if kinetics.strip() != '':
                        kineticsList.append(kinetics)
                        commentsList.append(comments)

                    if kineticsList[0] == '' and commentsList[-1] == '':
                        # True for mechanism files generated from RMG-Py
                        kineticsList.pop(0)
                        commentsList.pop(-1)
                    elif kineticsList[0] == '' and commentsList[0] == '':
                        # True for mechanism files generated from RMG-Java
                        kineticsList.pop(0)
                        commentsList.pop(0)
                    else:
                        # In reality, comments can occur anywhere in the mechanism
                        # file (e.g. either or both of before and after the
                        # reaction equation)
                        # If we can't tell what semantics we are using, then just
                        # throw the comments away
                        # (This is better than failing to load the mechanism file at
                        # all, which would likely occur otherwise)
                        if kineticsList[0] == '':
                            kineticsList.pop(0)
                        if len(kineticsList) != len(commentsList):
                            commentsList = ['' for kinetics in kineticsList]

                    for kinetics, comments in zip(kineticsList, commentsList):
                        reaction,revReaction = self.readKineticsEntry(kinetics)
                        self.reactions.append(reaction)
                        if revReaction is not None:
                            self.reactions.append(revReaction)

                elif 'TRAN' in line:
                    line = f.readline()
                    while 'END' not in line:
                        transportLines.append(line)

                line = f.readline()

        # Check for marked (and unmarked!) duplicate reactions
        # Raise exception for unmarked duplicate reactions
        for index1 in range(len(self.reactions)):
            reaction1 = self.reactions[index1]
            for index2 in range(index1+1, len(self.reactions)):
                reaction2 = self.reactions[index2]
                if reaction1.reactants == reaction2.reactants and reaction1.products == reaction2.products:
                    if reaction1.duplicate and reaction2.duplicate:
                        pass
                    elif reaction1.kinetics.isPressureDependent() == reaction2.kinetics.isPressureDependent():
                        # If both reactions are pressure-independent or both are pressure-dependent, then they need duplicate tags
                        # pdep and non-pdep reactions are treated as different, so those are okay
                        raise InputParseError('Encountered unmarked duplicate reaction {0}.'.format(reaction1))

        index = 0
        for reaction in self.reactions:
            index += 1
            reaction.index = index

        if transportLines:
            self.parseTransportData(transportLines)

    def parseTransportData(self, lines):
        """
        Parse the Chemkin-format transport data in ``lines`` (a list of strings)
        and add that transport data to the previously-loaded species.
        """

        for line in lines:
            line = line.strip()
            if not line or line.startswith('!'):
                continue
            if line.startswith('END'):
                break

            data = line.split()
            if len(data) < 7:
                raise InputParseError('Unable to parse transport data: not enough parameters')
            if len(data) >= 8:
                # comment may contain spaces. Rejoin into a single field.
                comment = ''.join(data[7:]).lstrip('!')
                data = data[:7] + [comment]

            speciesName = data[0]
            if speciesName in self.speciesDict:
                if self.speciesDict[speciesName].transport is None:
                    self.speciesDict[speciesName].transport = TransportData(*data)
                else:
                    self.warn('Ignoring duplicate transport data'
                         ' for species "{0}".'.format(speciesName))

    def writeCTI(self, header=None, name='gas', transportModel='Mix',
                 outName='mech.cti'):

        delimiterLine = '#' + '-'*79
        haveTransport = True
        speciesNameLength = 1
        elementsFromSpecies = set()
        for s in self.speciesList:
            if not s.transport:
                haveTransport = False
            if s.composition is None:
                raise InputParseError('No thermo data found for species: {0!r}'.format(s.label))
            elementsFromSpecies.update(s.composition)
            speciesNameLength = max(speciesNameLength, len(s.label))

        # validate list of elements
        missingElements = elementsFromSpecies - set(self.elements)
        if missingElements:
            raise InputParseError('Undefined elements: ' + str(missingElements))

        speciesNames = ['']
        for i,s in enumerate(self.speciesList):
            if i and not i % 5:
                speciesNames.append(' '*21)
            speciesNames[-1] += '{0:{1}s}'.format(s.label, speciesNameLength+2)

        speciesNames = '\n'.join(speciesNames).strip()

        lines = []
        if header:
            lines.extend(header)

        # Write the gas definition
        lines.append("units(length='cm', time='s', quantity={0!r}, act_energy={1!r})".format(self.quantity_units, self.energy_units))
        lines.append('')
        lines.append('ideal_gas(name={0!r},'.format(name))
        lines.append('          elements="{0}",'.format(' '.join(self.elements)))
        lines.append('          species="""{0}""",'.format(speciesNames))
        if self.reactions:
            lines.append("          reactions='all',")
        if haveTransport:
            lines.append("          transport={0!r},".format(transportModel))
        lines.append('          initial_state=state(temperature=300.0, pressure=OneAtm))')
        lines.append('')

        # Write the individual species data
        lines.append(delimiterLine)
        lines.append('# Species data')
        lines.append(delimiterLine)
        lines.append('')

        for s in self.speciesList:
            lines.append(s.to_cti())

        # Write the reactions
        lines.append(delimiterLine)
        lines.append('# Reaction data')
        lines.append(delimiterLine)

        for i,r in enumerate(self.reactions):
            lines.append('\n# Reaction {0}'.format(i+1))
            lines.append(r.to_cti())

        lines.append('')

        f = open(outName, 'w')
        f.write('\n'.join(lines))

    def showHelp(self):
        print """
ck2cti.py: Convert Chemkin-format mechanisms to Cantera input files (.cti)

Usage:
    ck2cti --input=<filename>
           [--thermo=<filename>]
           [--transport=<filename>]
           [--id=<phase-id>]
           [--output=<filename>]
           [--permissive]
           [-d | --debug]

Example:
    ck2cti --input=chem.inp --thermo=therm.dat --transport=tran.dat

If the output file name is not given, an output file with the same name as the
input file, with the extension changed to '.cti'.

The '--permissive' option allows certain recoverable parsing errors (e.g.
duplicate transport data) to be ignored.

"""

    def convertMech(self, inputFile, thermoFile=None,
                    transportFile=None, phaseName='gas',
                    outName=None, quiet=False, permissive=None):
        if quiet:
            logging.basicConfig(level=logging.ERROR)

        if permissive is not None:
            self.warning_as_error = not permissive

        # Read input mechanism files
        self.loadChemkinFile(inputFile)

        if thermoFile:
            self.loadChemkinFile(thermoFile)

        if transportFile:
            lines = open(transportFile).readlines()
            self.parseTransportData(lines)

            # Transport validation: make sure all species have transport data
            for s in self.speciesList:
                if s.transport is None:
                    raise InputParseError("No transport data for species '{0}'.".format(s))

        if not outName:
            outName = os.path.splitext(inputFile)[0] + '.cti'

        # Write output file
        self.writeCTI(name=phaseName, outName=outName)
        if not quiet:
            print 'Wrote CTI mechanism file to {0!r}.'.format(outName)
            print 'Mechanism contains {0} species and {1} reactions.'.format(len(self.speciesList), len(self.reactions))


if __name__ == '__main__':
    import getopt
    import sys

    longOptions = ['input=', 'thermo=', 'transport=', 'id=', 'output=',
                   'permissive', 'help', 'debug']

    try:
        optlist, args = getopt.getopt(sys.argv[1:], 'dh', longOptions)
        options = dict()
        for o,a in optlist:
            options[o] = a

        if args:
            raise getopt.GetoptError('Unexpected command line option: ' +
                                     repr(' '.join(args)))

    except getopt.GetoptError as e:
        print 'ck2cti.py: Error parsing arguments:'
        print e
        print 'Run "ck2cti.py --help" to see usage help.'
        sys.exit(1)

    parser = Parser()

    if not options or '-h' in options or '--help' in options:
        parser.showHelp()
        sys.exit(0)

    if '--input' in options:
        inputFile = options['--input']
    else:
        print 'Error: no mechanism input file specified'
        sys.exit(1)

    if '--output' in options:
        outName = options['--output']
        if not outName.endswith('.cti'):
            outName += '.cti'
    else:
        outName = None

    permissive = '--permissive' in options
    thermoFile = options.get('--thermo')
    transportFile = options.get('--transport')
    phaseName = options.get('--id', 'gas')

    parser.convertMech(inputFile, thermoFile, transportFile, phaseName,
                       outName, permissive=permissive)
