/**
 *  @file MolalityVPSSTP.cpp
 *   Definitions for intermediate ThermoPhase object for phases which
 *   employ molality based activity coefficient formulations
 *  (see \ref thermoprops
 * and class \link Cantera::MolalityVPSSTP MolalityVPSSTP\endlink).
 *
 * Header file for a derived class of ThermoPhase that handles
 * variable pressure standard state methods for calculating
 * thermodynamic properties that are further based upon activities
 * based on the molality scale.  These include most of the methods for
 * calculating liquid electrolyte thermodynamics.
 */
/*
 * Copyright (2005) Sandia Corporation. Under the terms of
 * Contract DE-AC04-94AL85000 with Sandia Corporation, the
 * U.S. Government retains certain rights in this software.
 */
#include "cantera/thermo/MolalityVPSSTP.h"
#include "cantera/base/stringUtils.h"

#include <iomanip>
#include <cstdio>
#include <fstream>

using namespace std;

namespace Cantera
{

MolalityVPSSTP::MolalityVPSSTP() :
    VPStandardStateTP(),
    m_indexSolvent(0),
    m_pHScalingType(PHSCALE_PITZER),
    m_indexCLM(npos),
    m_weightSolvent(18.01528),
    m_xmolSolventMIN(0.01),
    m_Mnaught(18.01528E-3)
{
    /*
     * Change the default to be that charge neutrality in the
     * phase is necessary condition for the proper specification
     * of thermodynamic functions within the phase
     */
    m_chargeNeutralityNecessary = true;
}

MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP& b) :
    VPStandardStateTP(),
    m_indexSolvent(b.m_indexSolvent),
    m_pHScalingType(b.m_pHScalingType),
    m_indexCLM(b.m_indexCLM),
    m_xmolSolventMIN(b.m_xmolSolventMIN),
    m_Mnaught(b.m_Mnaught),
    m_molalities(b.m_molalities)
{
    *this = operator=(b);
}

MolalityVPSSTP& MolalityVPSSTP::
operator=(const MolalityVPSSTP& b)
{
    if (&b != this) {
        VPStandardStateTP::operator=(b);
        m_indexSolvent     = b.m_indexSolvent;
        m_pHScalingType    = b.m_pHScalingType;
        m_indexCLM         = b.m_indexCLM;
        m_weightSolvent    = b.m_weightSolvent;
        m_xmolSolventMIN   = b.m_xmolSolventMIN;
        m_Mnaught          = b.m_Mnaught;
        m_molalities       = b.m_molalities;
    }
    return *this;
}

ThermoPhase*
MolalityVPSSTP::duplMyselfAsThermoPhase() const
{
    return new MolalityVPSSTP(*this);
}

/*
 *  -------------- Utilities -------------------------------
 */

int MolalityVPSSTP::eosType() const
{
    return 0;
}

void MolalityVPSSTP::setpHScale(const int pHscaleType)
{
    m_pHScalingType = pHscaleType;
    if (pHscaleType !=  PHSCALE_PITZER && pHscaleType !=  PHSCALE_NBS) {
        throw CanteraError("MolalityVPSSTP::setpHScale",
                           "Unknown scale type: " + int2str(pHscaleType));
    }
}

int MolalityVPSSTP::pHScale() const
{
    return m_pHScalingType;
}

void MolalityVPSSTP::setSolvent(size_t k)
{
    if (k >= m_kk) {
        throw CanteraError("MolalityVPSSTP::setSolute ",
                           "bad value");
    }
    m_indexSolvent = k;
    AssertThrowMsg(m_indexSolvent==0, "MolalityVPSSTP::setSolvent",
                   "Molality-based methods limit solvent id to being 0");
    m_weightSolvent = molecularWeight(k);
    m_Mnaught = m_weightSolvent / 1000.;
}

size_t MolalityVPSSTP::solventIndex() const
{
    return m_indexSolvent;
}

void  MolalityVPSSTP::
setMoleFSolventMin(doublereal xmolSolventMIN)
{
    if (xmolSolventMIN <= 0.0) {
        throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
    } else if (xmolSolventMIN > 0.9) {
        throw CanteraError("MolalityVPSSTP::setSolute ", "trouble");
    }
    m_xmolSolventMIN = xmolSolventMIN;
}

