<filename>pyjac/core/mech_interpret.py
# -*- coding: utf-8 -*-
"""Chemkin-format mechanism interpreter module.
"""

# Python 2 compatibility
from __future__ import division

# Standard libraries
import sys
import math
import re
from copy import deepcopy
import logging

import numpy as np

# Local imports
from pyjac import utils
from pyjac.core import chem_model as chem

# Related module
CANTERA_FLAG = False
try:
    import cantera as ct
    version = ct.__version__.split('.')
    logger = logging.getLogger(__name__)
    if int(version[0]) < 2 or int(version[1]) < 3:
        logger.error('Parsing of Cantera mechanisms requires at least version 2.3.0'
                     'in order to access species thermo properties. Detected version'
                     ' is only {}'.format(ct.__version__))
        sys.exit(1)
    CANTERA_FLAG = True
except ImportError:
    CANTERA_FLAG = False

pre_units = ['moles', 'molecules']
"""list(`str`): Supported units list for pre-exponential factor"""

act_energy_units = ['kelvins', 'evolts', 'cal/mole', 'joules/kmole',
                    'kcal/mole', 'joules/mole', 'kjoules/mole'
                    ]
"""list(`str`): Supported units list for activation energy"""

act_energy_fact = dict({'kelvins': 1.0,
                        'evolts': 11595.,
                        'cal/mole': 4.184 / chem.RU_JOUL,
                        'kcal/mole': 4184. / chem.RU_JOUL,
                        'joules/mole': 1. / chem.RU_JOUL,
                        'kjoules/mole': 1000.0 / chem.RU_JOUL,
                        'joules/kmole': 1. / (chem.RU_JOUL * 1000.)
                        })
"""dict: Activation energy conversion factor"""

# get local element atomic weight dict
elem_wt = chem.get_elem_wt()


def sort_reactions(reacs, sort_type, return_order=False):
    """
    reacs : list of `ReacInfo`
        List of reactions in mechanism to be sorted.
    sort_type : :class:`enum_types.reaction_sorting`
        If not None, sort the reactions in this mechanism according to the supplied
        scheme
    return_order : bool [False]
        If True, return new order of sorted reactions, if false, return the sorted
        reactions themselves

    """

    ordering = np.arange(len(reacs))
    from pyjac.core.enum_types import reaction_sorting
    if sort_type is not None and sort_type != reaction_sorting.none:
        from pyjac.core.enum_types import (
            reaction_type, falloff_form, reversible_type, thd_body_type)

        # we consider the following enums:
        # column #1, reaction type
        # column #2, falloff form
        # column #3, third body form
        # column #4, reversible type

        enum_order = (reaction_type, falloff_form, thd_body_type, reversible_type)
        sort_matrix = np.empty((len(reacs), 4))
        for i, rxn in enumerate(reacs):
            for j, enum in enumerate(enum_order):
                e_val = rxn.get_type(enum)
                assert len(e_val) == 1
                sort_matrix[i, j] = int(e_val[0])

        ordering = np.arange(len(reacs))

        # now sort by column -- https://stackoverflow.com/a/38194077
        for i in reversed(range(len(enum_order))):
            inds = sort_matrix[:, i].argsort(kind='mergesort')
            ordering = ordering[inds]
            sort_matrix = sort_matrix[inds]

    if return_order:
        return ordering

    # sort reactions
    return [reacs[i] for i in ordering]


