/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::pressureDependentArrheniusReactionRate

Description
    Arrhenius pressure dependent reaction rate analog to PLOG in Chemkin.
    Where the parameters are interpolated lineraly between the logarithmic
    pressures (see e.g. cantera)

    The parameters in the reaction file are stored as following:

    un-named-reaction
    {
        type       ...PressureDependentArrhenius
        reaction   "A + B = C + D";
        A           table (
                    (4   -17)
                    (10  7)
                    (18  75)
        );
        beta        table (
                    (4    2)
                    (10  -5)
                    (18  -6)
        );
        Ta          table (
                    (4   4400)
                    (10  12000)
                    (18  14000)
        );
    }

    the pressure is stored as log(p)
    also the values for A are stored as log(A)
    to make use of the available interpolation Table in OpenFOAM

SourceFiles
    pressureDependentArrheniusReactionRateI.H

\*---------------------------------------------------------------------------*/

#ifndef pressureDependentArrheniusReactionRate_H
#define pressureDependentArrheniusReactionRate_H

#include "scalarField.H"
#include "typeInfo.H"
#include "Table.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{

// Forward declaration of friend functions and operators

class pressureDependentArrheniusReactionRate;

Ostream& operator<<(Ostream&, const pressureDependentArrheniusReactionRate&);


/*---------------------------------------------------------------------------*\
               Class pressureDependentArrheniusReactionRate Declaration
\*---------------------------------------------------------------------------*/

class pressureDependentArrheniusReactionRate
{

    // Private Data
    Function1s::Table<scalar> A_;
    Function1s::Table<scalar> beta_;
    Function1s::Table<scalar> Ta_;


public:

    // Constructors

        //- Construct from components
        inline pressureDependentArrheniusReactionRate
        (
            Function1s::Table<scalar> A,
            Function1s::Table<scalar> beta,
            Function1s::Table<scalar> Ta
        );

        //- Construct from dictionary
        inline pressureDependentArrheniusReactionRate
        (
            const speciesTable& species,
            const dictionary& dict
        );


    // Member Functions

        //- Return the type name
        static word type()
        {
            return "pressureDependentArrhenius";
        }

        //- Pre-evaluation hook
        inline void preEvaluate() const;

        //- Post-evaluation hook
        inline void postEvaluate() const;

        //- Return the rate
        inline scalar operator()
        (
            const scalar p,
            const scalar T,
            const scalarField& c,
            const label li
        ) const;

        //- The derivative of the rate w.r.t. temperature
        inline scalar ddT
        (
            const scalar p,
            const scalar T,
            const scalarField& c,
            const label li
        ) const;

        //- Third-body efficiencies (beta = 1-alpha)
        //  non-empty only for third-body reactions
        //  with enhanced molecularity (alpha != 1)
        inline const List<Tuple2<label, scalar>>& beta() const;

        //- Species concentration derivative of the pressure dependent term
        //  By default this value is 1 as it multiplies the third-body term
        inline void dcidc
        (
            const scalar p,
            const scalar T,
            const scalarField& c,
            const label li,
            scalarField& dcidc
        ) const;

        //- Temperature derivative of the pressure dependent term
        //  By default this value is 0 since ddT of molecularity is approx.0
        inline scalar dcidT
        (
            const scalar p,
            const scalar T,
            const scalarField& c,
            const label li
        ) const;


        //- Write to stream
        inline void write(Ostream& os) const;


    // Ostream Operator

        inline friend Ostream& operator<<
        (
            Ostream&,
            const pressureDependentArrheniusReactionRate&
        );
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#include "pressureDependentArrheniusReactionRateI.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

// * * * * * * * * * * * * * * * * Constructors  * * * * * * * * * * * * * * //

inline Foam::pressureDependentArrheniusReactionRate::pressureDependentArrheniusReactionRate
(
    Function1s::Table<scalar> A,
    Function1s::Table<scalar> beta,
    Function1s::Table<scalar> Ta
)
:
    A_(A),
    beta_(beta),
    Ta_(Ta)
{}


inline Foam::pressureDependentArrheniusReactionRate::pressureDependentArrheniusReactionRate
(
    const speciesTable& species,
    const dictionary& dict
)
:
    A_("A", dict),
    beta_("beta", dict),
    Ta_("Ta", dict)
{}


// * * * * * * * * * * * * * * * Member Functions  * * * * * * * * * * * * * //

inline void Foam::pressureDependentArrheniusReactionRate::preEvaluate() const
{}


inline void Foam::pressureDependentArrheniusReactionRate::postEvaluate() const
{}


inline Foam::scalar Foam::pressureDependentArrheniusReactionRate::operator()
(
    const scalar p,
    const scalar T,
    const scalarField&,
    const label
) const
{
    const scalar logP = log(p);
    scalar ak = exp(A_.value(logP));

    if (mag(beta_.value(logP)) > vSmall)
    {
        ak *= pow(T, beta_.value(logP));
    }

    if (mag(Ta_.value(logP)) > vSmall)
    {
        ak *= exp(-Ta_.value(logP)/T);
    }

    return ak;
}


inline Foam::scalar Foam::pressureDependentArrheniusReactionRate::ddT
(
    const scalar p,
    const scalar T,
    const scalarField&,
    const label
) const
{
    const scalar logP = log(p);
    scalar ak = exp(A_.value(logP));

    if (mag(beta_.value(logP)) > vSmall)
    {
        ak *= pow(T, beta_.value(logP));
    }

    if (mag(Ta_.value(logP)) > vSmall)
    {
        ak *= exp(-Ta_.value(logP)/T);
    }

    return ak*(beta_.value(logP)+Ta_.value(logP)/T)/T;
}


inline const Foam::List<Foam::Tuple2<Foam::label, Foam::scalar>>&
Foam::pressureDependentArrheniusReactionRate::beta() const
{
    return NullObjectRef<List<Tuple2<label, scalar>>>();
}



inline void Foam::pressureDependentArrheniusReactionRate::dcidc
(
    const scalar p,
    const scalar T,
    const scalarField& c,
    const label li,
    scalarField& dcidc
) const
{}

inline Foam::scalar Foam::pressureDependentArrheniusReactionRate::dcidT
(
    const scalar p,
    const scalar T,
    const scalarField& c,
    const label li
) const
{
    return 0;
}

inline void Foam::pressureDependentArrheniusReactionRate::write(Ostream& os) const
{
    //convert to scalarFields to be able to loop over values
    scalarField pressureTable(A_.x());
    scalarField ATable(A_.y());
    scalarField betaTable(beta_.y());
    scalarField TATable(Ta_.y());
    os  << indent << "A               table ("  << nl ; 
    os << incrIndent;
    forAll(pressureTable, i)
    {
        os << indent << "                  (" 
           << pressureTable[i] << "  " << ATable[i] << " )" << nl;
    }
    os << indent << ");" << nl
       << decrIndent
       << indent << "beta            table (" << nl; 
    os << incrIndent;
    forAll(pressureTable, i)
    {
        os << indent << "                  (" 
           << pressureTable[i] << "  " << betaTable[i] << " )" << nl;
    }
    os << indent << ");" << nl 
       << decrIndent 
       << indent << "Ta            table (" << nl;
    os << incrIndent;
    forAll(pressureTable, i)
    {
        os << indent
            << "                  (" 
            << pressureTable[i] << "  " << TATable[i] << " )" << nl;
    }
    os << indent << ");" << nl
       << decrIndent;
}


inline Foam::Ostream& Foam::operator<<
(
    Ostream& os,
    const pressureDependentArrheniusReactionRate& arr
)
{
    arr.write(os);
    return os;
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2021 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "makeReaction.H"

#include "pressureDependentArrheniusReactionRate.H"

#include "forGases.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{
    forCoeffGases(makeIRReactions, pressureDependentArrheniusReactionRate);
}

// ************************************************************************* //

/*--------------------------------*- C++ -*----------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

%{

#undef yyFlexLexer

#include "chemkinIIReader.H"

#include "error.H"
#include "IStringStream.H"
// For EOF only
#include <cstdio>

// flex input buffer size
int Foam::chemkinIIReader::yyBufSize = YY_BUF_SIZE;

// Dummy yyFlexLexer::yylex() to keep the linker happy. It is not called
//! \cond dummy
int yyFlexLexer::yylex()
{
    FatalErrorInFunction
        << "should not have called this function"
        << abort(Foam::FatalError);

    return 0;
}
//! \endcond


// Dummy yywrap to keep yylex happy at compile time.
// It is called by yylex but is not used as the mechanism to change file.
// See <<EOF>>
//! \cond dummy
#if YY_FLEX_MINOR_VERSION < 6 && YY_FLEX_SUBMINOR_VERSION < 34
extern "C" int yywrap()
#else
int yyFlexLexer::yywrap()
#endif
{
    return 1;
}
//! \endcond


Foam::string foamSpecieString(const char* YYText)
{
    Foam::string specieString(YYText);
    return specieString;
}

class foamSpecieName
{
public:

    //- Is this character valid for a foamSpecieName
    inline static bool valid(char c)
    {
        return
        (
            !isspace(c)
            && c != '"'   // string quote
            && c != '\''  // string quote
            && c != '/'   // slash
            && c != ';'   // end statement
            && c != '{'   // beg subdict
            && c != '}'   // end subdict
        );
    }
};


Foam::word foamName(const char* YYText)
{
    Foam::string fn(YYText);
    Foam::string::stripInvalid<foamSpecieName>(fn);
    return fn;
}


Foam::word foamName(const Foam::string& s)
{
    Foam::string fn(s);
    Foam::string::stripInvalid<foamSpecieName>(fn);
    return fn;
}


 /* ------------------------------------------------------------------------- *\
   ------ cppLexer::yylex()
 \* ------------------------------------------------------------------------- */

#define YY_DECL int Foam::chemkinIIReader::lex()

%}

one_space             [ \t\f\r]
space                 {one_space}*
some_space            {one_space}+
newline               [ \n]

alpha                 [_A-Za-z]
digit                 [0-9]

exponent_part         [eEdD][-+]?{digit}+
fractional_constant   [-+]?(({digit}*"."{digit}+)|({digit}+"."?))

floatNum              (({fractional_constant}{exponent_part}?)|({digit}+{exponent_part}))


elements              {space}("ELEMENTS"|"ELEM"){space}
species               {space}("SPECIES"|"SPECIE"|"SPEC"){space}
thermoAll             {space}"THERMO"{some_space}"ALL"{space}
thermo                {space}"THERMO"{space}
reactions             {space}("REACTIONS"|"REAC"){space}
end                   {space}"END"{space}

elementName           {space}([A-Z]|([A-Z][A-Z])){space}
startIsotopeMolW      {space}"/"{space}
isotopeMolW           {space}{floatNum}{space}"/"{space}

