'''
Author: Vegard G. Jervell
Date: December 2020
Purpose: Implementation of L. J. T. M. Kempers' 1989 model for prediction of Soret-coefficients in multicomponent
        liquid mixtures, derived in "A thermodynamic theory of the Soret effect in a multicomponent liquid",
        Journal of Chemical Physics, 1989.
        doi : http://dx.doi.org/10.1063/1.456321
Requires: numpy, scipy, pandas, ThermoPack
Note: pandas is used to get ideal-gas state enthalpies from the file 'ideal_gas_enthalpies.xlsx'. It should
    be possible to aquire these from ThermoPack, thereby removing the need for a separate data-file of
    standard-state enthalpies.
'''

from pycThermopack.pyctp import cubic, cpa, lee_kesler
from models.kempers import Kempers
import numpy as np
from scipy.optimize import root
import pandas as pd
import warnings
import os

class Kempers89(Kempers):
    def __init__(self, comps, eos, x=[0.5, 0.5], temp=300, pres=1e5, phase=1):
        '''
        :param comps (str): Comma separated list of components, as per ThermoPack convention
        :param x (array-like): list of component mole fractions
        :param eos (ThermoPack): Initialized equation of state, with same components as Kempers
        :param temp (float): Temperature [K]
        :param pres (float): Pressure [Pa]
        :param phase (int): Phase identifier to be used by ThermoPack
        '''
        super().__init__(comps, eos, x=x, temp=temp, pres=pres, phase=phase)

        if self.eos.guess_phase(self.temp, self.pres, self.x) != 1:
            warnings.warn('This model may work poorly for phases other than liquid,'
                          ' ThermoPack predicted phase is '+
                          eos.get_phase_type(eos.guess_phase(self.temp, self.pres, self.x)))

    def get_standard_enthalpies(self):
        '''
        Get ideal-gas standard state enthalpies from 'ideal_gas_enthalpies.xlsx'
        :return: 1darray: standard state enthalpies of components
        '''
        file_path = os.path.dirname(__file__)
        data = pd.read_excel(file_path+'/ideal_gas_enthalpies.xlsx')
        return np.array([data.loc[data['Comp'] == comp].values.flatten()[1] for comp in self.comps.split(',')])

    def get_soret_cov(self, kin=False):
        '''
        Get soret coefficients at current settings
        :return: 1darray: soret coefficients of components
        '''

        if kin is True:
            raise ValueError('Kempers89 has no kinetic contribution!')

        self.eos.set_tmin(1)
        V, v = self.eos.specific_volume(self.temp, self.pres, self.x, self.phase, dvdn=True)
        H, h = self.eos.enthalpy(self.temp, self.pres, self.x, self.phase, dhdn=True)
        H0, h0 = self.eos.enthalpy(self.temp, 1e-5, self.x, 2, dhdn=True)
        H_molar_ideal = self.get_standard_enthalpies()

        h = h - h0 + H_molar_ideal
        dmudx = self.dmudx_TP()

        N = len(self.x)

        def eq_set(alpha):
            dxdT = - (alpha * self.x * (1 - self.x))
            eqs = np.zeros(N)
            for i in range(0, N - 1):
                eqs[i] = -((h[i] / v[i]) - (h[N - 1] / v[N - 1])) \
                         + sum([((dmudx[i, j] / v[i]) - (dmudx[N - 1, j] / v[N - 1])) * dxdT[j] for j in range(0, N - 1)])

            eqs[N - 1] = sum(alpha)

            return eqs

        alpha = root(eq_set, [1 for i in range(N)]).x
        soret = alpha / self.temp

        return soret

    def get_soret_com(self, kin=False):
        '''
        Get soret coefficients at current settings
        :return: 1darray: soret coefficients of components
        '''

        if kin is True:
            raise ValueError('Kempers89 has no kinetic contribution!')

        M = self.mole_weights = np.array([self.eos.compmoleweight(self.eos.getcompindex(comp))
                                          for comp in self.comps.split(',')])
        _, h = self.eos.enthalpy(self.temp, self.pres, self.x, self.phase, dhdn=True)
        _, h0 = self.eos.enthalpy(self.temp, 1e-5, self.x, 2, dhdn=True)
        H_molar_ideal = self.get_standard_enthalpies()

        h = h - h0 + H_molar_ideal
        dmudx = self.dmudx_TP()

        N = len(self.x)

        def eq_set(alpha):
            dxdT = - (alpha * self.x * (1 - self.x))
            eqs = np.zeros(N)
            for i in range(0, N - 1):
                eqs[i] = -((h[i] / M[i]) - (h[N - 1] / M[N - 1])) \
                         + sum([((dmudx[i, j] / M[i]) - (dmudx[N - 1, j] / M[N - 1])) * dxdT[j] for j in range(0, N - 1)])

            eqs[N - 1] = sum(alpha)

            return eqs

        alpha = root(eq_set, [1 for i in range(N)]).x
        soret = alpha / self.temp

        return soret