doublereal MolalityVPSSTP::moleFSolventMin() const
{
    return m_xmolSolventMIN;
}

void MolalityVPSSTP::calcMolalities() const
{
    getMoleFractions(DATA_PTR(m_molalities));
    double xmolSolvent = m_molalities[m_indexSolvent];
    if (xmolSolvent < m_xmolSolventMIN) {
        xmolSolvent = m_xmolSolventMIN;
    }
    double denomInv = 1.0/ (m_Mnaught * xmolSolvent);
    for (size_t k = 0; k < m_kk; k++) {
        m_molalities[k] *= denomInv;
    }
}

void MolalityVPSSTP::getMolalities(doublereal* const molal) const
{
    calcMolalities();
    for (size_t k = 0; k < m_kk; k++) {
        molal[k] = m_molalities[k];
    }
}

void MolalityVPSSTP::setMolalities(const doublereal* const molal)
{
    double Lsum = 1.0 / m_Mnaught;
    for (size_t k = 1; k < m_kk; k++) {
        m_molalities[k] = molal[k];
        Lsum += molal[k];
    }
    double tmp = 1.0 / Lsum;
    m_molalities[m_indexSolvent] = tmp / m_Mnaught;
    double sum = m_molalities[m_indexSolvent];
    for (size_t k = 1; k < m_kk; k++) {
        m_molalities[k] = tmp * molal[k];
        sum += m_molalities[k];
    }
    if (sum != 1.0) {
        tmp = 1.0 / sum;
        for (size_t k = 0; k < m_kk; k++) {
            m_molalities[k] *= tmp;
        }
    }
    setMoleFractions(DATA_PTR(m_molalities));
    /*
     * Essentially we don't trust the input: We calculate
     * the molalities from the mole fractions that we
     * just obtained.
     */
    calcMolalities();
}

void MolalityVPSSTP::setMolalitiesByName(compositionMap& mMap)
{
    /*
     *  HKM -> Might need to be more complicated here, setting
     *         neutrals so that the existing mole fractions are
     *         preserved.
     */
    size_t kk = nSpecies();
    doublereal x;
    /*
     * Get a vector of mole fractions
     */
    vector_fp mf(kk, 0.0);
    getMoleFractions(DATA_PTR(mf));
    double xmolS = mf[m_indexSolvent];
    double xmolSmin = std::max(xmolS, m_xmolSolventMIN);
    compositionMap::iterator p;
    for (size_t k = 0; k < kk; k++) {
        p = mMap.find(speciesName(k));
        if (p != mMap.end()) {
            x = mMap[speciesName(k)];
            if (x > 0.0) {
                mf[k] = x * m_Mnaught * xmolSmin;
            }
        }
    }
    /*
     * check charge neutrality
     */
    size_t largePos = npos;
    double cPos = 0.0;
    size_t largeNeg = npos;
    double cNeg = 0.0;
    double sum = 0.0;
    for (size_t k = 0; k < kk; k++) {
        double ch = charge(k);
        if (mf[k] > 0.0) {
            if (ch > 0.0) {
                if (ch * mf[k] > cPos) {
                    largePos = k;
                    cPos = ch * mf[k];
                }
            }
            if (ch < 0.0) {
                if (fabs(ch) * mf[k] > cNeg) {
                    largeNeg = k;
                    cNeg = fabs(ch) * mf[k];
                }
            }
        }
        sum += mf[k] * ch;
    }
    if (sum != 0.0) {
        if (sum > 0.0) {
            if (cPos > sum) {
                mf[largePos] -= sum / charge(largePos);
            } else {
                throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
                                   "unbalanced charges");
            }
        } else {
            if (cNeg > (-sum)) {
                mf[largeNeg] -= (-sum) / fabs(charge(largeNeg));
            } else {
                throw CanteraError("MolalityVPSSTP:setMolalitiesbyName",
                                   "unbalanced charges");
            }
        }

    }
    sum = 0.0;
    for (size_t k = 0; k < kk; k++) {
        sum += mf[k];
    }
    sum = 1.0/sum;
    for (size_t k = 0; k < kk; k++) {
        mf[k] *= sum;
    }
    setMoleFractions(DATA_PTR(mf));
    /*
     * After we formally set the mole fractions, we
     * calculate the molalities again and store it in
     * this object.
     */
    calcMolalities();
}

