// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.

#ifndef CT_EXAMPLE_UTILS_H
#define CT_EXAMPLE_UTILS_H

#include "cantera/base/Array.h"
#include <fstream>

// Save the temperature, density, pressure, and mole fractions at one
// time
template<class G, class A>
void saveSoln(int i, double time, const G& gas, A& soln)
{
    soln(0,i) = time;
    soln(1,i) = gas.temperature();
    soln(2,i) = gas.density();
    soln(3,i) = gas.pressure();
    gas.getMoleFractions(&soln(4,i));
}

template<class G, class A>
void saveSoln(double time, const G& gas, A& soln)
{
    soln.resize(soln.nRows(), soln.nColumns() + 1);
    int back = soln.nColumns() - 1;
    soln(0,back) = time;
    soln(1,back) = gas.temperature();
    soln(2,back) = gas.density();
    soln(3,back) = gas.pressure();
    int nsp = gas.nSpecies();
    for (int k = 0; k < nsp; k++) {
        soln(4+k,back) = gas.moleFraction(k);
    }
}

void writeCsv(const std::string& fname, const Cantera::ThermoPhase& gas,
              const Cantera::Array2D& data)
{
    std::ofstream s(fname);
    // Write labels
    s << "time (s),Temperature (K),Density (kg/m3),Pressure (Pa),";

    for (size_t k = 0; k < gas.nSpecies(); k++) {
        s << gas.speciesName(k);
        if (k != gas.nSpecies()-1) {
            s << ",";
        }
    }
    s << std::endl;
    for (size_t i = 0; i < data.nColumns(); i++) {
        for (size_t j = 0; j < data.nRows(); j++) {
            s << data(j,i);
            if (j != data.nRows()-1) {
                s << ",";
            }
        }
        s << std::endl;
    }
}

#endif

/*!
 * @file kinetics1.cpp
 *
 * Zero-dimensional kinetics
 *
 * This example simulates autoignition of hydrogen in a constant pressure
 * reactor and saves the time history to files that can be used for plotting.
 *
 * Keywords: combustion, reactor network, ignition delay, saving output
 */

// This file is part of Cantera. See License.txt in the top-level directory or
// at https://cantera.org/license.txt for license and copyright information.

#include "cantera/zerodim.h"
#include "cantera/numerics/Integrator.h"
#include "example_utils.h"

using namespace Cantera;
using std::cout;
using std::endl;

int kinetics1(int np, void* p)
{
    cout << "Constant-pressure ignition of a "
         << "hydrogen/oxygen/nitrogen"
         " mixture \nbeginning at T = 1001 K and P = 1 atm." << endl;

    // create an ideal gas mixture that corresponds to OH submech from GRI-Mech 3.0
    auto sol = newSolution("h2o2.yaml", "ohmech", "none");
    auto gas = sol->thermo();

    // set the state
    gas->setState_TPX(1001.0, OneAtm, "H2:2.0, O2:1.0, N2:4.0");
    int nsp = gas->nSpecies();

    // create a reactor
    IdealGasConstPressureReactor r;

    // 'insert' the gas into the reactor and environment.  Note
    // that it is ok to insert the same gas object into multiple
    // reactors or reservoirs. All this means is that this object
    // will be used to evaluate thermodynamic or kinetic
    // quantities needed.
    r.insert(sol);

    double dt = 1.e-5; // interval at which output is written
    int nsteps = 100; // number of intervals

    // create a 2D array to hold the output variables,
    // and store the values for the initial state
    Array2D states(nsp+4, 1);
    saveSoln(0, 0.0, *(sol->thermo()), states);

    // create a container object to run the simulation
    // and add the reactor to it
    ReactorNet sim;
    sim.addReactor(r);

    // main loop
    clock_t t0 = clock(); // save start time
    for (int i = 1; i <= nsteps; i++) {
        double tm = i*dt;
        sim.advance(tm);
        cout << "time = " << tm << " s" << endl;
        saveSoln(tm, *(sol->thermo()), states);
    }
    clock_t t1 = clock(); // save end time


    // make a CSV output file
    writeCsv("kin1.csv", *sol->thermo(), states);

    // print final temperature and timing data
    double tmm = 1.0*(t1 - t0)/CLOCKS_PER_SEC;
    cout << " Tfinal = " << r.temperature() << endl;
    cout << " time = " << tmm << endl;
    cout << " number of residual function evaluations = "
         << sim.integrator().nEvals() << endl;
    cout << " time per evaluation = " << tmm/sim.integrator().nEvals()
         << endl << endl;
    cout << "Output files:" << endl
         << "  kin1.csv    (Excel CSV file)" << endl;

    return 0;
}


int main()
{
    try {
        int retn = kinetics1(0, 0);
        appdelete();
        return retn;
    } catch (CanteraError& err) {
        // handle exceptions thrown by Cantera
        std::cout << err.what() << std::endl;
        cout << " terminating... " << endl;
        appdelete();
        return -1;
    }
}