def read_mech(mech_filename, therm_filename, sort_type=None):
    """Read and interpret mechanism file for elements, species, and reactions.

    Parameters
    ----------
    mech_filename : str
        Reaction mechanism filename (e.g. 'mech.dat')
    therm_filename : str, optional
        Thermodynamic database filename (e.g., 'therm.dat')
    sort_type : :class:`enum_types.reaction_sorting`
        If not None, sort the reactions in this mechanism according to the supplied
        scheme

    Returns
    -------
    elems : list of str
        List of elements in mechanism.
    specs : list of `SpecInfo`
        List of species in mechanism.
    reacs : list of `ReacInfo`
        List of reactions in mechanism.
    units : str
        Units of reactions' Arrhenius coefficients

    Notes
    -----
    Doesn't support element names with digits.

    """

    elems = []
    reacs = []
    specs = []
    key = ''

    with open(mech_filename, 'r') as file:
        # start line reading loop
        while True:
            line = file.readline()

            # end of file
            if not line:
                break

            # skip blank or commented lines
            if re.search(r'^\s*$', line) or re.search(r'^\s*!', line):
                continue

            # don't convert to lowercase, since thermo
            # needs to match (for Chemkin)

            # Remove trailing and leading whitespace, tabs, newline
            line = line.strip()

            # remove any comments from end of line
            ind = line.find('!')
            if ind > 0:
                line = line[0:ind]

            # now determine key
            if line[0:4].lower() == 'elem':
                key = 'elem'

                # check for any entries on this line
                line_split = line.split()
                if len(line_split) > 1:
                    ind = line.index(line_split[1])
                    line = line[ind:]
                else:
                    continue

            elif line[0:4].lower() == 'spec':
                key = 'spec'

                # check for any entries on this line
                line_split = line.split()
                if len(line_split) > 1:
                    ind = line.index(line_split[1])
                    line = line[ind:]
                else:
                    continue

            elif line[0:4].lower() == 'reac':
                key = 'reac'

                # default units from Chemkin
                units_E = 'cal/mole'
                units_A = 'moles'

                # get Arrhenius coefficient units (if specified)
                for unit in line.split()[1:]:
                    if unit.lower() in pre_units:
                        units_A = unit.lower()
                    elif unit.lower() in act_energy_units:
                        units_E = unit.lower()
                    else:
                        logger = logging.getLogger(__name__)
                        logger.error('Unsupported units on REACTION line. For '
                                     'pre-exponential factor, choose from: ' +
                                     ', '.join(pre_units) + '\nFor activation energy'
                                     ', choose from: ' + ', '.join(act_energy_units)
                                     + '\nOtherwise leave blank for moles and '
                                     'cal/mole.')
                        sys.exit(1)

                if units_A == 'molecules':
                    raise NotImplementedError('Molecules units not '
                                              'supported, sorry.'
                                              )
                continue
            elif line[0:4].lower() == 'ther':
                # thermo data is in mechanism file
                read_thermo(mech_filename, elems, specs)
                continue
            elif line[0:3].lower() == 'end':
                key = ''
                continue

            if key == 'elem':
                # if any atomic weight declarations, replace / with spaces
                line = line.replace('/', ' ')

                line_split = line.split()
                e_last = ''
                for e in line_split:
                    if e.isalpha():
                        if e[0:3] == 'end':
                            continue
                        if e not in elems:
                            elems.append(e)
                        e_last = e
                    else:
                        # either add new element or update existing
                        # atomic weight
                        elem_wt[e_last.lower()] = float(e)

            elif key == 'spec':
                line_split = line.split()
                for s in line_split:
                    if s[0:3] == 'end':
                        continue
                    if not next((sp for sp in specs if sp.name == s), None):
                        specs.append(chem.SpecInfo(s))

            elif key == 'reac':
                # determine if reaction or auxiliary info line

                if '=' in line:
                    # new reaction

                    cheb_flag = False

                    # get Arrhenius coefficients
                    line_split = line.split()
                    n = len(line_split)
                    reac_A = float(line_split[n - 3])
                    reac_b = float(line_split[n - 2])
                    reac_E = float(line_split[n - 1])

                    ind = line.index(line_split[n - 3])
                    line = line[0:ind].strip()

                    if '<=>' in line:
                        ind = line.index('<=>')
                        reac_rev = True
                        reac_str = line[0:ind].strip()
                        prod_str = line[ind + 3:].strip()
                    elif '=>' in line:
                        ind = line.index('=>')
                        reac_rev = False
                        reac_str = line[0:ind].strip()
                        prod_str = line[ind + 2:].strip()
                    else:
                        ind = line.index('=')
                        reac_rev = True
                        reac_str = line[0:ind].strip()
                        prod_str = line[ind + 1:].strip()

                    thd = False
                    pdep = False
                    pdep_sp = ''

                    reac_spec = []
                    reac_nu = []
                    prod_spec = []
                    prod_nu = []

                    # reactants

                    # look for third-body species
                    sub_str = reac_str
                    while '(' in sub_str:
                        ind1 = sub_str.find('(')
                        ind2 = sub_str.find(')')

                        # Need to check if '+' is first character inside
                        # parentheses and not embedded within parentheses
                        # (e.g., '(+)').
                        # If not, part of species name.
                        inParen = sub_str[ind1 + 1: ind2].strip()
                        if inParen is '+':
                            # '+' embedded within parentheses
                            sub_str = sub_str[ind2 + 1:]
                        elif inParen[0] is '+':
                            pdep = True

                            # either 'm' or a specific species
                            pdep_sp = sub_str[ind1 + 1: ind2].replace('+', ' ')
                            pdep_sp = pdep_sp.strip()

                            if pdep_sp.lower() == 'm':
                                thd = True
                                pdep_sp = ''

                            # now remove from string
                            ind = reac_str.find(sub_str)
                            reac_str = (reac_str[0: ind1 + ind] +
                                        reac_str[ind2 + ind + 1:]
                                        )
                            break
                        else:
                            # Part of species name, remove from substring
                            # and look at rest of reactant line.
                            sub_str = sub_str[ind2 + 1:]

                    reac_list = reac_str.split('+')

                    # Check for empty list elements, meaning there were
                    # multiple '+' in a row, which indicates species'
                    # name ended in '+'.
                    while '' in reac_list:
                        ind = reac_list.index('')
                        reac_list[ind - 1] = reac_list[ind - 1] + '+'
                        del reac_list[ind]

                    # check for any species with '(+)' that was split apart
                    for sp in reac_list:
                        ind = reac_list.index(sp)

                        # ensure not last entry
                        if (ind < len(reac_list) - 1):
                            spNext = reac_list[ind + 1]
                            if sp[len(sp) - 1] is '(' and spNext[0] is ')':
                                reac_list[ind] = sp + '+' + spNext
                                del reac_list[ind + 1]

                    for sp in reac_list:

                        sp = sp.strip()

                        # look for coefficient
                        if sp[0:1].isdigit():
                            # starts with number (coefficient)

                            # search for first letter
                            for i in range(len(sp)):
                                if sp[i: i + 1].isalpha():
                                    break

                            nu = sp[0:i]
                            if '.' in nu:
                                # float
                                nu = float(nu)
                            else:
                                # integer
                                nu = int(nu)

                            sp = sp[i:].strip()
                        else:
                            # no coefficient given
                            nu = 1

                        # check for third body
                        if sp.lower() == 'm':
                            thd = True
                            continue

                        # check if species already in reaction
                        if sp not in reac_spec:
                            # new reactant
                            reac_spec.append(sp)
                            reac_nu.append(nu)
                        else:
                            # existing reactant
                            i = reac_spec.index(sp)
                            reac_nu[i] += nu

                    # products

                    # look for third-body species
                    sub_str = prod_str
                    while '(' in sub_str:
                        ind1 = sub_str.find('(')
                        ind2 = sub_str.find(')')

                        # Need to check if '+' is first character inside
                        # parentheses and not embedded within parentheses
                        # (e.g., '(+)'). If not, part of species name.
                        inParen = sub_str[ind1 + 1: ind2].strip()
                        if inParen is '+':
                            # '+' embedded within parentheses
                            sub_str = sub_str[ind2 + 1:]
                        elif inParen[0] is '+':
                            pdep = True

                            # either 'm' or a specific species
                            pdep_sp = sub_str[ind1 + 1: ind2].replace('+', ' ')
                            pdep_sp = pdep_sp.strip()

                            if pdep_sp.lower() == 'm':
                                thd = True
                                pdep_sp = ''

                            # now remove from string
                            ind = prod_str.find(sub_str)
                            prod_str = (prod_str[0: ind1 + ind] +
                                        prod_str[ind2 + ind + 1:]
                                        )
                            break
                        else:
                            # Part of species name, remove from substring and
                            # look at rest of product line.
                            sub_str = sub_str[ind2 + 1:]

                    prod_list = prod_str.split('+')

                    # Check for empty list elements, meaning there were
                    # multiple '+' in a row, which indicates species
                    # name ended in '+'.
                    while '' in prod_list:
                        ind = prod_list.index('')
                        prod_list[ind - 1] = prod_list[ind - 1] + '+'
                        del prod_list[ind]

                    # check for any species with '(+)' that was split apart
                    for sp in prod_list:
                        ind = prod_list.index(sp)

                        # ensure not last entry
                        if (ind < len(prod_list) - 1):
                            spNext = prod_list[ind + 1]
                            if sp[len(sp) - 1] is '(' and spNext[0] is ')':
                                prod_list[ind] = sp + '+' + spNext
                                del prod_list[ind + 1]

                    for sp in prod_list:

                        sp = sp.strip()

                        # look for coefficient
                        if sp[0:1].isdigit():
                            # starts with number (coefficient)

                            # search for first letter
                            for i in range(len(sp)):
                                if sp[i: i + 1].isalpha():
                                    break

                            nu = sp[0:i]
                            if '.' in nu:
                                # float
                                nu = float(nu)
                            else:
                                # integer
                                nu = int(nu)

                            sp = sp[i:].strip()
                        else:
                            # no coefficient given
                            nu = 1

                        # check for third body
                        if sp in ['m', 'M']:
                            thd = True
                            continue

                        # check if species already in reaction
                        if sp not in prod_spec:
                            # new product
                            prod_spec.append(sp)
                            prod_nu.append(nu)
                        else:
                            # existing product
                            i = prod_spec.index(sp)
                            prod_nu[i] += nu

                    # Don't want to confuse third-body and pressure-dependent
                    # reactions... they are different!
                    if pdep:
                        thd = False

                    # Convert given activation energy units to internal units
                    reac_E *= act_energy_fact[units_E]

                    # Convert given pre-exponential units to internal units
                    if units_A == 'moles':
                        reac_ord = sum(reac_nu)
                        if thd:
                            reac_A /= 1000. ** reac_ord
                        elif pdep:
                            # Low- (chemically activated bimolecular reaction) or
                            # high-pressure (fall-off reaction) limit
                            # parameters
                            reac_A /= 1000. ** (reac_ord - 1.)
                        else:
                            # Elementary reaction
                            reac_A /= 1000. ** (reac_ord - 1.)

                    # add reaction to list
                    reac = chem.ReacInfo(reac_rev, reac_spec, reac_nu,
                                         prod_spec, prod_nu, reac_A, reac_b,
                                         reac_E
                                         )
                    reac.thd_body = thd
                    reac.pdep = pdep
                    if pdep:
                        reac.pdep_sp = pdep_sp

                    reacs.append(reac)

                else:
                    # auxiliary reaction info

                    aux = line[0:3].lower()
                    if aux == 'dup':
                        reacs[-1].dup = True

                    elif aux == 'rev':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        # Convert reverse activation energy units
                        par3 *= act_energy_fact[units_E]

                        # Convert reverse pre-exponential factor
                        if units_A == 'moles':
                            reac_ord = sum(reacs[-1].prod_nu)
                            if thd:
                                par1 /= 1000. ** reac_ord
                            elif pdep:
                                # Low- (chemically activated bimolecular reaction) or
                                # high-pressure (fall-off reaction) limit
                                # parameters
                                par1 /= 1000. ** (reac_ord - 1.)
                            else:
                                # Elementary reaction
                                par1 /= 1000. ** (reac_ord - 1.)

                        # Ensure nonzero reverse coefficients
                        if par1 != 0.0:
                            reacs[-1].rev_par.append(par1)
                            reacs[-1].rev_par.append(par2)
                            reacs[-1].rev_par.append(par3)
                        else:
                            reacs[-1].rev = False

                    elif aux == 'low':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        # Convert low-pressure activation energy units
                        par3 *= act_energy_fact[units_E]

                        # Convert low-pressure pre-exponential factor
                        if units_A == 'moles':
                            par1 /= 1000. ** sum(reacs[-1].reac_nu)

                        reacs[-1].low.append(par1)
                        reacs[-1].low.append(par2)
                        reacs[-1].low.append(par3)

                    elif aux == 'hig':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        # Convert high-pressure activation energy units
                        par3 *= act_energy_fact[units_E]

                        # Convert high-pressure pre-exponential factor
                        if units_A == 'moles':
                            par1 /= 1000. ** (sum(reacs[-1].reac_nu) - 2.)

                        reacs[-1].high.append(par1)
                        reacs[-1].high.append(par2)
                        reacs[-1].high.append(par3)

                    elif aux == 'tro':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        reacs[-1].troe = True
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        do_warn = False
                        if par2 == 0:
                            do_warn = True
                            par2 = 1e-30
                        if par3 == 0:
                            do_warn = True
                            par3 = 1e-30
                        if do_warn:
                            logger = logging.getLogger(__name__)
                            logger.warn(
                                'Troe parameters in reaction {} modified to avoid'
                                ' division by zero!.'.format(len(reacs)))

                        reacs[-1].troe_par.append(par1)
                        reacs[-1].troe_par.append(par2)
                        reacs[-1].troe_par.append(par3)

                        # optional fourth parameter
                        if len(line_split) > 4:
                            par4 = float(line_split[4])
                            reacs[-1].troe_par.append(par4)

                    elif aux == 'sri':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        reacs[-1].sri = True
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])
                        reacs[-1].sri_par.append(par1)
                        reacs[-1].sri_par.append(par2)
                        reacs[-1].sri_par.append(par3)

                        # optional fourth and fifth parameters
                        if len(line_split) > 4:
                            par4 = float(line_split[4])
                            par5 = float(line_split[5])
                            reacs[-1].sri_par.append(par4)
                            reacs[-1].sri_par.append(par5)
                    elif aux == 'che':
                        reacs[-1].cheb = True
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        if cheb_flag:
                            for par in line_split[1:]:
                                reacs[-1].cheb_par.append(float(par))
                        else:
                            # first CHEB line
                            cheb_flag = True
                            # Don't want Cheb reactions lumped in with
                            # standard falloff.
                            reacs[-1].pdep = False
                            reacs[-1].cheb_n_temp = int(line_split[1])
                            reacs[-1].cheb_n_pres = int(line_split[2])
                            reacs[-1].cheb_par = []
                            for par in line_split[3:]:
                                reacs[-1].cheb_par.append(float(par))
                    elif aux == 'pch':
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        # Convert pressure from atm to Pa
                        reacs[-1].cheb_plim = [float(line_split[1]) * chem.PA,
                                               float(line_split[2]) * chem.PA
                                               ]

                        # Look for temperature limits on same line:
                        if line_split[3].lower() == 'tcheb':
                            reacs[-1].cheb_tlim = [float(line_split[4]),
                                                   float(line_split[5])
                                                   ]
                    elif aux == 'tch':
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        reacs[-1].cheb_tlim = [float(line_split[1]),
                                               float(line_split[2])
                                               ]
                        # Look for pressure limits on same line:
                        if line_split[3].lower() == 'pcheb':
                            reacs[-1].cheb_plim = [float(line_split[4]) * chem.PA,
                                                   float(
                                                       line_split[5]) * chem.PA
                                                   ]
                    elif aux == 'plo':
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        if not reacs[-1].plog:
                            reacs[-1].plog = True
                            # Don't want Plog reactions lumped in with
                            # standard falloff.
                            reacs[-1].pdep = False
                            reacs[-1].plog_par = []
                        pars = [float(n) for n in line_split[1:5]]

                        # Convert pressure from atm to Pa
                        pars[0] *= 101325.0

                        # Convert given activation energy units to internal
                        # units
                        pars[3] *= act_energy_fact[units_E]

                        # Convert given pre-exponential units to internal units
                        if units_A == 'moles':
                            reac_ord = sum(reacs[-1].reac_nu)
                            # Looks like elementary reaction
                            pars[1] /= 1000. ** (reac_ord - 1.)

                        reacs[-1].plog_par.append(pars)
                    else:
                        # enhanced third body efficiencies
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        for i in range(0, len(line_split), 2):
                            pair = [line_split[i], float(line_split[i + 1])]
                            reacs[-1].thd_body_eff.append(pair)

    # process some reaction auxiliary info
    for idx, reac in enumerate(reacs):
        if reac.cheb:
            # check for correct number
            n = reac.cheb_n_temp
            m = reac.cheb_n_pres
            if len(reac.cheb_par) != n * m:
                logger = logging.getLogger(__name__)
                logging.error('Incorrect number of CHEB coefficients in reaction '
                              + repr(idx))
                sys.exit(1)
            else:
                # Convert units of first Chebyshev parameter
                order = sum(reac.reac_nu)
                if units_A == 'moles':
                    reac.cheb_par[0] += math.log10(0.001 ** (order - 1.))

                reacs[idx].cheb_par = np.reshape(reac.cheb_par, (n, m))

    # check that all species in reactions correspond to a known species
    spec_names = set(spec.name for spec in specs)
    for idx, reac in enumerate(reacs):
        in_rxn = set(reac.reac + reac.prod)
        for spec in in_rxn:
            if spec not in spec_names:
                logger = logging.getLogger(__name__)
                logger.error('Reaction {} contains unknown species {}'.format(
                    idx, spec))
                sys.exit(-1)

    # Split reversible reactions with explicit reverse parameters into
    # two irreversible reactions to match Cantera's behavior
    for reac in reacs[:]:
        if reac.rev_par:
            new_reac = deepcopy(reac)

            idx = reacs.index(reac)
            reacs[idx].rev = False
            reacs[idx].rev_par = []

            new_reac.A = new_reac.rev_par[0]
            new_reac.b = new_reac.rev_par[1]
            new_reac.E = new_reac.rev_par[2]
            new_reac.rev = False
            new_reac.rev_par = []
            new_reac.prod = reac.reac[:]
            new_reac.prod_nu = reac.reac_nu[:]
            new_reac.reac = reac.prod[:]
            new_reac.reac_nu = reac.prod_nu[:]

            reacs.insert(idx + 1, new_reac)

    # Read seperate thermo file if present and needed
    if any([not spec.mw for spec in specs]):
        if therm_filename:
            read_thermo(therm_filename, elems, specs)
        else:
            logger = logging.getLogger(__name__)
            logger.error(
                'No thermo file specified, but species missing \n'
                'data. Either specify file, or ensure complete data in\n'
                'mechanism file with THERMO option.')
            sys.exit(1)

    # Check for missing thermo data again
    missing_mw = [spec.name for spec in specs if not spec.mw]
    if missing_mw:
        logger = logging.getLogger(__name__)
        logger.error('Missing thermo data for ' + ', '.join(missing_mw))
        sys.exit(1)

    # determine reaction type enums
    for reac in reacs:
        reac.finalize(len(specs))
    for spec in specs:
        spec.finalize()

    if sort_type:
        reacs = sort_reactions(reacs, sort_type)

    return (elems, specs, reacs)