void MolalityVPSSTP::setMolalitiesByName(const std::string& x)
{
    compositionMap xx = parseCompString(x, speciesNames());
    setMolalitiesByName(xx);
}

/*
 * - Activities, Standard States, Activity Concentrations -----------
 */

int MolalityVPSSTP::activityConvention() const
{
    return cAC_CONVENTION_MOLALITY;
}

void MolalityVPSSTP::getActivityConcentrations(doublereal* c) const
{
    err("getActivityConcentrations");
}

doublereal MolalityVPSSTP::standardConcentration(size_t k) const
{
    err("standardConcentration");
    return -1.0;
}

doublereal MolalityVPSSTP::logStandardConc(size_t k) const
{
    err("logStandardConc");
    return -1.0;
}

void MolalityVPSSTP::getActivities(doublereal* ac) const
{
    err("getActivities");
}

void MolalityVPSSTP::getActivityCoefficients(doublereal* ac) const
{
    getMolalityActivityCoefficients(ac);
    AssertThrow(m_indexSolvent==0, "MolalityVPSSTP::getActivityCoefficients");
    double xmolSolvent = moleFraction(m_indexSolvent);
    if (xmolSolvent < m_xmolSolventMIN) {
        xmolSolvent = m_xmolSolventMIN;
    }
    for (size_t k = 1; k < m_kk; k++) {
        ac[k] /= xmolSolvent;
    }
}

void MolalityVPSSTP::getMolalityActivityCoefficients(doublereal* acMolality) const
{
    getUnscaledMolalityActivityCoefficients(acMolality);
    applyphScale(acMolality);
}

doublereal MolalityVPSSTP::osmoticCoefficient() const
{
    /*
     * First, we calculate the activities all over again
     */
    vector_fp act(m_kk);
    getActivities(DATA_PTR(act));
    /*
     * Then, we calculate the sum of the solvent molalities
     */
    double sum = 0;
    for (size_t k = 1; k < m_kk; k++) {
        sum += std::max(m_molalities[k], 0.0);
    }
    double oc = 1.0;
    double lac = log(act[m_indexSolvent]);
    if (sum > 1.0E-200) {
        oc = - lac / (m_Mnaught * sum);
    }
    return oc;
}

void MolalityVPSSTP::getElectrochemPotentials(doublereal* mu) const
{
    getChemPotentials(mu);
    double ve = Faraday * electricPotential();
    for (size_t k = 0; k < m_kk; k++) {
        mu[k] += ve*charge(k);
    }
}

doublereal MolalityVPSSTP::err(const std::string& msg) const
{
    throw CanteraError("MolalityVPSSTP","Base class method "
                       +msg+" called. Equation of state type: "+int2str(eosType()));
    return 0;
}

void MolalityVPSSTP::getUnitsStandardConc(double* uA, int k, int sizeUA) const
{
    for (int i = 0; i < sizeUA; i++) {
        if (i == 0) {
            uA[0] = 1.0;
        }
        if (i == 1) {
            uA[1] = -int(nDim());
        }
        if (i == 2) {
            uA[2] = 0.0;
        }
        if (i == 3) {
            uA[3] = 0.0;
        }
        if (i == 4) {
            uA[4] = 0.0;
        }
        if (i == 5) {
            uA[5] = 0.0;
        }
    }
}

void MolalityVPSSTP::setToEquilState(const doublereal* lambda_RT)
{
    updateStandardStateThermo();
    err("setToEquilState");
}

void MolalityVPSSTP::setStateFromXML(const XML_Node& state)
{
    VPStandardStateTP::setStateFromXML(state);
    string comp = ctml::getChildValue(state,"soluteMolalities");
    if (comp != "") {
        setMolalitiesByName(comp);
    }
    if (state.hasChild("pressure")) {
        double p = ctml::getFloat(state, "pressure", "pressure");
        setPressure(p);
    }
}

void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p,
                                  const doublereal* const molalities)
{
    setMolalities(molalities);
    setState_TP(t, p);
}