specieName            {space}[A-Za-z](([A-Za-z0-9)*+-])|("("[^+]))*{space}
Y                     {space}{floatNum}{space}

thermoTemp            .{10}

thermoSpecieName      .{18}
thermoDate            .{6}
thermoFormula         .{20}
thermoPhase           S|L|G|C
thermoLowTemp         .{10}
thermoHighTemp        .{10}
thermoCommonTemp      .{8}
thermoFormula2        .{5}
thermoCoeff           .{15}
thermoLineLabel1      " "1{space}
thermoLineLabel2      " "{4}2{space}
thermoLineLabel3      " "{4}3{space}
thermoLineLabel4      " "{4}4{space}

specieDelimiter       {space}"+"{space}
reversibleReactionDelimiter {space}("="|"<=>"){space}
irreversibleReactionDelimiter {space}"=>"{space}
startPDependentSpecie {space}"("{space}"+"{space}
pDependentSpecie      {specieName}")"{space}
reactionCoeffs        {space}{floatNum}{some_space}{floatNum}{some_space}{floatNum}{space}
reactionKeyword       {space}[A-Za-z](([A-Za-z0-9)*-])|("("[^+]))*{space}
reactionKeywordSlash  {reactionKeyword}"/"{space}
thirdBodyEfficiency   {space}{floatNum}{space}"/"{space}
startReactionCoeffs   {space}"/"{space}
endReactionCoeffs     {space}"/"{space}
endPReactionCoeffs    {space}"/"{newline}
reactionCoeff         {space}{floatNum}{space}

calPerMol             {space}"CAL/MOLE"{space}
kcalPerMol            {space}"KCAL/MOLE"{space}
joulePerMol           {space}"JOULES/MOLE"{space}
kelvins               {space}"KELVINS"{space}
otherReactionsUnit    {space}("EVOLTS"|"MOLES"|"MOLECULES"){space}

cal                   {space}"CAL"{space}
kcal                  {space}"KCAL"{space}
joule                 {space}"JOUL"{space}
kjoule                {space}"KJOU"{space}
kelvin                {space}("KELV"|"KELVIN"){space}
otherReactionUnit     {space}("MOLE"|"MOLECULE"|"EVOL"|"EVOLTS"){space}


 /* ------------------------------------------------------------------------- *\
                      -----  Exclusive start states -----
 \* ------------------------------------------------------------------------- */

%option stack

%x readElements
%x readIsotopeMolW
%x readSpecies
%x readThermoAll
%x readThermoSpecieName
%x readThermoDate
%x readThermoFormula
%x readThermoPhase
%x readThermoTemps
%x readThermoFormula2
%x readThermoLineLabel1
%x readThermoCoeff1
%x readThermoLineLabel2
%x readThermoCoeff2
%x readThermoLineLabel3
%x readThermoCoeff3
%x readThermoLineLabel4
%x readReactionsUnits
%x readReactionKeyword
%x readSpecieNamePlus
%x readReactionDelimiter
%x readPDependentSpecie
%x readThirdBodyEfficiency
%x readReactionCoeffs
%x readPdepReactionCoeffs
%x readTdepSpecie
%x readReactionOrderSpecie
%x readReactionOrder
%x readReactionUnit
%x CHEMKINError

%%

%{

static const char* stateNames[30] =
{
    "reading CHEMKIN III file",
    "reading elements",
    "reading isotope molecular weight",
    "reading species",
    "reading all thermodynamic data temperatures",
    "reading thermodynamic specie name",
    "reading thermodynamic data date",
    "reading thermodynamic data specie formula",
    "reading thermodynamic data specie phase",
    "reading thermodynamic data temperatures",
    "reading thermodynamic data specie formula (supplement)",
    "reading thermodynamic data line label 1",
    "reading thermodynamic data coefficient set 1",
    "reading thermodynamic data line label 2",
    "reading thermodynamic data coefficient set 2",
    "reading thermodynamic data line label 3",
    "reading thermodynamic data coefficient set 3",
    "reading thermodynamic data line label 4",
    "reading reaction units",
    "reading reaction specie/keyword",
    "reading reaction specie",
    "reading reaction delimiter",
    "reading reaction pressure dependent specie",
    "reading third-body efficiency",
    "reading reaction coeff",
    "reading temperature dependent specie name",
    "reading reaction order specie name",
    "reading reaction order",
    "reading reaction unit",
    "error"
};

static const char* stateExpects[30] =
{
    "'ELEMENTS' or 'ELEM', 'SPECIES' or 'SPEC', 'THERMO' or 'THERMO ALL', 'REACTIONS'",
    "<elementName>, <isotopeName> / or 'END'",
    "<scalar>/",
    "<specieName> or 'END'",
    "<scalar><scalar><scalar> (3F10.0)",
    "<word> (18A1)",
    "<word> (6A1)",
    "<word><label><word><label><word><label><word><label> (4(2A1,I3))",
    "<char> (A1)",
    "<scalar><scalar><scalar> (E10.0E10.0E8.0)",
    "<word><label> (2A1,I3)",
    "'1' (I1)",
    "<scalar><scalar><scalar><scalar><scalar> (5(E15.0))",
    "'2' (I1)",
    "<scalar><scalar><scalar><scalar><scalar> (5(E15.0))",
    "'3' (I1)",
    "<scalar><scalar><scalar><scalar> (4(E15.0))",
    "'4' (I1)",
    "'CAL/MOLE', 'KCAL/MOLE', 'JOULES/MOLE', 'KELVINS', 'EVOLTS', 'MOLES' or 'MOLECULES'",
    "<word> or <label>",
    "'+'",
    "'+', '=', '<=>', '=>', '(+<word>', '<scalar> <scalar> <scalar>'",
    "<word>')'",
    "<scalar>'/'",
    "<scalar>",
    "<word>",
    "<word>",
    "<scalar>",
    "'MOLE', 'MOLECULE', 'CAL', 'KCAL', 'JOUL', 'KJOU', 'KELV', 'KELVIN', 'EVOL' or 'EVOLTS'",
    ""
};

string startError;

static const scalar RRjoule = 8.31451; // J/kg-mol-K
static const scalar RRcal = 1.987316;  // cal/g-mol-K

scalar RRreactions = RRcal;
scalar RRreaction = RRcal;

scalar allCommonT = 1000.0;

word currentElementName;
label currentElementIndex = 0;

word currentSpecieName;
label currentSpecieIndex = 0;
label nSpecieElements = 0;
List<specieElement> currentSpecieComposition(5);

scalar currentLowT = 0;
scalar currentHighT = 0;
scalar currentCommonT = 0;
chemkinIIReader::thermoPhysics::coeffArray highCpCoeffs(scalarList(7));
chemkinIIReader::thermoPhysics::coeffArray lowCpCoeffs(scalarList(7));

specieCoeffs currentSpecieCoeff;

DynamicList<specieCoeffs> lhs;
DynamicList<specieCoeffs> rhs;

scalarList ArrheniusCoeffs(3);
DynamicList<scalar> reactionCoeffs;
DynamicList<scalar> pDepReactionCoeffs;
scalarList thirdBodyEfficiencies;
label currentThirdBodyIndex = -1;

word reactionCoeffsName = word::null;
HashTable<scalarList> reactionCoeffsTable;

DynamicList<specieCoeffs> *lrhsPtr = &lhs;

reactionType rType = unknownReactionType;
reactionRateType rrType = Arrhenius;
fallOffFunctionType fofType = unknownFallOffFunctionType;
word pDependentSpecieName = word::null;
label lhsThirdBodyCounter = 0;
label rhsThirdBodyCounter = 0;

bool finishReaction = false;

%}

 /* ------------------------------------------------------------------------- *\
                            ------ Start Lexing ------
 \* ------------------------------------------------------------------------- */

 /* ------------------------------------------------------------------------- *\
    ------ Discard comments being careful to count comments
 \* ------------------------------------------------------------------------- */

<*>{space}"!".* { // Remove one line comments
    }

 /* ------------------------------------------------------------------------- *\
    ------ Read elements
 \* ------------------------------------------------------------------------- */

{elements} {
        BEGIN(readElements);
    }

<readElements>{elementName} {
        currentElementName = foamName(YYText());
        correctElementName(currentElementName);

        if (!elementIndices_.found(currentElementName))
        {
            elementIndices_.insert(currentElementName, currentElementIndex++);
            elementNames_.append(currentElementName);
        }
        else
        {
            WarningInFunction
                << "element " << currentElementName
                << " already in table." << endl;
        }
    }

<readElements>{startIsotopeMolW} {
        BEGIN(readIsotopeMolW);
    }

<readIsotopeMolW>{isotopeMolW} {
        isotopeAtomicWts_.insert(currentElementName, stringToScalar(YYText()));
        BEGIN(readElements);
    }

<readElements>{end} {
        BEGIN(INITIAL);
    }

 /* ------------------------------------------------------------------------- *\
    ------ Read species
 \* ------------------------------------------------------------------------- */

<INITIAL,readElements>{species} {
        BEGIN(readSpecies);
    }

<readSpecies>{specieName} {
        word specieName(foamName(foamSpecieString(YYText())));

        if (specieName == "THERMO")
        {
            specieNames_.shrink();
            speciesTable_ = specieNames_;
            thirdBodyEfficiencies.setSize(specieNames_.size());
            thirdBodyEfficiencies = 1.0;
            BEGIN(readThermoSpecieName);
        }
        else if (specieName == "END")
        {
            specieNames_.shrink();
            speciesTable_ = specieNames_;
            thirdBodyEfficiencies.setSize(specieNames_.size());
            thirdBodyEfficiencies = 1.0;
            BEGIN(INITIAL);
        }
        else
        {
            if (!specieIndices_.found(specieName))
            {
                specieNames_.append(specieName);
                specieIndices_.insert(specieName, currentSpecieIndex++);
            }
            else
            {
                WarningInFunction
                    << "specie " << specieName
                    << " already in table." << endl;
            }
        }
    }

 /* ------------------------------------------------------------------------- *\
    ------ Read thermo
 \* ------------------------------------------------------------------------- */

<INITIAL,readSpecies>{thermo} {
        BEGIN(readThermoSpecieName);
    }

<INITIAL,readSpecies>{thermoAll} {
        BEGIN(readThermoAll);
    }

<readThermoAll>{thermoTemp}{thermoTemp}{thermoTemp} {
        string temperaturesString(YYText());
        // scalar lowestT(stringToScalar(temperaturesString(0, 10)));
        allCommonT = stringToScalar(temperaturesString(10, 10));
        // scalar highestT(stringToScalar(temperaturesString(20, 10)));
        BEGIN(readThermoSpecieName);
    }