def read_thermo(filename, elems, specs):
    """Read and interpret thermodynamic database for species data.

    Reads the thermodynamic file and returns the species thermodynamic
    coefficients as well as the species-specific temperature range
    values (if given).

    Parameters
    ----------
    filename : str
        Name of thermo database file.
    elems : list of str
        List of element names in mechanism.
    specs : list of `SpecInfo`
        List of species in mechanism.

    Returns
    -------
    None

    """

    # choose our element names to enable comparison between Chemkin & Cantera
    elem_name_map = {e.lower(): e for e in elem_wt}

    with open(filename, 'r') as file:

        # loop through intro lines
        while True:
            line = file.readline()

            # skip blank or commented lines
            if re.search(r'^\s*$', line) or re.search(r'^\s*!', line):
                continue

            # skip 'thermo' at beginning
            if 'thermo' in line.lower():
                break

        # next line either has common temperature ranges or first species
        last_line = file.tell()
        line = file.readline()

        line_split = line.split()
        if line_split[0][0:1].isdigit():
            T_ranges = utils.read_str_num(line)
        else:
            # no common temperature info
            file.seek(last_line)
            # default

        # now start reading species thermo info
        while True:
            # first line of species info
            line = file.readline()

            # don't convert to lowercase, needs to match thermo for Chemkin

            # break if end of file
            if line is None or line[0:3].lower() == 'end':
                break

            # skip blank/commented line
            if re.search(r'^\s*$', line) or re.search(r'^\s*!', line):
                continue

            # species name, columns 0:18
            spec = line[0:18].strip()

            # Apparently, in some cases, notes are in the
            # columns of shorter species names, so make
            # sure no spaces.
            if spec.find(' ') > 0:
                spec = spec[0: spec.find(' ')]

            # now need to determine if this species is in mechanism
            if next((sp for sp in specs if sp.name == spec), None):
                sp_ind = next(i for i in range(len(specs))
                              if specs[i].name == spec
                              )
            else:
                # not in mechanism, read next three lines and continue
                line = file.readline()
                line = file.readline()
                line = file.readline()
                continue

            # set species to the one matched
            spec = specs[sp_ind]

            # ensure not reading the same species more than once...
            if spec.mw:
                # already done! skip next three lines
                line = file.readline()
                line = file.readline()
                line = file.readline()
                continue

            # now get element composition of species, columns 24:44
            # each piece of data is 5 characters long (2 for element, 3 for #)
            elem_str = utils.split_str(line[24:44], 5)

            for e_str in elem_str:
                e = e_str[0:2].strip()
                # skip if blank
                if e == '' or e == '0':
                    continue
                # may need to convert to float first, in case of e.g. "1."
                e_num = float(e_str[2:].strip())
                e_num = int(e_num)

                spec.elem.append([elem_name_map[e.lower()], e_num])

                # calculate molecular weight
                spec.mw += e_num * elem_wt[e.lower()]

            # temperatures for species
            T_spec = utils.read_str_num(line[45:74])
            T_low = T_spec[0]
            T_high = T_spec[1]
            if len(T_spec) == 3:
                T_com = T_spec[2]
            else:
                T_com = T_ranges[1]

            spec.Trange = [T_low, T_com, T_high]

            # second species line
            line = file.readline()
            coeffs = utils.split_str(line[0:75], 15)
            spec.hi[0] = float(coeffs[0])
            spec.hi[1] = float(coeffs[1])
            spec.hi[2] = float(coeffs[2])
            spec.hi[3] = float(coeffs[3])
            spec.hi[4] = float(coeffs[4])

            # third species line
            line = file.readline()
            coeffs = utils.split_str(line[0:75], 15)
            spec.hi[5] = float(coeffs[0])
            spec.hi[6] = float(coeffs[1])
            spec.lo[0] = float(coeffs[2])
            spec.lo[1] = float(coeffs[3])
            spec.lo[2] = float(coeffs[4])

            # fourth species line
            line = file.readline()
            coeffs = utils.split_str(line[0:75], 15)
            spec.lo[3] = float(coeffs[0])
            spec.lo[4] = float(coeffs[1])
            spec.lo[5] = float(coeffs[2])
            spec.lo[6] = float(coeffs[3])

            # stop reading if all species in mechanism accounted for
            if not next((sp for sp in specs if sp.mw == 0.0), None):
                break

    return None