void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, compositionMap& m)
{
    setMolalitiesByName(m);
    setState_TP(t, p);
}

void MolalityVPSSTP::setState_TPM(doublereal t, doublereal p, const std::string& m)
{
    setMolalitiesByName(m);
    setState_TP(t, p);
}


void MolalityVPSSTP::initThermo()
{
    initLengths();
    VPStandardStateTP::initThermo();

    /*
     * The solvent defaults to species 0
     */
    setSolvent(0);
    /*
     * Find the Cl- species
     */
    m_indexCLM = findCLMIndex();
}

void MolalityVPSSTP::getUnscaledMolalityActivityCoefficients(doublereal* acMolality) const
{
    err("getUnscaledMolalityActivityCoefficients");
}

void MolalityVPSSTP::applyphScale(doublereal* acMolality) const
{
    err("applyphScale");
}

size_t MolalityVPSSTP::findCLMIndex() const
{
    size_t indexCLM = npos;
    size_t eCl = npos;
    size_t eE = npos;
    size_t ne = nElements();
    string sn;
    for (size_t e = 0; e < ne; e++) {
        sn = elementName(e);
        if (sn == "Cl" || sn == "CL") {
            eCl = e;
            break;
        }
    }
    // We have failed if we can't find the Cl element index
    if (eCl == npos) {
        return npos;
    }
    for (size_t e = 0; e < ne; e++) {
        sn = elementName(e);
        if (sn == "E" || sn == "e") {
            eE = e;
            break;
        }
    }
    // We have failed if we can't find the E element index
    if (eE == npos) {
        return npos;
    }
    for (size_t k = 1; k < m_kk; k++) {
        doublereal nCl = nAtoms(k, eCl);
        if (nCl != 1.0) {
            continue;
        }
        doublereal nE = nAtoms(k, eE);
        if (nE != 1.0) {
            continue;
        }
        for (size_t e = 0; e < ne; e++) {
            if (e != eE && e != eCl) {
                doublereal nA = nAtoms(k, e);
                if (nA != 0.0) {
                    continue;
                }
            }
        }
        sn = speciesName(k);
        if (sn != "Cl-" && sn != "CL-") {
            continue;
        }

        indexCLM = k;
        break;
    }
    return indexCLM;
}

//   Initialize lengths of local variables after all species have
//   been identified.
void  MolalityVPSSTP::initLengths()
{
    m_kk = nSpecies();
    m_molalities.resize(m_kk);
}

void MolalityVPSSTP::initThermoXML(XML_Node& phaseNode, const std::string& id_)
{

    initLengths();
    /*
     * The solvent defaults to species 0
     */
    setSolvent(0);

    VPStandardStateTP::initThermoXML(phaseNode, id_);
}

/**
  * Format a summary of the mixture state for output.
  */