<readThermoSpecieName>{thermoSpecieName} {
        string specieString(foamSpecieString(YYText()));
        if (newFormat_)
        {
            specieString.replaceAll(" ", "_");
            size_t strEnd = specieString.find_last_not_of('_');
            currentSpecieName = specieString.substr(0, strEnd + 1);
        }
        else
        {
            size_t spacePos = specieString.find(' ');
            if (spacePos != string::npos)
            {
                currentSpecieName = specieString(0, spacePos);
            }
            else
            {
                currentSpecieName = specieString;
            }
        }

        BEGIN(readThermoDate);
    }

<readThermoDate>{thermoDate} {
        // string thermoDate(YYText());
        // Date is not currently used
        BEGIN(readThermoFormula);
    }

<readThermoFormula>{thermoFormula} {
        string thermoFormula(YYText());

        nSpecieElements = 0;
        currentSpecieComposition.setSize(5);

        for (int i=0; i<4; i++)
        {
            word elementName(foamName(thermoFormula(5*i, 2)));
            label nAtoms = atoi(thermoFormula(5*i + 2, 3).c_str());

            if (elementName.size() && nAtoms)
            {
                correctElementName(elementName);
                currentSpecieComposition[nSpecieElements].name() =
                    elementName;
                currentSpecieComposition[nSpecieElements++].nAtoms() = nAtoms;
            }
        }

        BEGIN(readThermoPhase);
    }

<readThermoPhase>{thermoPhase} {
        char phaseChar = YYText()[0];

        switch (phaseChar)
        {
            case 'S':
                speciePhase_.insert(currentSpecieName, solid);
            break;

            case 'L':
                speciePhase_.insert(currentSpecieName, liquid);
            break;

            case 'G':
                speciePhase_.insert(currentSpecieName, gas);
            break;
        }

        BEGIN(readThermoTemps);
    }

<readThermoTemps>{thermoLowTemp}{thermoHighTemp}{thermoCommonTemp} {
        string temperaturesString(YYText());
        currentLowT = stringToScalar(temperaturesString(0, 10));
        currentHighT = stringToScalar(temperaturesString(10, 10));
        currentCommonT = stringToScalar(temperaturesString(20, 8));

        if (currentCommonT < SMALL)
        {
            currentCommonT = allCommonT;
        }

        BEGIN(readThermoFormula2);
    }

<readThermoFormula2>{thermoFormula2} {
        string thermoFormula(YYText());
        word elementName(foamName(thermoFormula(0, 2)));
        label nAtoms = atoi(thermoFormula(2, 3).c_str());

        if
        (
            elementName.size()
         && elementName.find('0') == string::npos
         && nAtoms
        )
        {
            correctElementName(elementName);
            currentSpecieComposition[nSpecieElements].name() =
                elementName;
            currentSpecieComposition[nSpecieElements++].nAtoms() = nAtoms;
        }

        currentSpecieComposition.setSize(nSpecieElements);

        speciesCompositionTable::iterator specieCompositionIter
        (
            speciesComposition_.find(currentSpecieName)
        );

        if (specieCompositionIter != speciesComposition_.end())
        {
            speciesComposition_.erase(specieCompositionIter);
        }

        speciesComposition_.insert
        (
            currentSpecieName,
            currentSpecieComposition
        );

        BEGIN(readThermoLineLabel1);
    }

<readThermoLineLabel1>{thermoLineLabel1} {
        BEGIN(readThermoCoeff1);
    }

<readThermoCoeff1>{thermoCoeff}{thermoCoeff}{thermoCoeff}{thermoCoeff}{thermoCoeff} {
        string thermoCoeffString(YYText());

        highCpCoeffs[0] = stringToScalar(thermoCoeffString(0, 15));
        highCpCoeffs[1] = stringToScalar(thermoCoeffString(15, 15));
        highCpCoeffs[2] = stringToScalar(thermoCoeffString(30, 15));
        highCpCoeffs[3] = stringToScalar(thermoCoeffString(45, 15));
        highCpCoeffs[4] = stringToScalar(thermoCoeffString(60, 15));

        BEGIN(readThermoLineLabel2);
    }

<readThermoLineLabel2>{thermoLineLabel2} {
        BEGIN(readThermoCoeff2);
    }

<readThermoCoeff2>{thermoCoeff}{thermoCoeff}{thermoCoeff}{thermoCoeff}{thermoCoeff} {
        string thermoCoeffString(YYText());

        highCpCoeffs[5] = stringToScalar(thermoCoeffString(0, 15));
        highCpCoeffs[6] = stringToScalar(thermoCoeffString(15, 15));

        lowCpCoeffs[0] = stringToScalar(thermoCoeffString(30, 15));
        lowCpCoeffs[1] = stringToScalar(thermoCoeffString(45, 15));
        lowCpCoeffs[2] = stringToScalar(thermoCoeffString(60, 15));

        BEGIN(readThermoLineLabel3);
    }

<readThermoLineLabel3>{thermoLineLabel3} {
        BEGIN(readThermoCoeff3);
    }

<readThermoCoeff3>{thermoCoeff}{thermoCoeff}{thermoCoeff}{thermoCoeff}{thermoCoeff} {
        string thermoCoeffString(YYText());

        lowCpCoeffs[3] = stringToScalar(thermoCoeffString(0, 15));
        lowCpCoeffs[4] = stringToScalar(thermoCoeffString(15, 15));
        lowCpCoeffs[5] = stringToScalar(thermoCoeffString(30, 15));
        lowCpCoeffs[6] = stringToScalar(thermoCoeffString(45, 15));

        BEGIN(readThermoLineLabel4);
    }

<readThermoLineLabel4>{thermoLineLabel4} {

        HashPtrTable<chemkinIIReader::thermoPhysics>::iterator specieThermoIter
        (
            speciesThermo_.find(currentSpecieName)
        );

        if (specieThermoIter != speciesThermo_.end())
        {
            speciesThermo_.erase(specieThermoIter);
        }

        speciesThermo_.insert
        (
            currentSpecieName,
            new chemkinIIReader::thermoPhysics
            (
                janafThermo<perfectGas<specie>>
                (
                    specie
                    (
                        currentSpecieName,
                        1.0,
                        molecularWeight(currentSpecieComposition)
                    ),
                    currentLowT,
                    currentHighT,
                    currentCommonT,
                    highCpCoeffs,
                    lowCpCoeffs,
                    true
                ),
                transportDict_.subDict(currentSpecieName)
            )
        );

        BEGIN(readThermoSpecieName);
    }

<readThermoSpecieName>{end} {
        Reaction<chemkinIIReader::thermoPhysics>::TlowDefault = max
        (
            Reaction<chemkinIIReader::thermoPhysics>::TlowDefault,
            currentLowT
        );

        Reaction<chemkinIIReader::thermoPhysics>::ThighDefault = min
        (
            Reaction<chemkinIIReader::thermoPhysics>::ThighDefault,
            currentHighT
        );

        BEGIN(INITIAL);
    }

 /* ------------------------------------------------------------------------- *\
    ------ Read reactions
 \* ------------------------------------------------------------------------- */

<INITIAL,readThermoSpecieName>{reactions} {
        currentSpecieCoeff.stoichCoeff = 1.0;
        currentSpecieCoeff.exponent = currentSpecieCoeff.stoichCoeff;
        BEGIN(readReactionsUnits);
    }

<readReactionsUnits>{calPerMol} {
        RRreactions = RRcal;
    }

<readReactionsUnits>{kcalPerMol} {
        RRreactions = RRcal/1000.0;
    }

<readReactionsUnits>{joulePerMol} {
        RRreactions = RRjoule;
    }

<readReactionsUnits>{kelvins} {
        RRreactions = 1.0;
    }

<readReactionsUnits>{otherReactionsUnit} {
    }

<readReactionsUnits>\n {
        lineNo_++;
        BEGIN(readReactionKeyword);
    }

<readReactionKeyword>{Y} {
        currentSpecieCoeff.stoichCoeff = atof(YYText());
        currentSpecieCoeff.exponent = currentSpecieCoeff.stoichCoeff;
    }