def read_mech_ct(filename=None, gas=None, sort_type=None):
    """Read and interpret Cantera-format mechanism file.

    Parameters
    ----------
    filename : str
        Reaction mechanism filename (e.g. 'mech.cti'). Optional.
    gas : `cantera.Solution` object
        Existing Cantera Solution object to be used. Optional.
    sort_type : :class:`enum_types.reaction_sorting`
        If not None, sort the reactions in this mechanism according to the supplied
        scheme

    Returns
    -------
    elems : list of str
        List of elements in mechanism.
    specs : list of `SpecInfo`
        List of species in mechanism.
    reacs : list of `ReacInfo`
        List of reactions in mechanism.
    units : str
        Units of reactions' Arrhenius coefficients

    """

    if not CANTERA_FLAG:
        logger = logging.getLogger(__name__)
        logger.error('Cantera not installed. Cannot interpret Cantera-format'
                     ' mechanism.')
        sys.exit(1)

    if filename:
        gas = ct.Solution(filename)
    elif not gas:
        logger = logging.getLogger(__name__)
        logger.error('Need either filename or Cantera Solution object.')
        sys.exit(1)

    # Elements
    elems = gas.element_names
    for e, wt in zip(elems, gas.atomic_weights):
        if e.lower() not in elem_wt:
            elem_wt[e.lower()] = wt

    # choose our element names to enable comparison between Chemkin & Cantera
    elem_name_map = {e.lower(): e for e in elem_wt}

    # Species
    specs = []
    for i, sp in enumerate(gas.species_names):
        # Get Species object
        species = gas.species(i)

        if isinstance(species.thermo, ct.NasaPoly2):
            spec = chem.SpecInfo(sp)
        elif isinstance(species.thermo, ct.Nasa9PolyMultiTempRegion):
            spec = chem.SpecInfoN9(sp)

        spec.mw = gas.molecular_weights[i]

        # Species elemental composition
        for e in species.composition:
            spec.elem.append([elem_name_map[e.lower()], int(species.composition[e])])

        # Species thermodynamic properties
        coeffs = species.thermo.coeffs
        if isinstance(species.thermo, ct.NasaPoly2):
            spec.Trange = [species.thermo.min_temp, coeffs[0],
                            species.thermo.max_temp
                            ]
            spec.hi = coeffs[1:8]
            spec.lo = coeffs[8:15]
        elif isinstance(species.thermo, ct.Nasa9PolyMultiTempRegion):
            spec.nzones = int(coeffs[0])
            spec.lo = np.resize(spec.lo, (spec.nzones, 9))
            spec.hi = np.resize(spec.hi, (spec.nzones, 9))
            # Set up temperature brackets.
            for i in range(0, spec.nzones):
                spec.Trange.append(coeffs[int(i*11) + 1])
                # Add the upper bound:
                if i == spec.nzones-1:
                    spec.Trange.append(coeffs[int(i*11) + 2])
            # Read in the coeffs.
            for i in range(0, spec.nzones-1):
                spec.lo[i][:] = coeffs[(int(i*11)+3):(int(i*11)+12)]
                spec.hi[i][:] = coeffs[(int(i*11)+14):(int(i*11)+23)]
            # Fill low range for the last zone as well
            spec.lo[int(spec.nzones-1)][:] = coeffs[((spec.nzones-1)*11+3):((spec.nzones-1)*11+12)]
        else:
            logger = logging.getLogger(__name__)
            logger.error('Unsupported thermo form for species ' + sp)
            sys.exit(1)

        specs.append(spec)

    # Reactions
    reacs = []

    # Cantera internally uses joules/kmol for activation energy
    E_fac = act_energy_fact['joules/kmole']

    def handle_effiencies(reac, ct_rxn):
        """Convert Cantera `cantera.Reaction`'s third body efficienicies
           to pyJac's internal format, and return updated reaction

        Parameters
        ----------
        reac : `ReacInfo`
            The pyJac reaction to update
        ct_rxn : `Reaction` object
            Corresponding cantera reaction to pull the third bodies from

        Returns
        -------
        updated_reac: `ReacInfo`
            The updated pyjac reaction with appropriate third body efficiencies
        """

        # See if single species acts as third body
        if rxn.default_efficiency == 0.0 \
                and len(ct_rxn.efficiencies.keys()) == 1\
                and list(ct_rxn.efficiencies.values())[0] == 1\
                and reac.pdep:
            reac.pdep_sp = list(rxn.efficiencies.keys())[0]
        else:
            for sp in gas.species_names:
                if sp in ct_rxn.efficiencies:
                    reac.thd_body_eff.append([sp, ct_rxn.efficiencies[sp]])
                elif ct_rxn.default_efficiency != 1.0:
                    reac.thd_body_eff.append([sp, ct_rxn.default_efficiency])
        return reac

    for rxn in gas.reactions():

        if isinstance(rxn, ct.ThreeBodyReaction):
            # Instantiate internal reaction based on Cantera Reaction data.
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.rate.pre_exponential_factor,
                                 rxn.rate.temperature_exponent,
                                 rxn.rate.activation_energy * E_fac
                                 )
            reac.thd_body = True
            reac = handle_effiencies(reac, rxn)

        elif isinstance(rxn, ct.FalloffReaction) and \
                not isinstance(rxn, ct.ChemicallyActivatedReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.high_rate.pre_exponential_factor,
                                 rxn.high_rate.temperature_exponent,
                                 rxn.high_rate.activation_energy * E_fac
                                 )
            reac.pdep = True
            reac = handle_effiencies(reac, rxn)

            reac.low = [rxn.low_rate.pre_exponential_factor,
                        rxn.low_rate.temperature_exponent,
                        rxn.low_rate.activation_energy * E_fac
                        ]

            if rxn.falloff.type == 'Troe':
                reac.troe = True
                reac.troe_par = rxn.falloff.parameters.tolist()
                do_warn = False
                if reac.troe_par[1] == 0:
                    reac.troe_par[1] = 1e-30
                    do_warn = True
                if reac.troe_par[2] == 0:
                    reac.troe_par[2] = 1e-30
                    do_warn = True
                if do_warn:
                    logger = logging.getLogger(__name__)
                    logger.warn('Troe parameters in reaction {} modified to avoid'
                                ' division by zero!.'.format(len(reacs)))
            elif rxn.falloff.type == 'SRI':
                reac.sri = True
                reac.sri_par = rxn.falloff.parameters.tolist()

        elif isinstance(rxn, ct.ChemicallyActivatedReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.low_rate.pre_exponential_factor,
                                 rxn.low_rate.temperature_exponent,
                                 rxn.low_rate.activation_energy * E_fac
                                 )
            reac.pdep = True
            reac = handle_effiencies(reac, rxn)

            reac.high = [rxn.high_rate.pre_exponential_factor,
                         rxn.high_rate.temperature_exponent,
                         rxn.high_rate.activation_energy * E_fac
                         ]

            if rxn.falloff.type == 'Troe':
                reac.troe = True
                reac.troe_par = rxn.falloff.parameters.tolist()
                do_warn = False
                if reac.troe_par[1] == 0:
                    reac.troe_par[1] = 1e-30
                    do_warn = True
                if reac.troe_par[2] == 0:
                    reac.troe_par[2] = 1e-30
                    do_warn = True
                if do_warn:
                    logging.warn('Troe parameters in reaction {} modified to avoid'
                                    ' division by zero!.'.format(len(reacs)))
            elif rxn.falloff.type == 'SRI':
                reac.sri = True
                reac.sri_par = rxn.falloff.parameters.tolist()

        elif isinstance(rxn, ct.PlogReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 0.0, 0.0, 0.0
                                 )
            reac.plog = True
            reac.plog_par = []
            for rate in rxn.rates:
                pars = [rate[0], rate[1].pre_exponential_factor,
                        rate[1].temperature_exponent,
                        rate[1].activation_energy * E_fac
                        ]
                reac.plog_par.append(pars)

        elif isinstance(rxn, ct.ChebyshevReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 0.0, 0.0, 0.0
                                 )
            reac.cheb = True
            reac.cheb_n_temp = rxn.nTemperature
            reac.cheb_n_pres = rxn.nPressure
            reac.cheb_plim = [rxn.Pmin, rxn.Pmax]
            reac.cheb_tlim = [rxn.Tmin, rxn.Tmax]
            reac.cheb_par = rxn.coeffs

        elif isinstance(rxn, ct.ElementaryReaction):
            # Instantiate internal reaction based on Cantera Reaction data.

            # Ensure no reactions with zero pre-exponential factor allowed
            if rxn.rate.pre_exponential_factor == 0.0:
                continue

            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.rate.pre_exponential_factor,
                                 rxn.rate.temperature_exponent,
                                 rxn.rate.activation_energy * E_fac
                                 )

        else:
            logger = logging.getLogger(__name__)
            logger.error('Unsupported reaction.')
            sys.exit(1)

        reac.dup = rxn.duplicate

        # No reverse reactions with explicit coefficients in Cantera.

        reacs.append(reac)

    # determine reaction type enums
    for reac in reacs:
        reac.finalize(len(specs))
    for spec in specs:
        spec.finalize()

    if sort_type:
        reacs = sort_reactions(reacs, sort_type)

    return (elems, specs, reacs)

<filename>pyjac/core/mech_interpret.py
# -*- coding: utf-8 -*-
"""Chemkin-format mechanism interpreter module.
"""

# Python 2 compatibility
from __future__ import division

# Standard libraries
import sys
import math
import re
from copy import deepcopy
import logging

import numpy as np

# Local imports
from pyjac import utils
from pyjac.core import chem_model as chem

# Related module
CANTERA_FLAG = False
try:
    import cantera as ct
    version = ct.__version__.split('.')
    logger = logging.getLogger(__name__)
    if int(version[0]) < 2 or int(version[1]) < 3:
        logger.error('Parsing of Cantera mechanisms requires at least version 2.3.0'
                     'in order to access species thermo properties. Detected version'
                     ' is only {}'.format(ct.__version__))
        sys.exit(1)
    CANTERA_FLAG = True
except ImportError:
    CANTERA_FLAG = False

pre_units = ['moles', 'molecules']
"""list(`str`): Supported units list for pre-exponential factor"""

act_energy_units = ['kelvins', 'evolts', 'cal/mole', 'joules/kmole',
                    'kcal/mole', 'joules/mole', 'kjoules/mole'
                    ]
"""list(`str`): Supported units list for activation energy"""

act_energy_fact = dict({'kelvins': 1.0,
                        'evolts': 11595.,
                        'cal/mole': 4.184 / chem.RU_JOUL,
                        'kcal/mole': 4184. / chem.RU_JOUL,
                        'joules/mole': 1. / chem.RU_JOUL,
                        'kjoules/mole': 1000.0 / chem.RU_JOUL,
                        'joules/kmole': 1. / (chem.RU_JOUL * 1000.)
                        })
"""dict: Activation energy conversion factor"""

# get local element atomic weight dict
elem_wt = chem.get_elem_wt()


def sort_reactions(reacs, sort_type, return_order=False):
    """
    reacs : list of `ReacInfo`
        List of reactions in mechanism to be sorted.
    sort_type : :class:`enum_types.reaction_sorting`
        If not None, sort the reactions in this mechanism according to the supplied
        scheme
    return_order : bool [False]
        If True, return new order of sorted reactions, if false, return the sorted
        reactions themselves

    """

    ordering = np.arange(len(reacs))
    from pyjac.core.enum_types import reaction_sorting
    if sort_type is not None and sort_type != reaction_sorting.none:
        from pyjac.core.enum_types import (
            reaction_type, falloff_form, reversible_type, thd_body_type)

        # we consider the following enums:
        # column #1, reaction type
        # column #2, falloff form
        # column #3, third body form
        # column #4, reversible type

        enum_order = (reaction_type, falloff_form, thd_body_type, reversible_type)
        sort_matrix = np.empty((len(reacs), 4))
        for i, rxn in enumerate(reacs):
            for j, enum in enumerate(enum_order):
                e_val = rxn.get_type(enum)
                assert len(e_val) == 1
                sort_matrix[i, j] = int(e_val[0])

        ordering = np.arange(len(reacs))

        # now sort by column -- https://stackoverflow.com/a/38194077
        for i in reversed(range(len(enum_order))):
            inds = sort_matrix[:, i].argsort(kind='mergesort')
            ordering = ordering[inds]
            sort_matrix = sort_matrix[inds]

    if return_order:
        return ordering

    # sort reactions
    return [reacs[i] for i in ordering]


