/*---------------------------------------------------------------------------*\
  =========                 |
  \\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
   \\    /   O peration     |
    \\  /    A nd           | Copyright (C) 2011-2016 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
    chemFoam

Description
    Solver for chemistry problems, designed for use on single cell cases to
    provide comparison against other chemistry solvers, that uses a single cell
    mesh, and fields created from the initial conditions.

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

#include <cantera/transport.h>
#include <cantera/IdealGasMix.h>

#include <fstream>

#include <cmath>

#include "cxxopts/cxxopts.hpp"

#include "ofFormats.h"

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

std::string stringNASACoefs (doublereal *coef)
{
    std::string arr;

    for (int j = 0; j < 7; j++)
    {
        arr += fmt::format(" {:15.10E} ", coef[j]);
    }

    return arr;
}

void calculateSutherland (Cantera::IdealGasMix *gas_, doublereal As[], doublereal C[])
{
    Cantera::Transport *tr_ = Cantera::newTransportMgr("Mix", gas_);

    doublereal mu0[gas_->nSpecies()];
    doublereal muRatio[gas_->nSpecies()];

    const doublereal T0 = gas_->minTemp();
    gas_->setState_TP(T0, Cantera::OneAtm);
    tr_->getSpeciesViscosities(mu0);

    size_t nInterval = 100;
    doublereal T[nInterval];
    doublereal TRatioPower[nInterval];
    doublereal dTemp = (gas_->maxTemp() - gas_->minTemp()) / doublereal(nInterval);

    for (size_t i = 0; i < nInterval; i++)
    {
        T[i] = gas_->minTemp() + (i+1) * dTemp;
    }

    for (size_t i = 0; i < nInterval; i++)
    {
        TRatioPower[i] = pow(T[i]/T0, 3./2.);
    }

    for (size_t j = 0; j < gas_->nSpecies(); j++)
    {
        C[j] = 0.0;
    }

    for (size_t i = 0; i < nInterval; i++)
    {
        doublereal Ti = T[i];

        gas_->setTemperature(Ti);

        tr_->getSpeciesViscosities(muRatio);
        for (size_t j = 0; j < gas_->nSpecies(); j++)
        {
            muRatio[j] /= mu0[j];
        }

        for (size_t j = 0; j < gas_->nSpecies(); j++)
        {
            C[j] += (muRatio[j] * Ti - TRatioPower[i] * T0) / (TRatioPower[i] - muRatio[j]);
        }
    }

    for (size_t j = 0; j < gas_->nSpecies(); j++)
    {
        C[j] /= doublereal(nInterval);
    }


    for (size_t j = 0; j < gas_->nSpecies(); j++)
    {
        As[j] = mu0[j] * (T0 + C[j]) / pow(T0, 3./2.);
    }

    delete tr_;

    return ;
}

int main(int argc, char *argv[])
{

    cxxopts::Options options("ct2foam", "Utility Converts Cantera formatted data in OpenFOAM format");

    std::string scti;
    std::string sgas("");
    std::string sthermo;
    std::string srxn;

    std::ofstream fthermo;
    std::ofstream frxn;


    options.add_options()
        ("name", "Name of ideal gas mix in CTI file", cxxopts::value(sgas))
        ("h,help", "print help message")
        ;

    options.add_options("internal")
        ("cti", "Input CTI File name", cxxopts::value(scti))
        ("reaction", "Output OpenFOAM Reaction File name", cxxopts::value(srxn)->default_value("reactions"))
        ("thermo", "Output OpenFOAM Thermo File name", cxxopts::value(sthermo)->default_value("thermos"))
        ("rest", "Input CTI File name", cxxopts::value<std::vector<std::string>>())
        ;

    options.parse_positional({"cti", "reaction", "thermo", "rest"});

    options.positional_help("CTI [REACTION [THERMO]]");

    cxxopts::ParseResult result = options.parse(argc, argv);

    if (result.count("help"))
    {
        std::cout << options.help();

        return 0;
    }

    if (!result["cti"].count())
    {
        std::cerr << "CTI file name must be provided" << std::endl << std::endl;

        std::cout << options.help();

        return -1;
    }

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

    Cantera::IdealGasMix gas_(scti, sgas);

    fthermo.open(sthermo);

    frxn.open(srxn);

    // Thermo Part =============================================================

    // Calculate Sutherland Coefficients
    doublereal As[gas_.nSpecies()] ;
    doublereal C[gas_.nSpecies()] ;

    calculateSutherland (&gas_, As, C);

    for (size_t n = 0; n < gas_.nSpecies(); n++)
    {
        std::cout << As[n] << std::endl;
    }

    std::cout << std::endl;

    for (size_t n = 0; n < gas_.nSpecies(); n++)
    {
        std::cout << C[n] << std::endl;
    }

    // these constants define the location of coefficient "a6" in the
    // cofficient array c. The c array contains Tmid in the first
    // location, followed by the 7 low-temperature coefficients, then
    // the seven high-temperature ones.
    for (size_t n = 0; n < gas_.nSpecies(); n++)
    {
        int type;
        doublereal c[15];
        doublereal minTemp, maxTemp, refPressure;

        // get the NASA coefficients in array c
        gas_.speciesThermo().reportParams(n, type, c, minTemp, maxTemp, refPressure);

        fthermo<<(
                fmt::format(
                    thermoFormat,
                    gas_.speciesName(n),
                    gas_.molecularWeight(n),
                    gas_.charge(n),
                    minTemp,
                    maxTemp,
                    c[0],
                    stringNASACoefs(&c[1]),
                    stringNASACoefs(&c[8]),
                    As[n],
                    C[n]
                    )
                );
    }

    fthermo.close();

    // End of Thermo Part ======================================================


    // Reaction Part

    // Write species name list
    frxn<<("species\n");
    frxn<<fmt::format("{}\n", gas_.nSpecies());
    frxn<<("(\n");
    for (size_t n = 0; n < gas_.nSpecies(); n++)
    {
        frxn<<fmt::format("{}\n", gas_.speciesName(n));
    }
    frxn<<(")\n");
    frxn<<(";\n\n");

    // Write reaction list
    frxn<<( "reactions\n{\n"); // begin "reactions" dictionary

    for (size_t k = 0; k < gas_.nReactions(); k++)
    {
        std::shared_ptr<Cantera::Reaction> r(gas_.reaction(k));

        std::string irn( r->reversible ? "reversible" : "irreversible" );

        std::string rate("Arrhenius");

        std::string rxn("Reaction");

        std::string ffn;

        double c[5] = {0};

        // reaction name
        frxn<<( fmt::format( "    un-named-reaction-{}\n", k)); // begin a reaction dictionary
        frxn<<( "    {\n" );

        switch (gas_.reactionType(k))
        {
            case Cantera::ELEMENTARY_RXN:

                frxn<<( fmt::format( rxnTypeFormat, irn + rate + rxn));

                frxn<<( fmt::format( rxnEqnFormat,
                            r->Reaction::reactantString() + " = " + r->Reaction::productString()));

                frxn<<(
                        fmt::format(
                            arrheniusFormat,
                            std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor(),
                            std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent(),
                            std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R()
                            )
                        );
                break;

            case Cantera::THREE_BODY_RXN:

                rate = "thirdBody" + rate;

                frxn<<( fmt::format( rxnTypeFormat, irn + rate + rxn));

                frxn<<( fmt::format( rxnEqnFormat,
                            r->Reaction::reactantString() + " = " + r->Reaction::productString()));

                frxn<<(
                        fmt::format(
                            arrheniusFormat,
                            std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.preExponentialFactor(),
                            std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.temperatureExponent(),
                            std::dynamic_pointer_cast<Cantera::ElementaryReaction>(r)->rate.activationEnergy_R()
                            )
                        );

                frxn<<("        coeffs          \n");

                frxn<<(fmt::format("{}\n(\n", gas_.nSpecies()));

                for (size_t l = 0; l < gas_.nSpecies(); l++)
                {
                    frxn<<(
                            fmt::format(
                                "({} {})\n",
                                gas_.speciesName(l),
                                std::dynamic_pointer_cast<Cantera::ThreeBodyReaction>(r)->third_body.efficiency(gas_.speciesName(l))
                                )
                            );
                }

                frxn<<(")\n;\n");

                break;

            case Cantera::FALLOFF_RXN:
            case Cantera::CHEMACT_RXN:

                rxn = (r->reaction_type == Cantera::FALLOFF_RXN ? "FallOff" : "ChemicallyActivated") + rxn;

                switch (std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getType())
                {
                    case Cantera::SIMPLE_FALLOFF:
                        rxn = "Lindemann" + rxn;
                        ffn = lindemannFormat;
                        break;

                    case Cantera::TROE_FALLOFF:
                        rxn = "Troe" + rxn;
                        std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getParameters(c);
                        ffn = fmt::format(troeFormat, c[0], c[1], c[2], c[3]);
                        break;

                    case Cantera::SRI_FALLOFF:
                        rxn = "SRI" + rxn;
                        std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->falloff->getParameters(c);
                        ffn = fmt::format(sriFormat, c[0], c[1], c[2], c[3], c[4]);
                        break;
                }

                frxn<<( fmt::format( rxnTypeFormat, irn + rate + rxn));

                frxn<<( fmt::format( rxnEqnFormat,
                            r->Reaction::reactantString() + " = " + r->Reaction::productString()));

                frxn<<(
                        fmt::format(
                            falloffFormat,
                            std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.preExponentialFactor(),
                            std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.temperatureExponent(),
                            std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->low_rate.activationEnergy_R(),
                            std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.preExponentialFactor(),
                            std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.temperatureExponent(),
                            std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->high_rate.activationEnergy_R()
                            )
                        );

                frxn<<( ffn );

                frxn<<( std::string() +
                        "        thirdBodyEfficiencies\n" +
                        "        {\n" +
                        "            coeffs          \n"
                        );

                frxn<<(fmt::format("{}\n(\n", gas_.nSpecies()));

                for (size_t l = 0; l < gas_.nSpecies(); l++)
                {
                    frxn<<(
                            fmt::format(
                                "({} {})\n",
                                gas_.speciesName(l),
                                std::dynamic_pointer_cast<Cantera::FalloffReaction>(r)->third_body.efficiency(gas_.speciesName(l))
                                )
                            );
                }

                frxn<<(")\n;\n");
                frxn<<("        }\n");

                break;

            case Cantera::PLOG_RXN:
                break;

            case Cantera::CHEBYSHEV_RXN:
                break;

            default:
                break;
        }

        frxn<<( "    }\n" ); // end of a reaction dictionary
    }

    frxn<<( "}\n"); // end of "reactions" dictionary

    frxn.close();

    return 0;
}


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