import numpy as np
import sympy as sp
from sympy.interactive import printing
printing.init_printing(use_latex = True)
Definition of symbolic variables
g_fo_o,g_fa_o,y_fo,y_fa,T,R,W = sp.symbols('g_fo^o,g_fa^o,y_fo,y_fa,T,R,W')
Definition of gibbs energy of the olivine solid solution
g = y_fo*g_fo_o + (1-y_fo)*g_fa_o + R*T*(2*y_fo*sp.log(y_fo) + 2*(1-y_fo)*sp.log(1-y_fo)) + W*y_fo*(1-y_fo)
display(g)
Definition of the gibss energy of forsterite in the olivine solid solution
g_fo_ss = g + (1-y_fo)*sp.diff(g,y_fo)
display(sp.simplify(g_fo_ss))
Activity of forsterite is defined as the gibbs energy of the endmember in the solid solution (g_fo_ss
a_fo = sp.exp((g_fo_ss-g_fo_o)/R/T)
sp.simplify(a_fo)
print(sp.simplify(a_fo))
y_fo**2*exp(W*(y_fo**2 - 2*y_fo + 1)/(R*T))
being awkwards the expression is usually expressed as
xxxxxxxxxx
RTlna_fo = R*T*sp.log(sp.exp((g_fo_ss-g_fo_o)/R/T))
xxxxxxxxxx
sp.simplify(RTlna_fo)
Just note that
xxxxxxxxxx
sp.simplify(g_fo_ss-g_fo_o)
xxxxxxxxxx
sp.collect(W*y_fo**2-2*W*y_fo+W,W)
and that
xxxxxxxxxx
sp.factor(W*y_fo**2-2*W*y_fo+W,y_fo)
so ,