def read_mech(mech_filename, therm_filename, sort_type=None):
    """Read and interpret mechanism file for elements, species, and reactions.

    Parameters
    ----------
    mech_filename : str
        Reaction mechanism filename (e.g. 'mech.dat')
    therm_filename : str, optional
        Thermodynamic database filename (e.g., 'therm.dat')
    sort_type : :class:`enum_types.reaction_sorting`
        If not None, sort the reactions in this mechanism according to the supplied
        scheme

    Returns
    -------
    elems : list of str
        List of elements in mechanism.
    specs : list of `SpecInfo`
        List of species in mechanism.
    reacs : list of `ReacInfo`
        List of reactions in mechanism.
    units : str
        Units of reactions' Arrhenius coefficients

    Notes
    -----
    Doesn't support element names with digits.

    """

    elems = []
    reacs = []
    specs = []
    key = ''

    with open(mech_filename, 'r') as file:
        # start line reading loop
        while True:
            line = file.readline()

            # end of file
            if not line:
                break

            # skip blank or commented lines
            if re.search(r'^\s*$', line) or re.search(r'^\s*!', line):
                continue

            # don't convert to lowercase, since thermo
            # needs to match (for Chemkin)

            # Remove trailing and leading whitespace, tabs, newline
            line = line.strip()

            # remove any comments from end of line
            ind = line.find('!')
            if ind > 0:
                line = line[0:ind]

            # now determine key
            if line[0:4].lower() == 'elem':
                key = 'elem'

                # check for any entries on this line
                line_split = line.split()
                if len(line_split) > 1:
                    ind = line.index(line_split[1])
                    line = line[ind:]
                else:
                    continue

            elif line[0:4].lower() == 'spec':
                key = 'spec'

                # check for any entries on this line
                line_split = line.split()
                if len(line_split) > 1:
                    ind = line.index(line_split[1])
                    line = line[ind:]
                else:
                    continue

            elif line[0:4].lower() == 'reac':
                key = 'reac'

                # default units from Chemkin
                units_E = 'cal/mole'
                units_A = 'moles'

                # get Arrhenius coefficient units (if specified)
                for unit in line.split()[1:]:
                    if unit.lower() in pre_units:
                        units_A = unit.lower()
                    elif unit.lower() in act_energy_units:
                        units_E = unit.lower()
                    else:
                        logger = logging.getLogger(__name__)
                        logger.error('Unsupported units on REACTION line. For '
                                     'pre-exponential factor, choose from: ' +
                                     ', '.join(pre_units) + '\nFor activation energy'
                                     ', choose from: ' + ', '.join(act_energy_units)
                                     + '\nOtherwise leave blank for moles and '
                                     'cal/mole.')
                        sys.exit(1)

                if units_A == 'molecules':
                    raise NotImplementedError('Molecules units not '
                                              'supported, sorry.'
                                              )
                continue
            elif line[0:4].lower() == 'ther':
                # thermo data is in mechanism file
                read_thermo(mech_filename, elems, specs)
                continue
            elif line[0:3].lower() == 'end':
                key = ''
                continue

            if key == 'elem':
                # if any atomic weight declarations, replace / with spaces
                line = line.replace('/', ' ')

                line_split = line.split()
                e_last = ''
                for e in line_split:
                    if e.isalpha():
                        if e[0:3] == 'end':
                            continue
                        if e not in elems:
                            elems.append(e)
                        e_last = e
                    else:
                        # either add new element or update existing
                        # atomic weight
                        elem_wt[e_last.lower()] = float(e)

            elif key == 'spec':
                line_split = line.split()
                for s in line_split:
                    if s[0:3] == 'end':
                        continue
                    if not next((sp for sp in specs if sp.name == s), None):
                        specs.append(chem.SpecInfo(s))

            elif key == 'reac':
                # determine if reaction or auxiliary info line

                if '=' in line:
                    # new reaction

                    cheb_flag = False

                    # get Arrhenius coefficients
                    line_split = line.split()
                    n = len(line_split)
                    reac_A = float(line_split[n - 3])
                    reac_b = float(line_split[n - 2])
                    reac_E = float(line_split[n - 1])

                    ind = line.index(line_split[n - 3])
                    line = line[0:ind].strip()

                    if '<=>' in line:
                        ind = line.index('<=>')
                        reac_rev = True
                        reac_str = line[0:ind].strip()
                        prod_str = line[ind + 3:].strip()
                    elif '=>' in line:
                        ind = line.index('=>')
                        reac_rev = False
                        reac_str = line[0:ind].strip()
                        prod_str = line[ind + 2:].strip()
                    else:
                        ind = line.index('=')
                        reac_rev = True
                        reac_str = line[0:ind].strip()
                        prod_str = line[ind + 1:].strip()

                    thd = False
                    pdep = False
                    pdep_sp = ''

                    reac_spec = []
                    reac_nu = []
                    prod_spec = []
                    prod_nu = []

                    # reactants

                    # look for third-body species
                    sub_str = reac_str
                    while '(' in sub_str:
                        ind1 = sub_str.find('(')
                        ind2 = sub_str.find(')')

                        # Need to check if '+' is first character inside
                        # parentheses and not embedded within parentheses
                        # (e.g., '(+)').
                        # If not, part of species name.
                        inParen = sub_str[ind1 + 1: ind2].strip()
                        if inParen is '+':
                            # '+' embedded within parentheses
                            sub_str = sub_str[ind2 + 1:]
                        elif inParen[0] is '+':
                            pdep = True

                            # either 'm' or a specific species
                            pdep_sp = sub_str[ind1 + 1: ind2].replace('+', ' ')
                            pdep_sp = pdep_sp.strip()

                            if pdep_sp.lower() == 'm':
                                thd = True
                                pdep_sp = ''

                            # now remove from string
                            ind = reac_str.find(sub_str)
                            reac_str = (reac_str[0: ind1 + ind] +
                                        reac_str[ind2 + ind + 1:]
                                        )
                            break
                        else:
                            # Part of species name, remove from substring
                            # and look at rest of reactant line.
                            sub_str = sub_str[ind2 + 1:]

                    reac_list = reac_str.split('+')

                    # Check for empty list elements, meaning there were
                    # multiple '+' in a row, which indicates species'
                    # name ended in '+'.
                    while '' in reac_list:
                        ind = reac_list.index('')
                        reac_list[ind - 1] = reac_list[ind - 1] + '+'
                        del reac_list[ind]

                    # check for any species with '(+)' that was split apart
                    for sp in reac_list:
                        ind = reac_list.index(sp)

                        # ensure not last entry
                        if (ind < len(reac_list) - 1):
                            spNext = reac_list[ind + 1]
                            if sp[len(sp) - 1] is '(' and spNext[0] is ')':
                                reac_list[ind] = sp + '+' + spNext
                                del reac_list[ind + 1]

                    for sp in reac_list:

                        sp = sp.strip()

                        # look for coefficient
                        if sp[0:1].isdigit():
                            # starts with number (coefficient)

                            # search for first letter
                            for i in range(len(sp)):
                                if sp[i: i + 1].isalpha():
                                    break

                            nu = sp[0:i]
                            if '.' in nu:
                                # float
                                nu = float(nu)
                            else:
                                # integer
                                nu = int(nu)

                            sp = sp[i:].strip()
                        else:
                            # no coefficient given
                            nu = 1

                        # check for third body
                        if sp.lower() == 'm':
                            thd = True
                            continue

                        # check if species already in reaction
                        if sp not in reac_spec:
                            # new reactant
                            reac_spec.append(sp)
                            reac_nu.append(nu)
                        else:
                            # existing reactant
                            i = reac_spec.index(sp)
                            reac_nu[i] += nu

                    # products

                    # look for third-body species
                    sub_str = prod_str
                    while '(' in sub_str:
                        ind1 = sub_str.find('(')
                        ind2 = sub_str.find(')')

                        # Need to check if '+' is first character inside
                        # parentheses and not embedded within parentheses
                        # (e.g., '(+)'). If not, part of species name.
                        inParen = sub_str[ind1 + 1: ind2].strip()
                        if inParen is '+':
                            # '+' embedded within parentheses
                            sub_str = sub_str[ind2 + 1:]
                        elif inParen[0] is '+':
                            pdep = True

                            # either 'm' or a specific species
                            pdep_sp = sub_str[ind1 + 1: ind2].replace('+', ' ')
                            pdep_sp = pdep_sp.strip()

                            if pdep_sp.lower() == 'm':
                                thd = True
                                pdep_sp = ''

                            # now remove from string
                            ind = prod_str.find(sub_str)
                            prod_str = (prod_str[0: ind1 + ind] +
                                        prod_str[ind2 + ind + 1:]
                                        )
                            break
                        else:
                            # Part of species name, remove from substring and
                            # look at rest of product line.
                            sub_str = sub_str[ind2 + 1:]

                    prod_list = prod_str.split('+')

                    # Check for empty list elements, meaning there were
                    # multiple '+' in a row, which indicates species
                    # name ended in '+'.
                    while '' in prod_list:
                        ind = prod_list.index('')
                        prod_list[ind - 1] = prod_list[ind - 1] + '+'
                        del prod_list[ind]

                    # check for any species with '(+)' that was split apart
                    for sp in prod_list:
                        ind = prod_list.index(sp)

                        # ensure not last entry
                        if (ind < len(prod_list) - 1):
                            spNext = prod_list[ind + 1]
                            if sp[len(sp) - 1] is '(' and spNext[0] is ')':
                                prod_list[ind] = sp + '+' + spNext
                                del prod_list[ind + 1]

                    for sp in prod_list:

                        sp = sp.strip()

                        # look for coefficient
                        if sp[0:1].isdigit():
                            # starts with number (coefficient)

                            # search for first letter
                            for i in range(len(sp)):
                                if sp[i: i + 1].isalpha():
                                    break

                            nu = sp[0:i]
                            if '.' in nu:
                                # float
                                nu = float(nu)
                            else:
                                # integer
                                nu = int(nu)

                            sp = sp[i:].strip()
                        else:
                            # no coefficient given
                            nu = 1

                        # check for third body
                        if sp in ['m', 'M']:
                            thd = True
                            continue

                        # check if species already in reaction
                        if sp not in prod_spec:
                            # new product
                            prod_spec.append(sp)
                            prod_nu.append(nu)
                        else:
                            # existing product
                            i = prod_spec.index(sp)
                            prod_nu[i] += nu

                    # Don't want to confuse third-body and pressure-dependent
                    # reactions... they are different!
                    if pdep:
                        thd = False

                    # Convert given activation energy units to internal units
                    reac_E *= act_energy_fact[units_E]

                    # Convert given pre-exponential units to internal units
                    if units_A == 'moles':
                        reac_ord = sum(reac_nu)
                        if thd:
                            reac_A /= 1000. ** reac_ord
                        elif pdep:
                            # Low- (chemically activated bimolecular reaction) or
                            # high-pressure (fall-off reaction) limit
                            # parameters
                            reac_A /= 1000. ** (reac_ord - 1.)
                        else:
                            # Elementary reaction
                            reac_A /= 1000. ** (reac_ord - 1.)

                    # add reaction to list
                    reac = chem.ReacInfo(reac_rev, reac_spec, reac_nu,
                                         prod_spec, prod_nu, reac_A, reac_b,
                                         reac_E
                                         )
                    reac.thd_body = thd
                    reac.pdep = pdep
                    if pdep:
                        reac.pdep_sp = pdep_sp

                    reacs.append(reac)

                else:
                    # auxiliary reaction info

                    aux = line[0:3].lower()
                    if aux == 'dup':
                        reacs[-1].dup = True

                    elif aux == 'rev':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        # Convert reverse activation energy units
                        par3 *= act_energy_fact[units_E]

                        # Convert reverse pre-exponential factor
                        if units_A == 'moles':
                            reac_ord = sum(reacs[-1].prod_nu)
                            if thd:
                                par1 /= 1000. ** reac_ord
                            elif pdep:
                                # Low- (chemically activated bimolecular reaction) or
                                # high-pressure (fall-off reaction) limit
                                # parameters
                                par1 /= 1000. ** (reac_ord - 1.)
                            else:
                                # Elementary reaction
                                par1 /= 1000. ** (reac_ord - 1.)

                        # Ensure nonzero reverse coefficients
                        if par1 != 0.0:
                            reacs[-1].rev_par.append(par1)
                            reacs[-1].rev_par.append(par2)
                            reacs[-1].rev_par.append(par3)
                        else:
                            reacs[-1].rev = False

                    elif aux == 'low':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        # Convert low-pressure activation energy units
                        par3 *= act_energy_fact[units_E]

                        # Convert low-pressure pre-exponential factor
                        if units_A == 'moles':
                            par1 /= 1000. ** sum(reacs[-1].reac_nu)

                        reacs[-1].low.append(par1)
                        reacs[-1].low.append(par2)
                        reacs[-1].low.append(par3)

                    elif aux == 'hig':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        # Convert high-pressure activation energy units
                        par3 *= act_energy_fact[units_E]

                        # Convert high-pressure pre-exponential factor
                        if units_A == 'moles':
                            par1 /= 1000. ** (sum(reacs[-1].reac_nu) - 2.)

                        reacs[-1].high.append(par1)
                        reacs[-1].high.append(par2)
                        reacs[-1].high.append(par3)

                    elif aux == 'tro':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        reacs[-1].troe = True
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])

                        do_warn = False
                        if par2 == 0:
                            do_warn = True
                            par2 = 1e-30
                        if par3 == 0:
                            do_warn = True
                            par3 = 1e-30
                        if do_warn:
                            logger = logging.getLogger(__name__)
                            logger.warn(
                                'Troe parameters in reaction {} modified to avoid'
                                ' division by zero!.'.format(len(reacs)))

                        reacs[-1].troe_par.append(par1)
                        reacs[-1].troe_par.append(par2)
                        reacs[-1].troe_par.append(par3)

                        # optional fourth parameter
                        if len(line_split) > 4:
                            par4 = float(line_split[4])
                            reacs[-1].troe_par.append(par4)

                    elif aux == 'sri':
                        line = line.replace('/', ' ')
                        line = line.replace(',', ' ')
                        line_split = line.split()
                        reacs[-1].sri = True
                        par1 = float(line_split[1])
                        par2 = float(line_split[2])
                        par3 = float(line_split[3])
                        reacs[-1].sri_par.append(par1)
                        reacs[-1].sri_par.append(par2)
                        reacs[-1].sri_par.append(par3)

                        # optional fourth and fifth parameters
                        if len(line_split) > 4:
                            par4 = float(line_split[4])
                            par5 = float(line_split[5])
                            reacs[-1].sri_par.append(par4)
                            reacs[-1].sri_par.append(par5)
                    elif aux == 'che':
                        reacs[-1].cheb = True
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        if cheb_flag:
                            for par in line_split[1:]:
                                reacs[-1].cheb_par.append(float(par))
                        else:
                            # first CHEB line
                            cheb_flag = True
                            # Don't want Cheb reactions lumped in with
                            # standard falloff.
                            reacs[-1].pdep = False
                            reacs[-1].cheb_n_temp = int(line_split[1])
                            reacs[-1].cheb_n_pres = int(line_split[2])
                            reacs[-1].cheb_par = []
                            for par in line_split[3:]:
                                reacs[-1].cheb_par.append(float(par))
                    elif aux == 'pch':
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        # Convert pressure from atm to Pa
                        reacs[-1].cheb_plim = [float(line_split[1]) * chem.PA,
                                               float(line_split[2]) * chem.PA
                                               ]

                        # Look for temperature limits on same line:
                        if line_split[3].lower() == 'tcheb':
                            reacs[-1].cheb_tlim = [float(line_split[4]),
                                                   float(line_split[5])
                                                   ]
                    elif aux == 'tch':
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        reacs[-1].cheb_tlim = [float(line_split[1]),
                                               float(line_split[2])
                                               ]
                        # Look for pressure limits on same line:
                        if line_split[3].lower() == 'pcheb':
                            reacs[-1].cheb_plim = [float(line_split[4]) * chem.PA,
                                                   float(
                                                       line_split[5]) * chem.PA
                                                   ]
                    elif aux == 'plo':
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        if not reacs[-1].plog:
                            reacs[-1].plog = True
                            # Don't want Plog reactions lumped in with
                            # standard falloff.
                            reacs[-1].pdep = False
                            reacs[-1].plog_par = []
                        pars = [float(n) for n in line_split[1:5]]

                        # Convert pressure from atm to Pa
                        pars[0] *= 101325.0

                        # Convert given activation energy units to internal
                        # units
                        pars[3] *= act_energy_fact[units_E]

                        # Convert given pre-exponential units to internal units
                        if units_A == 'moles':
                            reac_ord = sum(reacs[-1].reac_nu)
                            # Looks like elementary reaction
                            pars[1] /= 1000. ** (reac_ord - 1.)

                        reacs[-1].plog_par.append(pars)
                    else:
                        # enhanced third body efficiencies
                        line = line.replace('/', ' ')
                        line_split = line.split()
                        for i in range(0, len(line_split), 2):
                            pair = [line_split[i], float(line_split[i + 1])]
                            reacs[-1].thd_body_eff.append(pair)

    # process some reaction auxiliary info
    for idx, reac in enumerate(reacs):
        if reac.cheb:
            # check for correct number
            n = reac.cheb_n_temp
            m = reac.cheb_n_pres
            if len(reac.cheb_par) != n * m:
                logger = logging.getLogger(__name__)
                logging.error('Incorrect number of CHEB coefficients in reaction '
                              + repr(idx))
                sys.exit(1)
            else:
                # Convert units of first Chebyshev parameter
                order = sum(reac.reac_nu)
                if units_A == 'moles':
                    reac.cheb_par[0] += math.log10(0.001 ** (order - 1.))

                reacs[idx].cheb_par = np.reshape(reac.cheb_par, (n, m))

    # check that all species in reactions correspond to a known species
    spec_names = set(spec.name for spec in specs)
    for idx, reac in enumerate(reacs):
        in_rxn = set(reac.reac + reac.prod)
        for spec in in_rxn:
            if spec not in spec_names:
                logger = logging.getLogger(__name__)
                logger.error('Reaction {} contains unknown species {}'.format(
                    idx, spec))
                sys.exit(-1)

    # Split reversible reactions with explicit reverse parameters into
    # two irreversible reactions to match Cantera's behavior
    for reac in reacs[:]:
        if reac.rev_par:
            new_reac = deepcopy(reac)

            idx = reacs.index(reac)
            reacs[idx].rev = False
            reacs[idx].rev_par = []

            new_reac.A = new_reac.rev_par[0]
            new_reac.b = new_reac.rev_par[1]
            new_reac.E = new_reac.rev_par[2]
            new_reac.rev = False
            new_reac.rev_par = []
            new_reac.prod = reac.reac[:]
            new_reac.prod_nu = reac.reac_nu[:]
            new_reac.reac = reac.prod[:]
            new_reac.reac_nu = reac.prod_nu[:]

            reacs.insert(idx + 1, new_reac)

    # Read seperate thermo file if present and needed
    if any([not spec.mw for spec in specs]):
        if therm_filename:
            read_thermo(therm_filename, elems, specs)
        else:
            logger = logging.getLogger(__name__)
            logger.error(
                'No thermo file specified, but species missing \n'
                'data. Either specify file, or ensure complete data in\n'
                'mechanism file with THERMO option.')
            sys.exit(1)

    # Check for missing thermo data again
    missing_mw = [spec.name for spec in specs if not spec.mw]
    if missing_mw:
        logger = logging.getLogger(__name__)
        logger.error('Missing thermo data for ' + ', '.join(missing_mw))
        sys.exit(1)

    # determine reaction type enums
    for reac in reacs:
        reac.finalize(len(specs))
    for spec in specs:
        spec.finalize()

    if sort_type:
        reacs = sort_reactions(reacs, sort_type)

    return (elems, specs, reacs)


