#!/usr/bin/env python
# coding: utf-8

# In[1]:


import numpy as np
from   scipy.optimize import minimize
import pandas as pd
import matplotlib.pyplot as plt


# #### Theory

# ![image.png](attachment:image.png)

# #### Preliminary Calculation

# In[2]:


PmmHg = 440.8 
TC    = 9.158144271
zis    = [0.66,0.135,0.162,0.043]


# In[3]:


def AntoinePsat(coeff,TC):
    A,B,C = coeff
    return 10**(A-B/(TC + C))


# In[4]:


coeffs = [[7.94917,1657.462,227.02],
          [7.5788,863.35,273.15],
          [6.49454,255.6784,266.55],
          [6.69147,319.0117,266.7]]


# In[5]:


""" Calculate Psats """
Psats = [AntoinePsat(coeff,TC) for coeff in coeffs]
Psats  


# In[6]:


""" Calculate Kis"""
Kis = np.divide(Psats,PmmHg)
Kis


# In[7]:


""" solve for V"""

# define the mass balance
def vapor_phase_mass_balance(V,zis):
    """
    Equation 10.17 SVA
    """
    return np.abs(1-sum([(zi*Ki)/(1+V*(Ki-1)) for zi,Ki in list(zip(zis,Kis))]))

# minimize vapor_phase_mass_balance to solve for V

F_vapor_phase_mass_balance = lambda x : vapor_phase_mass_balance(x,zis)
res = minimize(F_vapor_phase_mass_balance, x0 = 0.5 ,method='Nelder-Mead', tol=1e-6)


# evaluate L and V
if res.success:
    V,L  = res.x[0],1-res.x[0]
else:
    print('Did not converge')
    raise


# In[8]:


""" evaluate yis"""
yis = [(zi*Ki)/(1+V*(Ki-1)) for zi,Ki in list(zip(zis,Kis))]
yis


# In[9]:


""" evaluate xis"""
xis = [yi/Ki for yi,Ki in list(zip(yis,Kis))]
xis


# In[10]:


L*xis[0],1*zis[0]*0.99


#  #### T99 Calculation

# In[11]:


zis    = [0.66,0.135,0.162,0.043]


# In[12]:


def calcT98(TC,zis):
    """ Calculate Psats """
    Psats = [AntoinePsat(coeff,TC) for coeff in coeffs]

    """ Calculate Kis"""
    Kis = np.divide(Psats,PmmHg)
    
    """ solve for the vapor fraction,V"""
    F_vapor_phase_mass_balance = lambda x : vapor_phase_mass_balance(x,zis)
    res = minimize(F_vapor_phase_mass_balance, x0 = 0.5 ,method='Nelder-Mead', tol=1e-6)
    
    # evaluate L and V
    if res.success:
        V,L  = res.x[0],1-res.x[0]
    else:
        print('Did not converge')
        raise
        
    """ evaluate yis"""
    yis = [(zi*Ki)/(1+V*(Ki-1)) for zi,Ki in list(zip(zis,Kis))]
    
    """ evaluate xis"""
    xis = [yi/Ki for yi,Ki in list(zip(yis,Kis))]
        
    return [np.abs(L*xis[0] - (1*zis[0]*0.98)),yis[0]]
                            

F_calcT98 = lambda x : calcT98(x,zis)[0]


# In[13]:


T98 = minimize(F_calcT98, x0 = 10 ,method='Nelder-Mead', tol=1e-6)


# In[14]:


T98 = T98.x[0]
T98


# #### Effect of Composition to T98

# In[15]:


"""water_to_carbondioxide_mass_ratios""" 
r1s = np.linspace(0.05,2,11)
"""wetcarbondioxide_to_air_mole_ratios"""
r2s = 1-np.linspace(0.00,0.40,11) 


# In[16]:


r1 = r1s[0]
r2 = r2s[0]


# In[17]:


results = []

for r1 in r1s:
    for r2 in r2s:

        mCO2   = 1
        mwater = r1*mCO2

        nCO2   = mCO2/ 44.01
        nwater = mwater / 18.01528

        z1_pure = nwater / (nwater + nCO2)
        z2_pure = nCO2   / (nwater + nCO2)

        zs  = [    r2*z1_pure,
                   r2*z2_pure,
                   0.79*(1 - r2*z1_pure - r2*z2_pure),
                   0.21*(1 - r2*z1_pure - r2*z2_pure)]

        zijs = [np.round(z,3) for z in zs]


        """define a lambda funtion at the calculated zijs"""
        F_calcT98 = lambda x : calcT98(x,zijs)[0]
        """solve for T98"""
        T98 = minimize(F_calcT98, x0 = 50 ,method='Nelder-Mead', tol=1e-6) 

        """append the result and the calculation parameter levels"""
        results.append({'zH2O' : zijs[0],
                        'zCO2' : zijs[1],
                        'zN2'  : zijs[2],
                        'zO2'  : zijs[3],
                        'r1'   : r1,
                        'r2'   : 1-r2,
                        'yH2O' : calcT98(T98.x[0],zijs)[1],
                        'T98'  : T98.x[0]})


# In[18]:


results = pd.DataFrame(results)
results.head()


# In[19]:


df=results.pivot('r1','r2','yH2O')
df


# In[20]:


plt.figure(dpi=100)
X=df.columns.values
Y=df.index.values
Z=df.values
x,y=np.meshgrid(X, Y)
plt.contourf(x, y, Z) #the NAN will be plotted as white spaces
plt.colorbar(label = 'y$_{water}$')
plt.xlabel('mole air / (mole CO$_2$ + H$_2$O)')
plt.ylabel('mass H$_2$O / mass CO$_2$')


# In[21]:


df=results.pivot('r1','r2','T98')
df


# In[22]:


plt.figure(dpi=100)
X=df.columns.values
Y=df.index.values
Z=df.values
x,y=np.meshgrid(X, Y)
plt.contourf(x, y, Z) #the NAN will be plotted as white spaces
plt.colorbar(label = 'T$_{98}$,[°C]')
plt.xlabel('mole air / (mole CO$_2$ + H$_2$O)')
plt.ylabel('mass H$_2$O / mass CO$_2$')
plt.show()


# //CBa20220220