std::string MolalityVPSSTP::report(bool show_thermo) const
{
    char p[800];
    string s = "";
    try {
        if (name() != "") {
            sprintf(p, " \n  %s:\n", name().c_str());
            s += p;
        }
        sprintf(p, " \n       temperature    %12.6g  K\n", temperature());
        s += p;
        sprintf(p, "          pressure    %12.6g  Pa\n", pressure());
        s += p;
        sprintf(p, "           density    %12.6g  kg/m^3\n", density());
        s += p;
        sprintf(p, "  mean mol. weight    %12.6g  amu\n", meanMolecularWeight());
        s += p;

        doublereal phi = electricPotential();
        sprintf(p, "         potential    %12.6g  V\n", phi);
        s += p;

        size_t kk = nSpecies();
        vector_fp x(kk);
        vector_fp molal(kk);
        vector_fp mu(kk);
        vector_fp muss(kk);
        vector_fp acMolal(kk);
        vector_fp actMolal(kk);
        getMoleFractions(&x[0]);
        getMolalities(&molal[0]);
        getChemPotentials(&mu[0]);
        getStandardChemPotentials(&muss[0]);
        getMolalityActivityCoefficients(&acMolal[0]);
        getActivities(&actMolal[0]);

        size_t iHp = speciesIndex("H+");
        if (iHp != npos) {
            double pH = -log(actMolal[iHp]) / log(10.0);
            sprintf(p, "                pH    %12.4g  \n", pH);
            s += p;
        }

        if (show_thermo) {
            sprintf(p, " \n");
            s += p;
            sprintf(p, "                          1 kg            1 kmol\n");
            s += p;
            sprintf(p, "                       -----------      ------------\n");
            s += p;
            sprintf(p, "          enthalpy    %12.6g     %12.4g     J\n",
                    enthalpy_mass(), enthalpy_mole());
            s += p;
            sprintf(p, "   internal energy    %12.6g     %12.4g     J\n",
                    intEnergy_mass(), intEnergy_mole());
            s += p;
            sprintf(p, "           entropy    %12.6g     %12.4g     J/K\n",
                    entropy_mass(), entropy_mole());
            s += p;
            sprintf(p, "    Gibbs function    %12.6g     %12.4g     J\n",
                    gibbs_mass(), gibbs_mole());
            s += p;
            sprintf(p, " heat capacity c_p    %12.6g     %12.4g     J/K\n",
                    cp_mass(), cp_mole());
            s += p;
            try {
                sprintf(p, " heat capacity c_v    %12.6g     %12.4g     J/K\n",
                        cv_mass(), cv_mole());
                s += p;
            } catch (CanteraError& e) {
                e.save();
                sprintf(p, " heat capacity c_v    <not implemented>       \n");
                s += p;
            }
        }

        sprintf(p, " \n");
        s += p;
        if (show_thermo) {
            sprintf(p, "                           X        "
                    "   Molalities         Chem.Pot.    ChemPotSS    ActCoeffMolal\n");
            s += p;
            sprintf(p, "                                    "
                    "                      (J/kmol)      (J/kmol)                 \n");
            s += p;
            sprintf(p, "                     -------------  "
                    "  ------------     ------------  ------------    ------------\n");
            s += p;
            for (size_t k = 0; k < kk; k++) {
                if (x[k] > SmallNumber) {
                    sprintf(p, "%18s  %12.6g     %12.6g     %12.6g   %12.6g   %12.6g\n",
                            speciesName(k).c_str(), x[k], molal[k], mu[k], muss[k], acMolal[k]);
                } else {
                    sprintf(p, "%18s  %12.6g     %12.6g          N/A      %12.6g   %12.6g \n",
                            speciesName(k).c_str(), x[k], molal[k], muss[k], acMolal[k]);
                }
                s += p;
            }
        } else {
            sprintf(p, "                           X"
                    "Molalities\n");
            s += p;
            sprintf(p, "                     -------------"
                    "     ------------\n");
            s += p;
            for (size_t k = 0; k < kk; k++) {
                sprintf(p, "%18s   %12.6g     %12.6g\n",
                        speciesName(k).c_str(), x[k], molal[k]);
                s += p;
            }
        }
    } catch (CanteraError& err) {
        err.save();
    }
    return s;
}

void MolalityVPSSTP::getCsvReportData(std::vector<std::string>& names,
                                      std::vector<vector_fp>& data) const
{
    names.clear();
    data.assign(10, vector_fp(nSpecies()));

    names.push_back("X");
    getMoleFractions(&data[0][0]);

    names.push_back("Molal");
    getMolalities(&data[1][0]);

    names.push_back("Chem. Pot. (J/kmol)");
    getChemPotentials(&data[2][0]);

    names.push_back("Chem. Pot. SS (J/kmol)");
    getStandardChemPotentials(&data[3][0]);

    names.push_back("Molal Act. Coeff.");
    getMolalityActivityCoefficients(&data[4][0]);

    names.push_back("Molal Activity");
    getActivities(&data[5][0]);

    names.push_back("Part. Mol Enthalpy (J/kmol)");
    getPartialMolarEnthalpies(&data[5][0]);

    names.push_back("Part. Mol. Entropy (J/K/kmol)");
    getPartialMolarEntropies(&data[6][0]);

    names.push_back("Part. Mol. Energy (J/kmol)");
    getPartialMolarIntEnergies(&data[7][0]);

    names.push_back("Part. Mol. Cp (J/K/kmol");
    getPartialMolarCp(&data[8][0]);

    names.push_back("Part. Mol. Cv (J/K/kmol)");
    getPartialMolarVolumes(&data[9][0]);
}

}