import numpy as np import sympy as spfrom 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
xxxxxxxxxxRTlna_fo = R*T*sp.log(sp.exp((g_fo_ss-g_fo_o)/R/T))xxxxxxxxxxsp.simplify(RTlna_fo)Just note that
xxxxxxxxxxsp.simplify(g_fo_ss-g_fo_o)xxxxxxxxxxsp.collect(W*y_fo**2-2*W*y_fo+W,W)and that
xxxxxxxxxxsp.factor(W*y_fo**2-2*W*y_fo+W,y_fo)so ,