def read_thermo(filename, elems, specs):
    """Read and interpret thermodynamic database for species data.

    Reads the thermodynamic file and returns the species thermodynamic
    coefficients as well as the species-specific temperature range
    values (if given).

    Parameters
    ----------
    filename : str
        Name of thermo database file.
    elems : list of str
        List of element names in mechanism.
    specs : list of `SpecInfo`
        List of species in mechanism.

    Returns
    -------
    None

    """

    # choose our element names to enable comparison between Chemkin & Cantera
    elem_name_map = {e.lower(): e for e in elem_wt}

    with open(filename, 'r') as file:

        # loop through intro lines
        while True:
            line = file.readline()

            # skip blank or commented lines
            if re.search(r'^\s*$', line) or re.search(r'^\s*!', line):
                continue

            # skip 'thermo' at beginning
            if 'thermo' in line.lower():
                break

        # next line either has common temperature ranges or first species
        last_line = file.tell()
        line = file.readline()

        line_split = line.split()
        if line_split[0][0:1].isdigit():
            T_ranges = utils.read_str_num(line)
        else:
            # no common temperature info
            file.seek(last_line)
            # default

        # now start reading species thermo info
        while True:
            # first line of species info
            line = file.readline()

            # don't convert to lowercase, needs to match thermo for Chemkin

            # break if end of file
            if line is None or line[0:3].lower() == 'end':
                break

            # skip blank/commented line
            if re.search(r'^\s*$', line) or re.search(r'^\s*!', line):
                continue

            # species name, columns 0:18
            spec = line[0:18].strip()

            # Apparently, in some cases, notes are in the
            # columns of shorter species names, so make
            # sure no spaces.
            if spec.find(' ') > 0:
                spec = spec[0: spec.find(' ')]

            # now need to determine if this species is in mechanism
            if next((sp for sp in specs if sp.name == spec), None):
                sp_ind = next(i for i in range(len(specs))
                              if specs[i].name == spec
                              )
            else:
                # not in mechanism, read next three lines and continue
                line = file.readline()
                line = file.readline()
                line = file.readline()
                continue

            # set species to the one matched
            spec = specs[sp_ind]

            # ensure not reading the same species more than once...
            if spec.mw:
                # already done! skip next three lines
                line = file.readline()
                line = file.readline()
                line = file.readline()
                continue

            # now get element composition of species, columns 24:44
            # each piece of data is 5 characters long (2 for element, 3 for #)
            elem_str = utils.split_str(line[24:44], 5)

            for e_str in elem_str:
                e = e_str[0:2].strip()
                # skip if blank
                if e == '' or e == '0':
                    continue
                # may need to convert to float first, in case of e.g. "1."
                e_num = float(e_str[2:].strip())
                e_num = int(e_num)

                spec.elem.append([elem_name_map[e.lower()], e_num])

                # calculate molecular weight
                spec.mw += e_num * elem_wt[e.lower()]

            # temperatures for species
            T_spec = utils.read_str_num(line[45:74])
            T_low = T_spec[0]
            T_high = T_spec[1]
            if len(T_spec) == 3:
                T_com = T_spec[2]
            else:
                T_com = T_ranges[1]

            spec.Trange = [T_low, T_com, T_high]

            # second species line
            line = file.readline()
            coeffs = utils.split_str(line[0:75], 15)
            spec.hi[0] = float(coeffs[0])
            spec.hi[1] = float(coeffs[1])
            spec.hi[2] = float(coeffs[2])
            spec.hi[3] = float(coeffs[3])
            spec.hi[4] = float(coeffs[4])

            # third species line
            line = file.readline()
            coeffs = utils.split_str(line[0:75], 15)
            spec.hi[5] = float(coeffs[0])
            spec.hi[6] = float(coeffs[1])
            spec.lo[0] = float(coeffs[2])
            spec.lo[1] = float(coeffs[3])
            spec.lo[2] = float(coeffs[4])

            # fourth species line
            line = file.readline()
            coeffs = utils.split_str(line[0:75], 15)
            spec.lo[3] = float(coeffs[0])
            spec.lo[4] = float(coeffs[1])
            spec.lo[5] = float(coeffs[2])
            spec.lo[6] = float(coeffs[3])

            # stop reading if all species in mechanism accounted for
            if not next((sp for sp in specs if sp.mw == 0.0), None):
                break

    return None


