Method init

src/kinetics/ReactionPath.cpp:601–718

}

int ReactionPathBuilder::init(ostream& logfile, Kinetics& kin)
{
    m_transfer.clear();
    m_elementSymbols.clear();
    findElements(kin);
    m_ns = kin.nTotalSpecies();
    m_nr = kin.nReactions();

    // all reactants / products, even ones appearing on both sides of the
    // reaction
    vector<vector<size_t>> allProducts(m_nr);
    vector<vector<size_t>> allReactants(m_nr);
    for (size_t i = 0; i < m_nr; i++) {
        for (size_t k = 0; k < m_ns; k++) {
            for (int n = 0; n < kin.reactantStoichCoeff(k, i); n++) {
                allReactants[i].push_back(k);
            }
            for (int n = 0; n < kin.productStoichCoeff(k, i); n++) {
                allProducts[i].push_back(k);
            }
        }
    }

    // m_reac and m_prod exclude indices for species that appear on
    // both sides of the reaction, so that the diagram contains no loops.
    m_reac.resize(m_nr);
    m_prod.resize(m_nr);
    m_ropf.resize(m_nr);
    m_ropr.resize(m_nr);
    m_determinate.resize(m_nr);
    m_x.resize(m_ns); // not currently used ?
    m_elatoms.resize(m_nel, m_nr);

    for (size_t i = 0; i < m_nr; i++) {
        // construct the lists of reactant and product indices, not including
        // molecules that appear on both sides.
        m_reac[i].clear();
        m_prod[i].clear();
        map<size_t, int> net;
        size_t nr = allReactants[i].size();
        size_t np = allProducts[i].size();
        for (size_t ir = 0; ir < nr; ir++) {
            net[allReactants[i][ir]]--;
        }
        for (size_t ip = 0; ip < np; ip++) {
            net[allProducts[i][ip]]++;
        }

        for (size_t k = 0; k < m_ns; k++) {
            if (net[k] < 0) {
                size_t nmol = -net[k];
                for (size_t jr = 0; jr < nmol; jr++) {
                    m_reac[i].push_back(k);
                }
            } else if (net[k] > 0) {
                size_t nmol = net[k];
                for (size_t jp = 0; jp < nmol; jp++) {
                    m_prod[i].push_back(k);