<readReactionKeyword>{reactionKeyword} {

        word keyword(foamName(YYText()));

        HashTable<int>::iterator reactionKeywordIter
        (
            reactionKeywordTable_.find(keyword)
        );

        if (reactionKeywordIter != reactionKeywordTable_.end())
        {
            switch(reactionKeywordIter())
            {
                case duplicateReactionType:
                {
                    BEGIN(readReactionKeyword);
                    break;
                }

                case thirdBodyReactionType:
                {
                    if (rrType != Arrhenius && rrType != thirdBodyArrhenius)
                    {
                        FatalErrorInFunction
                            << "Attempt to set reaction rate type to"
                               " thirdBodyArrhenius when it is already set to "
                            << reactionRateTypeNames[rrType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if (pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "A non-pressure dependent third-body appears in"
                               " the pressure dependent reaction on line "
                            << lineNo_
                            << exit(FatalError);
                    }

                    rrType = thirdBodyArrhenius;

                    if (lrhsPtr == &lhs)
                    {
                        lhsThirdBodyCounter++;
                    }
                    else
                    {
                        rhsThirdBodyCounter++;
                    }

                    BEGIN(readReactionDelimiter);
                    break;
                }

                case plasmaMomentumTransfer:
                {
                    FatalErrorInFunction
                        << "Plasma momentum-transfer in reaction on line "
                        << lineNo_ << "not yet supported"
                        << exit(FatalError);

                    BEGIN(readReactionKeyword);
                    break;
                }

                case collisionCrossSection:
                {
                    FatalErrorInFunction
                        << "Collision cross-section in reaction on line "
                        << lineNo_ << "not yet supported"
                        << exit(FatalError);

                    BEGIN(readReactionKeyword);
                    break;
                }

                case end:
                {
                    BEGIN(INITIAL);
                    addReaction
                    (
                        lhs,
                        rhs,
                        thirdBodyEfficiencies,
                        rType,
                        rrType,
                        fofType,
                        ArrheniusCoeffs,
                        reactionCoeffsTable,
                        RRreaction
                    );
                    finishReaction = false;
                    rType = unknownReactionType;
                    rrType = Arrhenius;
                    fofType = unknownFallOffFunctionType;
                    thirdBodyEfficiencies = 1.0;
                    pDependentSpecieName = word::null;
                    lrhsPtr = &lhs;
                    break;
                }

                default:
                {
                    FatalErrorInFunction
                        << "keyword " << keyword
                        << " should be followed by parameters"
                        << " on line " << lineNo_
                        << exit(FatalError);
                }
            }
        }
        else
        {
            currentSpecieName = foamName(foamSpecieString(YYText()));

            HashTable<label>::iterator specieIndexIter
            (
                specieIndices_.find(currentSpecieName)
            );

            if (specieIndexIter != specieIndices_.end())
            {
                if (finishReaction)
                {
                    
                    if (reactionCoeffsName == "pressureDependentArrhenius")
                    {
                        reactionCoeffsTable.insert(reactionCoeffsName, pDepReactionCoeffs.shrink());
                        pDepReactionCoeffs.clear();
                    }                 

                    addReaction
                    (
                        lhs,
                        rhs,
                        thirdBodyEfficiencies,
                        rType,
                        rrType,
                        fofType,
                        ArrheniusCoeffs,
                        reactionCoeffsTable,
                        RRreaction
                    );
                    finishReaction = false;
                    rType = unknownReactionType;
                    rrType = Arrhenius;
                    fofType = unknownFallOffFunctionType;
                    thirdBodyEfficiencies = 1.0;
                    pDependentSpecieName = word::null;
                    lrhsPtr = &lhs;
                }

                currentSpecieCoeff.index = specieIndexIter();
                lrhsPtr->append(currentSpecieCoeff);

                BEGIN(readReactionDelimiter);
            }
            else
            {
                BEGIN(readSpecieNamePlus);
            }
        }
    }

<readReactionKeyword>{reactionKeywordSlash} {

        word keyword(foamName(YYText()));

        HashTable<int>::iterator reactionKeywordIter
        (
            reactionKeywordTable_.find(keyword)
        );

        if (reactionKeywordIter != reactionKeywordTable_.end())
        {
            switch(reactionKeywordIter())
            {
                case unimolecularFallOffReactionType:
                {
                    if (!pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "LOW keyword given for a unimolecular fall-off"
                               " reaction which does not contain a pressure"
                               " dependent specie" << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if (rrType == Arrhenius)
                    {
                        rrType = unimolecularFallOff;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set reaction rate type to"
                               " unimolecularFallOff when it is already set to "
                            << reactionRateTypeNames[rrType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if (fofType == unknownFallOffFunctionType)
                    {
                        fofType = Lindemann;
                    }

                    reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case chemicallyActivatedBimolecularReactionType:
                {
                    if (!pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "HIGH keyword given for a chemically"
                               " activated bimolecular reaction which does not"
                               " contain a pressure dependent specie"
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if (rrType == Arrhenius)
                    {
                        rrType = chemicallyActivatedBimolecular;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set reaction rate type to"
                               " chemicallyActivatedBimolecular when it is"
                               " already set to "
                            << reactionRateTypeNames[rrType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if (fofType == unknownFallOffFunctionType)
                    {
                        fofType = Lindemann;
                    }

                    reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case TroeReactionType:
                {
                    if (!pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "TROE keyword given for a"
                               " reaction which does not contain a pressure"
                               " dependent specie" << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if
                    (
                        fofType == unknownFallOffFunctionType
                     || fofType == Lindemann
                    )
                    {
                        fofType = Troe;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set fall-off function type to Troe"
                               " when it is already set to "
                            << fallOffFunctionNames[fofType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    reactionCoeffsName = fallOffFunctionNames[fofType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case SRIReactionType:
                {
                    if (!pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "SRI keyword given for a"
                               " reaction which does not contain a pressure"
                               " dependent specie" << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    if
                    (
                        fofType == unknownFallOffFunctionType
                     || fofType == Lindemann
                    )
                    {
                        fofType = SRI;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set fall-off function type to SRI"
                               " when it is already set to "
                            << fallOffFunctionNames[fofType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    reactionCoeffsName = fallOffFunctionNames[fofType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case LandauTellerReactionType:
                {
                    if (pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "Landau-Teller reaction rate cannot be used"
                               " for the pressure-dependent reaction on line "
                            << lineNo_
                            << exit(FatalError);
                    }

                    if (rrType == Arrhenius)
                    {
                        rrType = LandauTeller;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set reaction rate type to"
                               " LandauTeller when it is already set to "
                            << reactionRateTypeNames[rrType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    rrType = LandauTeller;
                    reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case reverseLandauTellerReactionType:
                {
                    if (pDependentSpecieName.size())
                    {
                        FatalErrorInFunction
                            << "Non-equilibrium Landau-Teller reaction rate"
                               " cannot be used"
                               " for the pressure-dependent reaction on line "
                            << lineNo_
                            << exit(FatalError);
                    }

                    if (rType != nonEquilibriumReversible)
                    {
                        FatalErrorInFunction
                            << "Reverse reaction Arrhenius coefficients not"
                               " given for reverse LandauTeller reaction."
                               " Please reorder 'REV' keyword to precede 'RLT'"
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    rrType = LandauTeller;
                    reactionCoeffsName =
                        word(reactionTypeNames[rType])
                      + reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case JanevReactionType:
                {
                    if (rrType == Arrhenius)
                    {
                        rrType = Janev;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set reaction rate type to"
                               " Janev when it is already set to "
                            << reactionRateTypeNames[rrType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case powerSeriesReactionRateType:
                {
                    if (rrType == Arrhenius)
                    {
                        rrType = powerSeries;
                    }
                    else
                    {
                        FatalErrorInFunction
                            << "Attempt to set reaction rate type to"
                               " powerSeries when it is already set to "
                            << reactionRateTypeNames[rrType]
                            << " on line " << lineNo_
                            << exit(FatalError);
                    }

                    reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case radiationActivatedReactionType:
                {
                    FatalErrorInFunction
                        << "Radiation activated reaction on line "
                        << lineNo_ << "not yet supported"
                        << exit(FatalError);
                    // reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case energyLossReactionType:
                {
                    FatalErrorInFunction
                        << "Energy loss in reaction on line "
                        << lineNo_ << "not yet supported"
                        << exit(FatalError);
                    // reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case nonEquilibriumReversibleReactionType:
                {
                    rType = nonEquilibriumReversible;
                    reactionCoeffsName = reactionTypeNames[rType];
                    BEGIN(readReactionCoeffs);
                    break;
                }

                case speciesOrderForward:
                {
                    lrhsPtr = &lhs;
                    BEGIN(readReactionOrderSpecie);
                    break;
                }

                case speciesOrderReverse:
                {
                    lrhsPtr = &rhs;
                    BEGIN(readReactionOrderSpecie);
                    break;
                }

                case UnitsOfReaction:
                {
                    BEGIN(readReactionUnit);
                    break;
                }

                case pressureDependentReactionType:
                {
                    if (rrType == Arrhenius)
                    {
                        rrType = pressureDependentArrhenius;
                    }
                    reactionCoeffsName = reactionRateTypeNames[rrType];
                    BEGIN(readPdepReactionCoeffs);
                    break;
                }

                default:
                {
                    FatalErrorInFunction
                        << "unknown reaction keyword " << keyword
                        << " on line " << lineNo_
                        << exit(FatalError);
                }
            }
        }
        else
        {
            HashTable<label>::iterator specieIndexIter
            (
                specieIndices_.find(keyword)
            );

            if (specieIndexIter != specieIndices_.end())
            {
                currentThirdBodyIndex = specieIndexIter();
            }
            else
            {
                FatalErrorInFunction
                    << "unknown third-body specie " << keyword
                    << " on line " << lineNo_ << nl
                    << "Valid species are : " << nl
                    << specieIndices_.toc() << endl
                    << exit(FatalError);
            }

            BEGIN(readThirdBodyEfficiency);
        }
    }

<readSpecieNamePlus>"+" {
        currentSpecieName += "+";

        HashTable<label>::iterator specieIndexIter
        (
            specieIndices_.find(currentSpecieName)
        );

        if (specieIndexIter != specieIndices_.end())
        {
            if (finishReaction)
            {
                addReaction
                (
                    lhs,
                    rhs,
                    thirdBodyEfficiencies,
                    rType,
                    rrType,
                    fofType,
                    ArrheniusCoeffs,
                    reactionCoeffsTable,
                    RRreaction
                );
                finishReaction = false;
                rType = unknownReactionType;
                rrType = Arrhenius;
                fofType = unknownFallOffFunctionType;
                thirdBodyEfficiencies = 1.0;
                pDependentSpecieName = word::null;
                lrhsPtr = &lhs;
            }

            currentSpecieCoeff.index = specieIndexIter();
            lrhsPtr->append(currentSpecieCoeff);

            BEGIN(readReactionDelimiter);
        }
        else
        {
            FatalErrorInFunction
                << "unknown specie " << currentSpecieName
                << " on line " << lineNo_ << nl
                << "Valid species are : " << nl
                << specieIndices_.toc() << endl
                << exit(FatalError);
        }
    }

<readReactionDelimiter>{specieDelimiter} {
        currentSpecieCoeff.stoichCoeff = 1.0;
        currentSpecieCoeff.exponent = currentSpecieCoeff.stoichCoeff;
        BEGIN(readReactionKeyword);
    }

<readReactionDelimiter>{irreversibleReactionDelimiter} {
        currentSpecieCoeff.stoichCoeff = 1.0;
        currentSpecieCoeff.exponent = currentSpecieCoeff.stoichCoeff;
        rType = irreversible;
        lrhsPtr = &rhs;
        BEGIN(readReactionKeyword);
    }

<readReactionDelimiter>{reversibleReactionDelimiter} {
        currentSpecieCoeff.stoichCoeff = 1.0;
        currentSpecieCoeff.exponent = currentSpecieCoeff.stoichCoeff;
        rType = reversible;
        lrhsPtr = &rhs;
        BEGIN(readReactionKeyword);
    }

<readReactionDelimiter>& {
    }

<readReactionDelimiter>{reactionCoeffs} {
        string reactionCoeffsString(YYText());
        reactionCoeffsString.replaceAll("d", "e");
        reactionCoeffsString.replaceAll("D", "e");
        IStringStream reactionCoeffsStream(reactionCoeffsString);
        reactionCoeffsStream.lineNumber() = lineNo_;

        reactionCoeffsStream
            >> ArrheniusCoeffs[0]
            >> ArrheniusCoeffs[1]
            >> ArrheniusCoeffs[2];

        finishReaction = true;
        currentSpecieCoeff.stoichCoeff = 1.0;
        currentSpecieCoeff.exponent = currentSpecieCoeff.stoichCoeff;
        RRreaction = RRreactions;

        if (lhsThirdBodyCounter || rhsThirdBodyCounter)
        {
            if (!lhsThirdBodyCounter || !rhsThirdBodyCounter)
            {
                FatalErrorInFunction
                    << "Third body not present on both sides of reaction"
                       " on line " << lineNo_
                << exit(FatalError);
            }

            if (lhsThirdBodyCounter != 1)
            {
                FatalErrorInFunction
                    << "More than 1 third body present on l.h.s. side"
                       " of reaction on line " << lineNo_
                    << exit(FatalError);
            }

            if (rhsThirdBodyCounter != 1)
            {
                FatalErrorInFunction
                    << "More than 1 third body present on r.h.s. side"
                       " of reaction on line " << lineNo_
                    << exit(FatalError);
            }

            lhsThirdBodyCounter = 0;
            rhsThirdBodyCounter = 0;
        }

        BEGIN(readReactionKeyword);
    }


<readThirdBodyEfficiency>{thirdBodyEfficiency} {
        thirdBodyEfficiencies[currentThirdBodyIndex] = stringToScalar(YYText());
        BEGIN(readReactionKeyword);
    }

<readReactionDelimiter>{startPDependentSpecie} {
        BEGIN(readPDependentSpecie);
    }

<readPDependentSpecie>{pDependentSpecie} {

        word rhsPDependentSpecieName = pDependentSpecieName;
        pDependentSpecieName =
            foamName(foamSpecieString(YYText()));
        pDependentSpecieName =
            pDependentSpecieName(0, pDependentSpecieName.size() - 1);

        if (rrType == thirdBodyArrhenius)
        {
            FatalErrorInFunction
                << "The pressure-dependent third-body '"
                << pDependentSpecieName
                << "' is given in non-pressure-dependent third-body reaction"
                << " on line " << lineNo_
                << exit(FatalError);
        }

        if (lrhsPtr == &lhs)
        {
            lhsThirdBodyCounter++;
        }
        else
        {
            if (pDependentSpecieName != rhsPDependentSpecieName)
            {
                FatalErrorInFunction
                    << "The third-body reactant '"
                    << pDependentSpecieName
                    << "' is not the same as the third-body product '"
                    << rhsPDependentSpecieName
                    << "' in pressure-dependent reaction on line " << lineNo_
                    << exit(FatalError);
            }

            rhsThirdBodyCounter++;
        }

        if (pDependentSpecieName != "M")
        {
            HashTable<label>::iterator specieIndexIter
            (
                specieIndices_.find(pDependentSpecieName)
            );

            if (specieIndexIter != specieIndices_.end())
            {
                thirdBodyEfficiencies = 0.0;
                thirdBodyEfficiencies[specieIndexIter()] = 1.0;
            }
            else
            {
                FatalErrorInFunction
                    << "unknown third-body specie " << pDependentSpecieName
                    << " on line " << lineNo_ << nl
                    << "Valid species are : " << nl
                    << specieIndices_.toc() << endl
                    << exit(FatalError);
            }
        }

        BEGIN(readReactionDelimiter);
    }

<readReactionCoeffs>{reactionCoeff} {
        reactionCoeffs.append(stringToScalar(YYText()));
    }

<readReactionCoeffs>{endReactionCoeffs} {
        reactionCoeffsTable.insert(reactionCoeffsName, reactionCoeffs.shrink());
        reactionCoeffs.clear();
        BEGIN(readReactionKeyword);
    }

<readPdepReactionCoeffs>{reactionCoeff} {
        pDepReactionCoeffs.append(stringToScalar(YYText()));
    }

<readPdepReactionCoeffs>{endPReactionCoeffs} {
        BEGIN(readReactionKeyword);
    }

<readTdepSpecie>{specieName} {
        word specieName(foamName(foamSpecieString(YYText())));
        FatalErrorInFunction
            << "Temperature-dependent reaction on line "
            << lineNo_ << "not yet supported"
            << exit(FatalError);
        BEGIN(readReactionKeyword);
    }

<readReactionOrderSpecie>{specieName} {
        currentSpecieName =
            word(foamName(foamSpecieString(YYText())));

        HashTable<label>::iterator specieIndexIter
        (
            specieIndices_.find(currentSpecieName)
        );

        if (specieIndexIter != specieIndices_.end())
        {
            currentSpecieIndex = specieIndexIter();
        }
        else
        {
            FatalErrorInFunction
                << "unknown specie " << currentSpecieName
                << " given in reaction-order specification"
                << " on line " << lineNo_ << nl
                << "Valid species are : " << nl
                << specieIndices_.toc() << endl
                << exit(FatalError);
        }

        BEGIN(readReactionOrder);
    }

<readReactionOrder>{reactionCoeff}{endReactionCoeffs} {

        DynamicList<specieCoeffs>& lrhs = *lrhsPtr;

        bool found = false;

        forAll(lrhs, i)
        {
            if (lrhs[i].index == currentSpecieIndex)
            {
                lrhs[i].exponent = stringToScalar(YYText());
                found = true;
                break;
            }
        }

        if (!found)
        {
            word side("l.h.s.");

            if (lrhsPtr == &rhs)
            {
                side = "r.h.s.";
            }

            FatalErrorInFunction
                << "Specie " << currentSpecieName
                << " on line " << lineNo_
                << " not present in " << side << " of reaction " << nl << lrhs
                << exit(FatalError);
        }

        BEGIN(readReactionKeyword);
    }

<readReactionUnit>{cal} {
        RRreaction = RRcal;
    }

<readReactionUnit>{kcal} {
        RRreaction = RRcal/1000.0;
    }

<readReactionUnit>{joule} {
        RRreaction = RRjoule;
    }

<readReactionUnit>{kjoule} {
        RRreaction = RRjoule/1000.0;
    }

<readReactionUnit>{kelvin} {
        RRreaction = 1.0;
    }

<readReactionUnit>{otherReactionsUnit} {
    }

<readReactionUnit>{endReactionCoeffs} {
        BEGIN(readReactionKeyword);
    }


 /* ------------------ Ignore remaining space and \n s. --------------------- */

<*>\n      { lineNo_++; }

 /* ------ Ignore remaining space and \n s.  Any other characters are errors. */

<*>. {
        startError = YYText();
        yy_push_state(CHEMKINError);
    }

<CHEMKINError>.* {
        yy_pop_state();
        FatalErrorInFunction
            << "while " << stateNames[YY_START] << " on line " << lineNo_ << nl
            << "    expected " << stateExpects[YY_START]
            << " but found '" << startError << YYText() << "'"
            << exit(FatalError);
    }


 /*  ------------------------ On EOF terminate. ----------------------------  */

<<EOF>> {
        if (finishReaction)
        {
            addReaction
            (
                lhs,
                rhs,
                thirdBodyEfficiencies,
                rType,
                rrType,
                fofType,
                ArrheniusCoeffs,
                reactionCoeffsTable,
                RRreaction
            );
        }
        yyterminate();
    }
%%

 /* ------------------------------------------------------------------------- *\
    ------ End of chemkinIILexer.L
 \* ------------------------------------------------------------------------- */

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

\*---------------------------------------------------------------------------*/

#include "chemkinIIReader.H"
#include "IFstream.H"
#include "atomicWeights.H"
#include "ReactionProxy.H"
#include "IrreversibleReaction.H"
#include "ReversibleReaction.H"
#include "NonEquilibriumReversibleReaction.H"
#include "ArrheniusReactionRate.H"
#include "thirdBodyArrheniusReactionRate.H"
#include "pressureDependentArrheniusReactionRate.H"
#include "FallOffReactionRate.H"
#include "ChemicallyActivatedReactionRate.H"
#include "LindemannFallOffFunction.H"
#include "TroeFallOffFunction.H"
#include "SRIFallOffFunction.H"
#include "LandauTellerReactionRate.H"
#include "JanevReactionRate.H"
#include "powerSeriesReactionRate.H"
#include "addToRunTimeSelectionTable.H"


/* * * * * * * * * * * * * * Private Member Functions  * * * * * * * * * * * */

const char* Foam::chemkinIIReader::reactionTypeNames[4] =
{
    "irreversible",
    "reversible",
    "nonEquilibriumReversible",
    "unknownReactionType"
};

const char* Foam::chemkinIIReader::reactionRateTypeNames[9] =
{
    "Arrhenius",
    "pressureDependentArrhenius",
    "thirdBodyArrhenius",
    "unimolecularFallOff",
    "chemicallyActivatedBimolecular",
    "LandauTeller",
    "Janev",
    "powerSeries",
    "unknownReactionRateType"
};

const char* Foam::chemkinIIReader::fallOffFunctionNames[4] =
{
    "Lindemann",
    "Troe",
    "SRI",
    "unknownFallOffFunctionType"
};

void Foam::chemkinIIReader::initReactionKeywordTable()
{
    reactionKeywordTable_.insert("M", thirdBodyReactionType);
    reactionKeywordTable_.insert("LOW", unimolecularFallOffReactionType);
    reactionKeywordTable_.insert
    (
        "HIGH",
        chemicallyActivatedBimolecularReactionType
    );
    reactionKeywordTable_.insert("TROE", TroeReactionType);
    reactionKeywordTable_.insert("SRI", SRIReactionType);
    reactionKeywordTable_.insert("LT", LandauTellerReactionType);
    reactionKeywordTable_.insert("RLT", reverseLandauTellerReactionType);
    reactionKeywordTable_.insert("JAN", JanevReactionType);
    reactionKeywordTable_.insert("FIT1", powerSeriesReactionRateType);
    reactionKeywordTable_.insert("HV", radiationActivatedReactionType);
    reactionKeywordTable_.insert("TDEP", speciesTempReactionType);
    reactionKeywordTable_.insert("EXCI", energyLossReactionType);
    reactionKeywordTable_.insert("MOME", plasmaMomentumTransfer);
    reactionKeywordTable_.insert("XSMI", collisionCrossSection);
    reactionKeywordTable_.insert("REV", nonEquilibriumReversibleReactionType);
    reactionKeywordTable_.insert("DUPLICATE", duplicateReactionType);
    reactionKeywordTable_.insert("DUP", duplicateReactionType);
    reactionKeywordTable_.insert("FORD", speciesOrderForward);
    reactionKeywordTable_.insert("RORD", speciesOrderReverse);
    reactionKeywordTable_.insert("UNITS", UnitsOfReaction);
    reactionKeywordTable_.insert("PLOG", pressureDependentReactionType);
    reactionKeywordTable_.insert("END", end);
}


Foam::scalar Foam::chemkinIIReader::molecularWeight
(
    const List<specieElement>& specieComposition
) const
{
    scalar molWt = 0.0;

    forAll(specieComposition, i)
    {
        label nAtoms = specieComposition[i].nAtoms();
        const word& elementName = specieComposition[i].name();

        if (isotopeAtomicWts_.found(elementName))
        {
            molWt += nAtoms*isotopeAtomicWts_[elementName];
        }
        else if (atomicWeights.found(elementName))
        {
            molWt += nAtoms*atomicWeights[elementName];
        }
        else
        {
            FatalErrorInFunction
                << "Unknown element " << elementName
                << " on line " << lineNo_-1 << nl
                << "    specieComposition: " << specieComposition
                << exit(FatalError);
        }
    }

    return molWt;
}


void Foam::chemkinIIReader::checkCoeffs
(
    const scalarList& reactionCoeffs,
    const char* reactionRateName,
    const label nCoeffs
) const
{
    if (reactionCoeffs.size() != nCoeffs)
    {
        FatalErrorInFunction
            << "Wrong number of coefficients for the " << reactionRateName
            << " rate expression on line "
            << lineNo_-1 << ", should be "
            << nCoeffs << " but " << reactionCoeffs.size() << " supplied." << nl
            << "Coefficients are "
            << reactionCoeffs << nl
            << exit(FatalError);
    }
}

template<class ReactionRateType>
void Foam::chemkinIIReader::addReactionType
(
    const reactionType rType,
    DynamicList<specieCoeffs>& lhs,
    DynamicList<specieCoeffs>& rhs,
    const ReactionRateType& rr
)
{
    switch (rType)
    {
        case irreversible:
        {
            reactions_.append
            (
                new IrreversibleReaction<thermoPhysics, ReactionRateType>
                (
                    ReactionProxy<thermoPhysics>
                    (
                        speciesTable_,
                        lhs.shrink(),
                        rhs.shrink(),
                        speciesThermo_
                    ),
                    rr
                )
            );
        }
        break;

        case reversible:
        {
            reactions_.append
            (
                new ReversibleReaction<thermoPhysics, ReactionRateType>
                (
                    ReactionProxy<thermoPhysics>
                    (
                        speciesTable_,
                        lhs.shrink(),
                        rhs.shrink(),
                        speciesThermo_
                    ),
                    rr
                )
            );
        }
        break;

        default:

            if (rType < 3)
            {
                FatalErrorInFunction
                    << "Reaction type " << reactionTypeNames[rType]
                    << " on line " << lineNo_-1
                    << " not handled by this function"
                    << exit(FatalError);
            }
            else
            {
                FatalErrorInFunction
                    << "Unknown reaction type " << rType
                    << " on line " << lineNo_-1
                    << exit(FatalError);
            }
    }
}

template<template<class, class> class PressureDependencyType>
void Foam::chemkinIIReader::addPressureDependentReaction
(
    const reactionType rType,
    const fallOffFunctionType fofType,
    DynamicList<specieCoeffs>& lhs,
    DynamicList<specieCoeffs>& rhs,
    const scalarList& efficiencies,
    const scalarList& k0Coeffs,
    const scalarList& kInfCoeffs,
    const HashTable<scalarList>& reactionCoeffsTable,
    const scalar Afactor0,
    const scalar AfactorInf,
    const scalar RR
)
{
    checkCoeffs(k0Coeffs, "k0", 3);
    checkCoeffs(kInfCoeffs, "kInf", 3);

    switch (fofType)
    {
        case Lindemann:
        {
            addReactionType
            (
                rType,
                lhs, rhs,
                PressureDependencyType
                    <ArrheniusReactionRate, LindemannFallOffFunction>
                (
                    ArrheniusReactionRate
                    (
                        Afactor0*k0Coeffs[0],
                        k0Coeffs[1],
                        k0Coeffs[2]/RR
                    ),
                    ArrheniusReactionRate
                    (
                        AfactorInf*kInfCoeffs[0],
                        kInfCoeffs[1],
                        kInfCoeffs[2]/RR
                    ),
                    LindemannFallOffFunction(),
                    thirdBodyEfficiencies(speciesTable_, efficiencies)
                )
            );
            break;
        }
        case Troe:
        {
            scalarList TroeCoeffs
            (
                reactionCoeffsTable[fallOffFunctionNames[fofType]]
            );

            if (TroeCoeffs.size() != 4 && TroeCoeffs.size() != 3)
            {
                FatalErrorInFunction
                    << "Wrong number of coefficients for Troe rate expression"
                       " on line " << lineNo_-1 << ", should be 3 or 4 but "
                    << TroeCoeffs.size() << " supplied." << nl
                    << "Coefficients are "
                    << TroeCoeffs << nl
                    << exit(FatalError);
            }

            if (TroeCoeffs.size() == 3)
            {
                TroeCoeffs.setSize(4);
                TroeCoeffs[3] = great;
            }

            addReactionType
            (
                rType,
                lhs, rhs,
                PressureDependencyType
                    <ArrheniusReactionRate, TroeFallOffFunction>
                (
                    ArrheniusReactionRate
                    (
                        Afactor0*k0Coeffs[0],
                        k0Coeffs[1],
                        k0Coeffs[2]/RR
                    ),
                    ArrheniusReactionRate
                    (
                        AfactorInf*kInfCoeffs[0],
                        kInfCoeffs[1],
                        kInfCoeffs[2]/RR
                    ),
                    TroeFallOffFunction
                    (
                        TroeCoeffs[0],
                        TroeCoeffs[1],
                        TroeCoeffs[2],
                        TroeCoeffs[3]
                    ),
                    thirdBodyEfficiencies(speciesTable_, efficiencies)
                )
            );
            break;
        }
        case SRI:
        {
            scalarList SRICoeffs
            (
                reactionCoeffsTable[fallOffFunctionNames[fofType]]
            );

            if (SRICoeffs.size() != 5 && SRICoeffs.size() != 3)
            {
                FatalErrorInFunction
                    << "Wrong number of coefficients for SRI rate expression"
                       " on line " << lineNo_-1 << ", should be 3 or 5 but "
                    << SRICoeffs.size() << " supplied." << nl
                    << "Coefficients are "
                    << SRICoeffs << nl
                    << exit(FatalError);
            }

            if (SRICoeffs.size() == 3)
            {
                SRICoeffs.setSize(5);
                SRICoeffs[3] = 1.0;
                SRICoeffs[4] = 0.0;
            }

            addReactionType
            (
                rType,
                lhs, rhs,
                PressureDependencyType
                    <ArrheniusReactionRate, SRIFallOffFunction>
                (
                    ArrheniusReactionRate
                    (
                        Afactor0*k0Coeffs[0],
                        k0Coeffs[1],
                        k0Coeffs[2]/RR
                    ),
                    ArrheniusReactionRate
                    (
                        AfactorInf*kInfCoeffs[0],
                        kInfCoeffs[1],
                        kInfCoeffs[2]/RR
                    ),
                    SRIFallOffFunction
                    (
                        SRICoeffs[0],
                        SRICoeffs[1],
                        SRICoeffs[2],
                        SRICoeffs[3],
                        SRICoeffs[4]
                    ),
                    thirdBodyEfficiencies(speciesTable_, efficiencies)
                )
            );
            break;
        }
        default:
        {
            FatalErrorInFunction
                << "Fall-off function type "
                << fallOffFunctionNames[fofType]
                << " on line " << lineNo_-1
                << " not implemented"
                << exit(FatalError);
        }
    }
}


void Foam::chemkinIIReader::addReaction
(
    DynamicList<specieCoeffs>& lhs,
    DynamicList<specieCoeffs>& rhs,
    const scalarList& efficiencies,
    const reactionType rType,
    const reactionRateType rrType,
    const fallOffFunctionType fofType,
    const scalarList& ArrheniusCoeffs,
    HashTable<scalarList>& reactionCoeffsTable,
    const scalar RR
)
{
    checkCoeffs(ArrheniusCoeffs, "Arrhenius", 3);

    scalarList nAtoms(elementNames_.size(), 0.0);

    forAll(lhs, i)
    {
        const List<specieElement>& specieComposition =
            speciesComposition_[speciesTable_[lhs[i].index]];

        forAll(specieComposition, j)
        {
            label elementi = elementIndices_[specieComposition[j].name()];
            nAtoms[elementi] +=
                lhs[i].stoichCoeff*specieComposition[j].nAtoms();
        }
    }

    forAll(rhs, i)
    {
        const List<specieElement>& specieComposition =
            speciesComposition_[speciesTable_[rhs[i].index]];

        forAll(specieComposition, j)
        {
            label elementi = elementIndices_[specieComposition[j].name()];
            nAtoms[elementi] -=
                rhs[i].stoichCoeff*specieComposition[j].nAtoms();
        }
    }


    // Calculate the unit conversion factor for the A coefficient
    // for the change from mol/cm^3 to kmol/m^3 concentration units
    const scalar concFactor = 0.001;
    scalar sumExp = 0.0;
    forAll(lhs, i)
    {
        sumExp += lhs[i].exponent;
    }
    scalar Afactor = pow(concFactor, sumExp - 1.0);

    scalar AfactorRev = Afactor;

    if (rType == nonEquilibriumReversible)
    {
        sumExp = 0.0;
        forAll(rhs, i)
        {
            sumExp += rhs[i].exponent;
        }
        AfactorRev = pow(concFactor, sumExp - 1.0);
    }

    switch (rrType)
    {
        case Arrhenius:
        {
            if (rType == nonEquilibriumReversible)
            {
                const scalarList& reverseArrheniusCoeffs =
                    reactionCoeffsTable[reactionTypeNames[rType]];

                checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);

                reactions_.append
                (
                    new NonEquilibriumReversibleReaction
                    <
                        thermoPhysics,
                        ArrheniusReactionRate
                    >
                    (
                        ReactionProxy<thermoPhysics>
                        (
                            speciesTable_,
                            lhs.shrink(),
                            rhs.shrink(),
                            speciesThermo_
                        ),
                        ArrheniusReactionRate
                        (
                            Afactor*ArrheniusCoeffs[0],
                            ArrheniusCoeffs[1],
                            ArrheniusCoeffs[2]/RR
                        ),
                        ArrheniusReactionRate
                        (
                            AfactorRev*reverseArrheniusCoeffs[0],
                            reverseArrheniusCoeffs[1],
                            reverseArrheniusCoeffs[2]/RR
                        )
                    )
                );
            }
            else
            {
                addReactionType
                (
                    rType,
                    lhs, rhs,
                    ArrheniusReactionRate
                    (
                        Afactor*ArrheniusCoeffs[0],
                        ArrheniusCoeffs[1],
                        ArrheniusCoeffs[2]/RR
                    )
                );
            }
            break;
        }
        case thirdBodyArrhenius:
        {
            if (rType == nonEquilibriumReversible)
            {
                const scalarList& reverseArrheniusCoeffs =
                    reactionCoeffsTable[reactionTypeNames[rType]];

                checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);

                reactions_.append
                (
                    new NonEquilibriumReversibleReaction
                    <
                        thermoPhysics,
                        thirdBodyArrheniusReactionRate
                    >
                    (
                        ReactionProxy<thermoPhysics>
                        (
                            speciesTable_,
                            lhs.shrink(),
                            rhs.shrink(),
                            speciesThermo_
                        ),
                        thirdBodyArrheniusReactionRate
                        (
                            Afactor*concFactor*ArrheniusCoeffs[0],
                            ArrheniusCoeffs[1],
                            ArrheniusCoeffs[2]/RR,
                            thirdBodyEfficiencies(speciesTable_, efficiencies)
                        ),
                        thirdBodyArrheniusReactionRate
                        (
                            AfactorRev*concFactor*reverseArrheniusCoeffs[0],
                            reverseArrheniusCoeffs[1],
                            reverseArrheniusCoeffs[2]/RR,
                            thirdBodyEfficiencies(speciesTable_, efficiencies)
                        )
                    )
                );
            }
            else
            {
                addReactionType
                (
                    rType,
                    lhs, rhs,
                    thirdBodyArrheniusReactionRate
                    (
                        Afactor*concFactor*ArrheniusCoeffs[0],
                        ArrheniusCoeffs[1],
                        ArrheniusCoeffs[2]/RR,
                        thirdBodyEfficiencies(speciesTable_, efficiencies)
                    )
                );
            }
            break;
        }
        case unimolecularFallOff:
        {
            addPressureDependentReaction<FallOffReactionRate>
            (
                rType,
                fofType,
                lhs,
                rhs,
                efficiencies,
                reactionCoeffsTable[reactionRateTypeNames[rrType]],
                ArrheniusCoeffs,
                reactionCoeffsTable,
                concFactor*Afactor,
                Afactor,
                RR
            );
            break;
        }
        case chemicallyActivatedBimolecular:
        {
            addPressureDependentReaction<ChemicallyActivatedReactionRate>
            (
                rType,
                fofType,
                lhs,
                rhs,
                efficiencies,
                ArrheniusCoeffs,
                reactionCoeffsTable[reactionRateTypeNames[rrType]],
                reactionCoeffsTable,
                Afactor,
                Afactor/concFactor,
                RR
            );
            break;
        }
        case LandauTeller:
        {
            const scalarList& LandauTellerCoeffs =
                reactionCoeffsTable[reactionRateTypeNames[rrType]];
            checkCoeffs(LandauTellerCoeffs, "Landau-Teller", 2);

            if (rType == nonEquilibriumReversible)
            {
                const scalarList& reverseArrheniusCoeffs =
                    reactionCoeffsTable[reactionTypeNames[rType]];
                checkCoeffs(reverseArrheniusCoeffs, "reverse Arrhenius", 3);

                const scalarList& reverseLandauTellerCoeffs =
                    reactionCoeffsTable
                    [
                        word(reactionTypeNames[rType])
                      + reactionRateTypeNames[rrType]
                    ];
                checkCoeffs(LandauTellerCoeffs, "reverse Landau-Teller", 2);

                reactions_.append
                (
                    new NonEquilibriumReversibleReaction
                    <
                        thermoPhysics,
                        LandauTellerReactionRate
                    >
                    (
                        ReactionProxy<thermoPhysics>
                        (
                            speciesTable_,
                            lhs.shrink(),
                            rhs.shrink(),
                            speciesThermo_
                        ),
                        LandauTellerReactionRate
                        (
                            Afactor*ArrheniusCoeffs[0],
                            ArrheniusCoeffs[1],
                            ArrheniusCoeffs[2]/RR,
                            LandauTellerCoeffs[0],
                            LandauTellerCoeffs[1]
                        ),
                        LandauTellerReactionRate
                        (
                            AfactorRev*reverseArrheniusCoeffs[0],
                            reverseArrheniusCoeffs[1],
                            reverseArrheniusCoeffs[2]/RR,
                            reverseLandauTellerCoeffs[0],
                            reverseLandauTellerCoeffs[1]
                        )
                    )
                );
            }
            else
            {
                addReactionType
                (
                    rType,
                    lhs, rhs,
                    LandauTellerReactionRate
                    (
                        Afactor*ArrheniusCoeffs[0],
                        ArrheniusCoeffs[1],
                        ArrheniusCoeffs[2]/RR,
                        LandauTellerCoeffs[0],
                        LandauTellerCoeffs[1]
                    )
                );
            }
            break;
        }
        case Janev:
        {
            const scalarList& JanevCoeffs =
                reactionCoeffsTable[reactionRateTypeNames[rrType]];

            checkCoeffs(JanevCoeffs, "Janev", 9);

            addReactionType
            (
                rType,
                lhs, rhs,
                JanevReactionRate
                (
                    Afactor*ArrheniusCoeffs[0],
                    ArrheniusCoeffs[1],
                    ArrheniusCoeffs[2]/RR,
                    FixedList<scalar, 9>(JanevCoeffs)
                )
            );
            break;
        }
        case powerSeries:
        {
            const scalarList& powerSeriesCoeffs =
                reactionCoeffsTable[reactionRateTypeNames[rrType]];

            checkCoeffs(powerSeriesCoeffs, "power-series", 4);

            addReactionType
            (
                rType,
                lhs, rhs,
                powerSeriesReactionRate
                (
                    Afactor*ArrheniusCoeffs[0],
                    ArrheniusCoeffs[1],
                    ArrheniusCoeffs[2]/RR,
                    FixedList<scalar, 4>(powerSeriesCoeffs)
                )
            );
            break;
        }
        case pressureDependentArrhenius:
        {
            const scalarList& pressureCoeffs = reactionCoeffsTable[reactionRateTypeNames[rrType]];
            List<Tuple2<scalar, scalar>> pArrheniusA(0.0);
            List<Tuple2<scalar, scalar>> pArrheniusbeta(0.0);
            List<Tuple2<scalar, scalar>> pArrheniusTA(0.0);
            for(int p = 0; p < pressureCoeffs.size(); p += 4)
            {
                pArrheniusA.append
                (
                    Tuple2<scalar, scalar>
                    (
                        log(pressureCoeffs[p]*1.0e5), 
                        log(Afactor*pressureCoeffs[p+1])
                    )
                );
                pArrheniusbeta.append
                (
                    Tuple2<scalar, scalar>
                    (
                        log(pressureCoeffs[p]*1.0e5), 
                        pressureCoeffs[p+2]
                    )
                );
                pArrheniusTA.append
                (
                    Tuple2<scalar, scalar>
                    (
                        log(pressureCoeffs[p]*1.0e5), 
                        pressureCoeffs[p+3]/RR
                    )
                );
            }
            Function1s::Table<scalar> pFuncA
            (
                "A",
                Function1s::tableBase::boundsHandling::clamp,
                linearInterpolationWeights::typeName,
                pArrheniusA
            );

            Function1s::Table<scalar> pFuncBeta
            (
                "beta",
                Function1s::tableBase::boundsHandling::clamp,
                linearInterpolationWeights::typeName,
                pArrheniusbeta
            );

            Function1s::Table<scalar> pFuncTA
            (
                "TA",
                Function1s::tableBase::boundsHandling::clamp,
                linearInterpolationWeights::typeName,
                pArrheniusTA
            );
            addReactionType
            (
                rType,
                lhs, rhs,
                pressureDependentArrheniusReactionRate
                (
                    pFuncA,
                    pFuncBeta,
                    pFuncTA
                )
            );
            break;
        }
        case unknownReactionRateType:
        {
            FatalErrorInFunction
                << "Internal error on line " << lineNo_-1
                << ": reaction rate type has not been set"
                << exit(FatalError);
            break;
        }
        default:
        {
            FatalErrorInFunction
                << "Reaction rate type " << reactionRateTypeNames[rrType]
                << " on line " << lineNo_-1
                << " not implemented"
                << exit(FatalError);
        }
    }


    forAll(nAtoms, i)
    {
        if (mag(nAtoms[i]) > imbalanceTol_)
        {
            FatalErrorInFunction
                << "Elemental imbalance of " << mag(nAtoms[i])
                << " in " << elementNames_[i]
                << " in reaction" << nl
                << reactions_.last() << nl
                << " on line " << lineNo_-1
                << exit(FatalError);
        }
    }

    lhs.clear();
    rhs.clear();
    reactionCoeffsTable.clear();
}


void Foam::chemkinIIReader::read
(
    const fileName& CHEMKINFileName,
    const fileName& thermoFileName,
    const fileName& transportFileName
)
{
    Reaction<thermoPhysics>::TlowDefault = 0;
    Reaction<thermoPhysics>::ThighDefault = great;

    transportDict_.read(IFstream(transportFileName)());

    if (thermoFileName != fileName::null)
    {
        std::ifstream thermoStream(thermoFileName.c_str());

        if (!thermoStream)
        {
            FatalErrorInFunction
                << "file " << thermoFileName << " not found"
                << exit(FatalError);
        }

        yy_buffer_state* bufferPtr(yy_create_buffer(&thermoStream, yyBufSize));
        yy_switch_to_buffer(bufferPtr);

        while (lex() != 0)
        {}

        yy_delete_buffer(bufferPtr);

        lineNo_ = 1;
    }

    std::ifstream CHEMKINStream(CHEMKINFileName.c_str());

    if (!CHEMKINStream)
    {
        FatalErrorInFunction
            << "file " << CHEMKINFileName << " not found"
            << exit(FatalError);
    }

    yy_buffer_state* bufferPtr(yy_create_buffer(&CHEMKINStream, yyBufSize));
    yy_switch_to_buffer(bufferPtr);

    initReactionKeywordTable();

    while (lex() != 0)
    {}

    yy_delete_buffer(bufferPtr);
}


// * * * * * * * * * * * * * * * * Constructor * * * * * * * * * * * * * * * //

Foam::chemkinIIReader::chemkinIIReader
(
    speciesTable& species,
    const fileName& CHEMKINFileName,
    const fileName& transportFileName,
    const fileName& thermoFileName,
    const bool newFormat
)
:
    lineNo_(1),
    specieNames_(10),
    speciesTable_(species),
    newFormat_(newFormat),
    imbalanceTol_(rootSmall)
{
    read(CHEMKINFileName, thermoFileName, transportFileName);
}


Foam::chemkinIIReader::chemkinIIReader
(
    const dictionary& thermoDict,
    speciesTable& species
)
:
    lineNo_(1),
    specieNames_(10),
    speciesTable_(species),
    newFormat_(thermoDict.lookupOrDefault("newFormat", false)),
    imbalanceTol_(thermoDict.lookupOrDefault("imbalanceTolerance", rootSmall))
{
    if (newFormat_)
    {
        Info<< "Reading CHEMKIN thermo data in new file format" << endl;
    }

    fileName chemkinFile(fileName(thermoDict.lookup("CHEMKINFile")).expand());

    fileName thermoFile = fileName::null;

    if (thermoDict.found("CHEMKINThermoFile"))
    {
        thermoFile = fileName(thermoDict.lookup("CHEMKINThermoFile")).expand();
    }

    fileName transportFile
    (
        fileName(thermoDict.lookup("CHEMKINTransportFile")).expand()
    );

    // allow relative file names
    fileName relPath = thermoDict.name().path();
    if (relPath.size())
    {
        if (!chemkinFile.isAbsolute())
        {
            chemkinFile = relPath/chemkinFile;
        }

        if (thermoFile != fileName::null && !thermoFile.isAbsolute())
        {
            thermoFile = relPath/thermoFile;
        }

        if (!transportFile.isAbsolute())
        {
            transportFile = relPath/transportFile;
        }
    }

    read(chemkinFile, thermoFile, transportFile);
}


// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Class
    Foam::chemkinIIReader

Description
    Foam::chemkinIIReader

SourceFiles
    chemkinIIReader.C
    chemkinIILexer.C

\*---------------------------------------------------------------------------*/

#ifndef chemkinIIReader_H
#define chemkinIIReader_H

#include "fileName.H"
#include "typeInfo.H"
#include "Switch.H"
#include "HashPtrTable.H"
#include "ReactionList.H"
#include "DynamicList.H"
#include "labelList.H"
#include "speciesTable.H"
#include "specieElement.H"
#include "atomicWeights.H"
#include "Table.H"

#include "specie.H"
#include "perfectGas.H"
#include "janafThermo.H"
#include "sensibleEnthalpy.H"
#include "sutherlandTransport.H"
#include "thermo.H"

#include <FlexLexer.h>

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

namespace Foam
{



typedef HashTable<List<specieElement>> speciesCompositionTable;

/*---------------------------------------------------------------------------*\
                        Class chemkinIIReader Declaration
\*---------------------------------------------------------------------------*/

class chemkinIIReader
:
    public yyFlexLexer
{

public:

    // Public data types

        enum phase
        {
            solid,
            liquid,
            gas
        };

    // Public typedefs

        typedef
            sutherlandTransport
            <
                species::thermo
                <
                    janafThermo
                    <
                        perfectGas<specie>
                    >,
                    sensibleEnthalpy
                >
            > thermoPhysics;


private:

    // Private Data

        static int yyBufSize;
        label lineNo_;

        //- Table of reaction type keywords
        HashTable<int> reactionKeywordTable_;

        //- Currently supported reaction types
        enum reactionKeyword
        {
            thirdBodyReactionType,
            unimolecularFallOffReactionType,
            chemicallyActivatedBimolecularReactionType,
            TroeReactionType,
            SRIReactionType,
            LandauTellerReactionType,
            reverseLandauTellerReactionType,
            JanevReactionType,
            powerSeriesReactionRateType,
            radiationActivatedReactionType,
            speciesTempReactionType,
            energyLossReactionType,
            plasmaMomentumTransfer,
            collisionCrossSection,
            nonEquilibriumReversibleReactionType,
            duplicateReactionType,
            speciesOrderForward,
            speciesOrderReverse,
            UnitsOfReaction,
            pressureDependentReactionType,
            end
        };

        enum reactionType
        {
            irreversible,
            reversible,
            nonEquilibriumReversible,
            unknownReactionType
        };

        static const char* reactionTypeNames[4];

        enum reactionRateType
        {
            Arrhenius,
            pressureDependentArrhenius,
            thirdBodyArrhenius,
            unimolecularFallOff,
            chemicallyActivatedBimolecular,
            LandauTeller,
            Janev,
            powerSeries,
            unknownReactionRateType
        };

        static const char* reactionRateTypeNames[9];

        enum fallOffFunctionType
        {
            Lindemann,
            Troe,
            SRI,
            unknownFallOffFunctionType
        };

        static const char* fallOffFunctionNames[4];


        void initReactionKeywordTable();


        //- List of elements
        DynamicList<word> elementNames_;

        //- Element indices
        HashTable<label> elementIndices_;

        //- Isotope molecular weights
        HashTable<scalar> isotopeAtomicWts_;

        //- List of species
        DynamicList<word> specieNames_;

        //- Specie indices
        HashTable<label> specieIndices_;

        //- Table of species
        speciesTable& speciesTable_;

        //- Specie phase
        HashTable<phase> speciePhase_;

        //- Table of the thermodynamic data given in the CHEMKIN file
        HashPtrTable<thermoPhysics> speciesThermo_;

        //- Table of species composition
        speciesCompositionTable speciesComposition_;

        //- List of the reactions
        ReactionList<thermoPhysics> reactions_;

        //- Transport properties dictionary
        dictionary transportDict_;

        //- Flag to indicate that file is in new format
        Switch newFormat_;

        //- Tolerance for element imbalance in a reaction
        scalar imbalanceTol_;


    // Private Member Functions

        //- Flex lexer to read the CHEMKIN III file
        int lex();

        inline scalar stringToScalar(const string& s)
        {
            string& str = const_cast<string&>(s);
            str.replaceAll(" ", "");
            str.replaceAll("D", "e");
            str.replaceAll("d", "e");
            return atof(str.c_str());
        }

        inline scalar stringToScalar(const char* cstr)
        {
            return stringToScalar(string(cstr));
        }

        inline void correctElementName(word& elementName)
        {
            if (elementName.size() == 2)
            {
                elementName[1] = tolower(elementName[1]);
            }
            else if (elementName[0] == 'E')
            {
                elementName = 'e';
            }
        }

        scalar molecularWeight
        (
            const List<specieElement>& specieComposition
        ) const;

        void finishElements(labelList& currentAtoms);

        void checkCoeffs
        (
            const scalarList& reactionCoeffs,
            const char* reactionRateName,
            const label nCoeffs
        ) const;

        template<class ReactionRateType>
        void addReactionType
        (
            const reactionType rType,
            DynamicList<specieCoeffs>& lhs,
            DynamicList<specieCoeffs>& rhs,
            const ReactionRateType& rr
        );

        template<template<class, class> class PressureDependencyType>
        void addPressureDependentReaction
        (
            const reactionType rType,
            const fallOffFunctionType fofType,
            DynamicList<specieCoeffs>& lhs,
            DynamicList<specieCoeffs>& rhs,
            const scalarList& thirdBodyEfficiencies,
            const scalarList& k0Coeffs,
            const scalarList& kInfCoeffs,
            const HashTable<scalarList>& reactionCoeffsTable,
            const scalar Afactor0,
            const scalar AfactorInf,
            const scalar RR
        );

        void addReaction
        (
            DynamicList<specieCoeffs>& lhs,
            DynamicList<specieCoeffs>& rhs,
            const scalarList& thirdBodyEfficiencies,
            const reactionType rType,
            const reactionRateType rrType,
            const fallOffFunctionType fofType,
            const scalarList& ArrheniusReactionCoeffs,
            HashTable<scalarList>& reactionCoeffsTable,
            const scalar RR
        );

        // Read the CHEMKIN files
        void read
        (
            const fileName& CHEMKINFileName,
            const fileName& thermoFileName,
            const fileName& transportFileName
        );

        //- Disallow default bitwise copy construction
        chemkinIIReader(const chemkinIIReader&) = delete;

        //- Disallow default bitwise assignment
        void operator=(const chemkinIIReader&) = delete;


public:

    // Constructors

        //- Construct from CHEMKIN III file name
        chemkinIIReader
        (
            speciesTable& species,
            const fileName& CHEMKINFileName,
            const fileName& transportFileName,
            const fileName& thermoFileName = fileName::null,
            const bool newFormat = false
        );

        //- Construct by getting the CHEMKIN III file name from dictionary
        chemkinIIReader(const dictionary& thermoDict, speciesTable& species);


    //- Destructor
    virtual ~chemkinIIReader()
    {}


    // Member Functions

        //- List of elements
        const wordList& elementNames() const
        {
            return elementNames_;
        }

        //- Element indices
        const HashTable<label>& elementIndices() const
        {
            return elementIndices_;
        }

        //- Isotope molecular weights
        const HashTable<scalar>& isotopeAtomicWts() const
        {
            return isotopeAtomicWts_;
        }

        //- Table of species
        const speciesTable& species() const
        {
            return speciesTable_;
        }

        //- Table of species composition
        const speciesCompositionTable& specieComposition() const
        {
            return speciesComposition_;
        }

        //- Specie phase
        const HashTable<phase>& speciePhase() const
        {
            return speciePhase_;
        }

        //- Table of the thermodynamic data given in the CHEMKIN file
        const HashPtrTable<thermoPhysics>& speciesThermo() const
        {
            return speciesThermo_;
        }

        //- List of the reactions
        const ReactionList<thermoPhysics>& reactions() const
        {
            return reactions_;
        }
};


// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

} // End namespace Foam

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

#endif

// ************************************************************************* //

/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     | Website:  https://openfoam.org
    \\  /    A nd           | Copyright (C) 2011-2020 OpenFOAM Foundation
     \\/     M anipulation  |
-------------------------------------------------------------------------------
License
    This file is part of OpenFOAM.

    OpenFOAM is free software: you can redistribute it and/or modify it
    under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
    ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
    FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
    for more details.

    You should have received a copy of the GNU General Public License
    along with OpenFOAM.  If not, see <http://www.gnu.org/licenses/>.

Application
    chemkinIIToFoam

Description
    Converts CHEMKINIII thermodynamics and reaction data files into
    OpenFOAM format.
    Updated to cover also PLOG reactions

\*---------------------------------------------------------------------------*/

#include "argList.H"
#include "chemkinIIReader.H"
#include "OFstream.H"
#include "OStringStream.H"
#include "IStringStream.H"

using namespace Foam;

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
    #include "removeCaseOptions.H"

    // Increase the precision of the output for JANAF coefficients
    Ostream::defaultPrecision(10);

    argList::validArgs.append("CHEMKIN file");
    argList::validArgs.append("CHEMKIN thermodynamics file");
    argList::validArgs.append("CHEMKIN transport file");
    argList::validArgs.append("OpenFOAM chemistry file");
    argList::validArgs.append("OpenFOAM thermodynamics file");

    argList::addBoolOption
    (
        "newFormat",
        "read Chemkin thermo file in new format"
    );

    argList args(argc, argv);

    bool newFormat = args.optionFound("newFormat");

    speciesTable species;
    chemkinIIReader cr(species, args[1], args[3], args[2], newFormat);
    const HashPtrTable<chemkinIIReader::thermoPhysics>& speciesThermo =
        cr.speciesThermo();

    dictionary thermoDict;
    thermoDict.add("species", cr.species());

    // Add the species thermo formatted entries
    {
        OStringStream os;
        speciesThermo.write(os);
        dictionary speciesThermoDict(IStringStream(os.str())());
        thermoDict.merge(speciesThermoDict);
    }

    // Temporary hack to splice the specie composition data into the thermo file
    // pending complete integration into the thermodynamics structure

    // Add elements
    forAllConstIter
    (
        HashPtrTable<chemkinIIReader::thermoPhysics>,
        speciesThermo,
        iter
    )
    {
        const word specieName(iter.key());

        dictionary elementsDict("elements");
        forAll(cr.specieComposition()[specieName], ei)
        {
            elementsDict.add
            (
                cr.specieComposition()[specieName][ei].name(),
                cr.specieComposition()[specieName][ei].nAtoms()
            );
        }

        thermoDict.subDict(specieName).add("elements", elementsDict);
    }

    thermoDict.write(OFstream(args[5])(), false);
    cr.reactions().write(OFstream(args[4])());

    Info<< "End\n" << endl;

    return 0;
}


// ************************************************************************* //