def read_mech_ct(filename=None, gas=None, sort_type=None):
    """Read and interpret Cantera-format mechanism file.

    Parameters
    ----------
    filename : str
        Reaction mechanism filename (e.g. 'mech.cti'). Optional.
    gas : `cantera.Solution` object
        Existing Cantera Solution object to be used. Optional.
    sort_type : :class:`enum_types.reaction_sorting`
        If not None, sort the reactions in this mechanism according to the supplied
        scheme

    Returns
    -------
    elems : list of str
        List of elements in mechanism.
    specs : list of `SpecInfo`
        List of species in mechanism.
    reacs : list of `ReacInfo`
        List of reactions in mechanism.
    units : str
        Units of reactions' Arrhenius coefficients

    """

    if not CANTERA_FLAG:
        logger = logging.getLogger(__name__)
        logger.error('Cantera not installed. Cannot interpret Cantera-format'
                     ' mechanism.')
        sys.exit(1)

    if filename:
        gas = ct.Solution(filename)
    elif not gas:
        logger = logging.getLogger(__name__)
        logger.error('Need either filename or Cantera Solution object.')
        sys.exit(1)

    # Elements
    elems = gas.element_names
    for e, wt in zip(elems, gas.atomic_weights):
        if e.lower() not in elem_wt:
            elem_wt[e.lower()] = wt

    # choose our element names to enable comparison between Chemkin & Cantera
    elem_name_map = {e.lower(): e for e in elem_wt}

    # Species
    specs = []
    for i, sp in enumerate(gas.species_names):
        # Get Species object
        species = gas.species(i)

        if isinstance(species.thermo, ct.NasaPoly2):
            spec = chem.SpecInfo(sp)
        elif isinstance(species.thermo, ct.Nasa9PolyMultiTempRegion):
            spec = chem.SpecInfoN9(sp)

        spec.mw = gas.molecular_weights[i]

        # Species elemental composition
        for e in species.composition:
            spec.elem.append([elem_name_map[e.lower()], int(species.composition[e])])

        # Species thermodynamic properties
        coeffs = species.thermo.coeffs
        if isinstance(species.thermo, ct.NasaPoly2):
            spec.Trange = [species.thermo.min_temp, coeffs[0],
                            species.thermo.max_temp
                            ]
            spec.hi = coeffs[1:8]
            spec.lo = coeffs[8:15]
        elif isinstance(species.thermo, ct.Nasa9PolyMultiTempRegion):
            spec.nzones = int(coeffs[0])
            spec.lo = np.resize(spec.lo, (spec.nzones, 9))
            spec.hi = np.resize(spec.hi, (spec.nzones, 9))
            # Set up temperature brackets.
            for i in range(0, spec.nzones):
                spec.Trange.append(coeffs[int(i*11) + 1])
                # Add the upper bound:
                if i == spec.nzones-1:
                    spec.Trange.append(coeffs[int(i*11) + 2])
            # Read in the coeffs.
            for i in range(0, spec.nzones-1):
                spec.lo[i][:] = coeffs[(int(i*11)+3):(int(i*11)+12)]
                spec.hi[i][:] = coeffs[(int(i*11)+14):(int(i*11)+23)]
            # Fill low range for the last zone as well
            spec.lo[int(spec.nzones-1)][:] = coeffs[((spec.nzones-1)*11+3):((spec.nzones-1)*11+12)]
        else:
            logger = logging.getLogger(__name__)
            logger.error('Unsupported thermo form for species ' + sp)
            sys.exit(1)

        specs.append(spec)

    # Reactions
    reacs = []

    # Cantera internally uses joules/kmol for activation energy
    E_fac = act_energy_fact['joules/kmole']

    def handle_effiencies(reac, ct_rxn):
        """Convert Cantera `cantera.Reaction`'s third body efficienicies
           to pyJac's internal format, and return updated reaction

        Parameters
        ----------
        reac : `ReacInfo`
            The pyJac reaction to update
        ct_rxn : `Reaction` object
            Corresponding cantera reaction to pull the third bodies from

        Returns
        -------
        updated_reac: `ReacInfo`
            The updated pyjac reaction with appropriate third body efficiencies
        """

        # See if single species acts as third body
        if rxn.default_efficiency == 0.0 \
                and len(ct_rxn.efficiencies.keys()) == 1\
                and list(ct_rxn.efficiencies.values())[0] == 1\
                and reac.pdep:
            reac.pdep_sp = list(rxn.efficiencies.keys())[0]
        else:
            for sp in gas.species_names:
                if sp in ct_rxn.efficiencies:
                    reac.thd_body_eff.append([sp, ct_rxn.efficiencies[sp]])
                elif ct_rxn.default_efficiency != 1.0:
                    reac.thd_body_eff.append([sp, ct_rxn.default_efficiency])
        return reac

    for rxn in gas.reactions():

        if isinstance(rxn, ct.ThreeBodyReaction):
            # Instantiate internal reaction based on Cantera Reaction data.
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.rate.pre_exponential_factor,
                                 rxn.rate.temperature_exponent,
                                 rxn.rate.activation_energy * E_fac
                                 )
            reac.thd_body = True
            reac = handle_effiencies(reac, rxn)

        elif isinstance(rxn, ct.FalloffReaction) and \
                not isinstance(rxn, ct.ChemicallyActivatedReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.high_rate.pre_exponential_factor,
                                 rxn.high_rate.temperature_exponent,
                                 rxn.high_rate.activation_energy * E_fac
                                 )
            reac.pdep = True
            reac = handle_effiencies(reac, rxn)

            reac.low = [rxn.low_rate.pre_exponential_factor,
                        rxn.low_rate.temperature_exponent,
                        rxn.low_rate.activation_energy * E_fac
                        ]

            if rxn.falloff.type == 'Troe':
                reac.troe = True
                reac.troe_par = rxn.falloff.parameters.tolist()
                do_warn = False
                if reac.troe_par[1] == 0:
                    reac.troe_par[1] = 1e-30
                    do_warn = True
                if reac.troe_par[2] == 0:
                    reac.troe_par[2] = 1e-30
                    do_warn = True
                if do_warn:
                    logger = logging.getLogger(__name__)
                    logger.warn('Troe parameters in reaction {} modified to avoid'
                                ' division by zero!.'.format(len(reacs)))
            elif rxn.falloff.type == 'SRI':
                reac.sri = True
                reac.sri_par = rxn.falloff.parameters.tolist()

        elif isinstance(rxn, ct.ChemicallyActivatedReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.low_rate.pre_exponential_factor,
                                 rxn.low_rate.temperature_exponent,
                                 rxn.low_rate.activation_energy * E_fac
                                 )
            reac.pdep = True
            reac = handle_effiencies(reac, rxn)

            reac.high = [rxn.high_rate.pre_exponential_factor,
                         rxn.high_rate.temperature_exponent,
                         rxn.high_rate.activation_energy * E_fac
                         ]

            if rxn.falloff.type == 'Troe':
                reac.troe = True
                reac.troe_par = rxn.falloff.parameters.tolist()
                do_warn = False
                if reac.troe_par[1] == 0:
                    reac.troe_par[1] = 1e-30
                    do_warn = True
                if reac.troe_par[2] == 0:
                    reac.troe_par[2] = 1e-30
                    do_warn = True
                if do_warn:
                    logging.warn('Troe parameters in reaction {} modified to avoid'
                                    ' division by zero!.'.format(len(reacs)))
            elif rxn.falloff.type == 'SRI':
                reac.sri = True
                reac.sri_par = rxn.falloff.parameters.tolist()

        elif isinstance(rxn, ct.PlogReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 0.0, 0.0, 0.0
                                 )
            reac.plog = True
            reac.plog_par = []
            for rate in rxn.rates:
                pars = [rate[0], rate[1].pre_exponential_factor,
                        rate[1].temperature_exponent,
                        rate[1].activation_energy * E_fac
                        ]
                reac.plog_par.append(pars)

        elif isinstance(rxn, ct.ChebyshevReaction):
            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 0.0, 0.0, 0.0
                                 )
            reac.cheb = True
            reac.cheb_n_temp = rxn.nTemperature
            reac.cheb_n_pres = rxn.nPressure
            reac.cheb_plim = [rxn.Pmin, rxn.Pmax]
            reac.cheb_tlim = [rxn.Tmin, rxn.Tmax]
            reac.cheb_par = rxn.coeffs

        elif isinstance(rxn, ct.ElementaryReaction):
            # Instantiate internal reaction based on Cantera Reaction data.

            # Ensure no reactions with zero pre-exponential factor allowed
            if rxn.rate.pre_exponential_factor == 0.0:
                continue

            reac = chem.ReacInfo(rxn.reversible,
                                 list(rxn.reactants.keys()),
                                 list(rxn.reactants.values()),
                                 list(rxn.products.keys()),
                                 list(rxn.products.values()),
                                 rxn.rate.pre_exponential_factor,
                                 rxn.rate.temperature_exponent,
                                 rxn.rate.activation_energy * E_fac
                                 )

        else:
            logger = logging.getLogger(__name__)
            logger.error('Unsupported reaction.')
            sys.exit(1)

        reac.dup = rxn.duplicate

        # No reverse reactions with explicit coefficients in Cantera.

        reacs.append(reac)

    # determine reaction type enums
    for reac in reacs:
        reac.finalize(len(specs))
    for spec in specs:
        spec.finalize()

    if sort_type:
        reacs = sort_reactions(reacs, sort_type)

    return (elems, specs, reacs)

# -*- coding: utf-8 -*- Chemkin-format mechanism interpreter module. # Python 2 compatibility # Standard libraries # Local imports # Related module list(`str`): Supported units list for pre-exponential factor list(`str`): Supported units list for activation energy dict: Activation energy conversion factor # get local element atomic weight dict reacs : list of `ReacInfo`         List of reactions in mechanism to be sorted.     sort_type : :class:`enum_types.reaction_sorting`         If not None, sort the reactions in this mechanism according to the supplied         scheme     return_order : bool [False]         If True, return new order of sorted reactions, if false, return the sorted         reactions themselves # we consider the following enums: # column #1, reaction type # column #2, falloff form # column #3, third body form # column #4, reversible type # now sort by column -- https://stackoverflow.com/a/38194077 # sort reactions Read and interpret mechanism file for elements, species, and reactions.      Parameters     ----------     mech_filename : str         Reaction mechanism filename (e.g. 'mech.dat')     therm_filename : str, optional         Thermodynamic database filename (e.g., 'therm.dat')     sort_type : :class:`enum_types.reaction_sorting`         If not None, sort the reactions in this mechanism according to the supplied         scheme      Returns     -------     elems : list of str         List of elements in mechanism.     specs : list of `SpecInfo`         List of species in mechanism.     reacs : list of `ReacInfo`         List of reactions in mechanism.     units : str         Units of reactions' Arrhenius coefficients      Notes     -----     Doesn't support element names with digits. # start line reading loop # end of file # skip blank or commented lines # don't convert to lowercase, since thermo # needs to match (for Chemkin) # Remove trailing and leading whitespace, tabs, newline # remove any comments from end of line # now determine key # check for any entries on this line # check for any entries on this line # default units from Chemkin # get Arrhenius coefficient units (if specified) # thermo data is in mechanism file # if any atomic weight declarations, replace / with spaces # either add new element or update existing # atomic weight # determine if reaction or auxiliary info line # new reaction # get Arrhenius coefficients # reactants # look for third-body species # Need to check if '+' is first character inside # parentheses and not embedded within parentheses # (e.g., '(+)'). # If not, part of species name. # '+' embedded within parentheses # either 'm' or a specific species # now remove from string # Part of species name, remove from substring # and look at rest of reactant line. # Check for empty list elements, meaning there were # multiple '+' in a row, which indicates species' # name ended in '+'. # check for any species with '(+)' that was split apart # ensure not last entry # look for coefficient # starts with number (coefficient) # search for first letter # float # integer # no coefficient given # check for third body # check if species already in reaction # new reactant # existing reactant # products # look for third-body species # Need to check if '+' is first character inside # parentheses and not embedded within parentheses # (e.g., '(+)'). If not, part of species name. # '+' embedded within parentheses # either 'm' or a specific species # now remove from string # Part of species name, remove from substring and # look at rest of product line. # Check for empty list elements, meaning there were # multiple '+' in a row, which indicates species # name ended in '+'. # check for any species with '(+)' that was split apart # ensure not last entry # look for coefficient # starts with number (coefficient) # search for first letter # float # integer # no coefficient given # check for third body # check if species already in reaction # new product # existing product # Don't want to confuse third-body and pressure-dependent # reactions... they are different! # Convert given activation energy units to internal units # Convert given pre-exponential units to internal units # Low- (chemically activated bimolecular reaction) or # high-pressure (fall-off reaction) limit # parameters # Elementary reaction # add reaction to list # auxiliary reaction info # Convert reverse activation energy units # Convert reverse pre-exponential factor # Low- (chemically activated bimolecular reaction) or # high-pressure (fall-off reaction) limit # parameters # Elementary reaction # Ensure nonzero reverse coefficients # Convert low-pressure activation energy units # Convert low-pressure pre-exponential factor # Convert high-pressure activation energy units # Convert high-pressure pre-exponential factor # optional fourth parameter # optional fourth and fifth parameters # first CHEB line # Don't want Cheb reactions lumped in with # standard falloff. # Convert pressure from atm to Pa # Look for temperature limits on same line: # Look for pressure limits on same line: # Don't want Plog reactions lumped in with # standard falloff. # Convert pressure from atm to Pa # Convert given activation energy units to internal # units # Convert given pre-exponential units to internal units # Looks like elementary reaction # enhanced third body efficiencies # process some reaction auxiliary info # check for correct number # Convert units of first Chebyshev parameter # check that all species in reactions correspond to a known species # Split reversible reactions with explicit reverse parameters into # two irreversible reactions to match Cantera's behavior # Read seperate thermo file if present and needed # Check for missing thermo data again # determine reaction type enums Read and interpret thermodynamic database for species data.      Reads the thermodynamic file and returns the species thermodynamic     coefficients as well as the species-specific temperature range     values (if given).      Parameters     ----------     filename : str         Name of thermo database file.     elems : list of str         List of element names in mechanism.     specs : list of `SpecInfo`         List of species in mechanism.      Returns     -------     None # choose our element names to enable comparison between Chemkin & Cantera # loop through intro lines # skip blank or commented lines # skip 'thermo' at beginning # next line either has common temperature ranges or first species # no common temperature info # default # now start reading species thermo info # first line of species info # don't convert to lowercase, needs to match thermo for Chemkin # break if end of file # skip blank/commented line # species name, columns 0:18 # Apparently, in some cases, notes are in the # columns of shorter species names, so make # sure no spaces. # now need to determine if this species is in mechanism # not in mechanism, read next three lines and continue # set species to the one matched # ensure not reading the same species more than once... # already done! skip next three lines # now get element composition of species, columns 24:44 # each piece of data is 5 characters long (2 for element, 3 for #) # skip if blank # may need to convert to float first, in case of e.g. "1." # calculate molecular weight # temperatures for species # second species line # third species line # fourth species line # stop reading if all species in mechanism accounted for Read and interpret Cantera-format mechanism file.      Parameters     ----------     filename : str         Reaction mechanism filename (e.g. 'mech.cti'). Optional.     gas : `cantera.Solution` object         Existing Cantera Solution object to be used. Optional.     sort_type : :class:`enum_types.reaction_sorting`         If not None, sort the reactions in this mechanism according to the supplied         scheme      Returns     -------     elems : list of str         List of elements in mechanism.     specs : list of `SpecInfo`         List of species in mechanism.     reacs : list of `ReacInfo`         List of reactions in mechanism.     units : str         Units of reactions' Arrhenius coefficients # Elements # choose our element names to enable comparison between Chemkin & Cantera # Species # Get Species object # Species elemental composition # Species thermodynamic properties # Set up temperature brackets. # Add the upper bound: # Read in the coeffs. # Fill low range for the last zone as well # Reactions # Cantera internally uses joules/kmol for activation energy Convert Cantera `cantera.Reaction`'s third body efficienicies            to pyJac's internal format, and return updated reaction          Parameters         ----------         reac : `ReacInfo`             The pyJac reaction to update         ct_rxn : `Reaction` object             Corresponding cantera reaction to pull the third bodies from          Returns         -------         updated_reac: `ReacInfo`             The updated pyjac reaction with appropriate third body efficiencies # See if single species acts as third body # Instantiate internal reaction based on Cantera Reaction data. # Instantiate internal reaction based on Cantera Reaction data. # Ensure no reactions with zero pre-exponential factor allowed # No reverse reactions with explicit coefficients in Cantera. # determine reaction type enums