#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Equilibrium Fe^3+ + SCN^- = FeSCN^2+ The equilibrium constant is taken equal to 10^2.3 and the volume equal to 100mL. We suppose that the activity is equal to the concentration of each species. Informations ------------ Author : Martin Vérot from the ENS de Lyon, France, with the help of some scripts to create the buttons and sliders taken from https://github.com/araoux/python_agregation (written by P Cladé, A Raoux and F Levrier) Licence : Creative Commons CC-BY-NC-SA 4.0 WARNING this program requires the widgets.py file to work """ #Libraries import numpy as np from scipy import optimize import scipy.constants as constants import math import matplotlib.pyplot as plt import matplotlib import matplotlib.patches as patches import widgets import time def Quotient(xieq,nFe0,nSCN0,nFeSCN0,V): K = 10**2.3 return (nFeSCN0+xieq)*V-K*((nFe0-xieq)*(nSCN0-xieq)) def plot_data(nFe0,nSCN0,nFeSCN0,xi,widget,axes): V=0.1 #finding extremums for xi ximin = min(0,-nFeSCN0) ximax = min(nFe0,nSCN0) #equilibrium constant K = 10**2.3 #starting value for optimization xieq = [(ximin+ximax)/2] #x for cure xis,dxi = np.linspace(ximin,ximax,1000,retstep=True) #updating slider extremal values widget['xi'].valmin = ximin widget['xi'].valmax = ximax widget['xi'].valstep = dxi #Computing values for the given xi nFe=nFe0-xi nSCN=nSCN0-xi nFeSCN=nFeSCN0+xi if nFe>0 and nSCN >0 and nFeSCN>=0: Q = nFeSCN*V/(nFe*nSCN) else: Q=np.nan nFe=np.nan nSCN=np.nan nFeSCN=np.nan #computing values for all possible xis nsFe = nFe0 - xis nsSCN = nSCN0 - xis nsFeSCN = nFeSCN0 + xis Qs = nsFeSCN*V/(nsFe*nsSCN) boolArr = (nsFe >0) & (nsSCN>0) & (nsFeSCN >= 0) Qdef = Qs[boolArr] xidef = xis[boolArr] #Finding root, the system is given in more soluble form to avoid numerical problems sol = optimize.root(Quotient, xieq,args=(nFe0,nSCN0,nFeSCN0,V)) xieq = sol.x[0] nFeeq=nFe0-xieq nSCNeq=nSCN0-xieq nFeSCNeq=nFeSCN0+xieq lines['Q'].set_data(xidef,Qdef) lines['eq'].set_data(sol.x,[K]) lines['Qsingle'].set_data([xi],[Q]) ax1.set_xlim(min(xis),max(xis)) ax1.set_ylim(0,1000) #Updating all text values texts['Fe0'].set_text('$n^0_\\mathrm{{Fe^{{3+}}}}= {:.3e}$'.format(nFe0)) texts['SCN0'].set_text('$n^0_\\mathrm{{SCN^{{-}}}}= {:.3e}$'.format(nSCN0)) texts['FeSCN0'].set_text('$n^0_\\mathrm{{FeSCN^{{2+}}}}= {:.3e}$'.format(nFeSCN0)) texts['Fe'].set_text('$n_\\mathrm{{Fe^{{3+}}}}= {:.3e}$'.format(nFe)) texts['SCN'].set_text('$n_\\mathrm{{SCN^{{-}}}}= {:.3e}$'.format(nSCN)) texts['FeSCN'].set_text('$n_\\mathrm{{FeSCN^{{2+}}}}= {:.3e}$'.format(nFeSCN)) texts['Feeq'].set_text('$n_\\mathrm{{Fe^{{3+}}}}= {:.3e}$'.format(nFe)) texts['SCNeq'].set_text('$n_\\mathrm{{SCN^{{-}}}}= {:.3e}$'.format(nSCN)) texts['Feeq'].set_text('$n^{{\\mathrm{{eq}}}}_\\mathrm{{Fe^{{3+}}}}= {:.3e}$'.format(nFeeq)) texts['SCNeq'].set_text('$n^{{\\mathrm{{eq}}}}_\\mathrm{{SCN^{{-}}}}= {:.3e}$'.format(nSCNeq)) texts['FeSCNeq'].set_text('$n^{{\\mathrm{{eq}}}}_\\mathrm{{FeSCN^{{2+}}}}= {:.3e}$'.format(nFeSCNeq)) texts['x'].set_text('$x= {:.3e}$'.format(xi)) texts['xeq'].set_text('$x_\\mathrm{{eq}}= {:.3e}$'.format(xieq)) texts['FeSCNeq'].set_text('$n^{{\\mathrm{{eq}}}}_\\mathrm{{FeSCN^{{2+}}}}= {:.3e}$'.format(nFeSCNeq)) texts['Q'].set_text('$Q={:.3e}$'.format(Q)) #updating rectangles rects['Fe0'].set_width(nFe0/nmax*0.12) rects['SCN0'].set_width(nSCN0/nmax*0.12) rects['FeSCN0'].set_width(nFeSCN0/nmax*0.12) rects['Fe'].set_width(nFe/nmax*0.12) rects['SCN'].set_width(nSCN/nmax*0.12) rects['FeSCN'].set_width(nFeSCN/nmax*0.12) rects['Feeq'].set_x(0.48+nFeeq/nmax*0.12) rects['SCNeq'].set_x(0.68+nSCNeq/nmax*0.12) rects['FeSCNeq'].set_x(0.88+nFeSCNeq/nmax*0.12) fig.canvas.draw_idle() if __name__ == "__main__": nmax = 1e-2 parameters = { 'nFe0' : widgets.FloatSlider(value=0.05, description='$n^0_{\\mathrm{Fe^{3+}}}$', min=0.0, max=nmax), 'nSCN0' : widgets.FloatSlider(value=0.05, description='$n^0_{\\mathrm{SCN^{-}}}$', min=0.0, max=nmax), 'nFeSCN0' : widgets.FloatSlider(value=0, description='$n^0_{\\mathrm{FeSCN^{2+}}}$', min=0.0, max=nmax), 'xi' : widgets.FloatSlider(value=0.00000001, description='$x$', min=-0.01, max=0.01), } fig = plt.figure(figsize=(12,6)) ax = fig.add_axes([0,0,1,1]) ax.xaxis.set_visible(False) ax.yaxis.set_visible(False) gs = fig.add_gridspec(1, 2,left=0.05, right=0.95, bottom=0.05, top=0.8) ax1 = fig.add_subplot(gs[0,0]) fig.text(0.52, 0.8, '$\\mathrm{Fe^{3+}}$', fontsize=12,ha='center') fig.text(0.62, 0.8, '$+$', fontsize=12,ha='center') fig.text(0.72, 0.8, '$\\mathrm{SCN^-}$', fontsize=12,ha='center') fig.text(0.82, 0.8, '$=$', fontsize=12,ha='center') fig.text(0.92, 0.8, '$\\mathrm{FeSCN^{2+}}$', fontsize=12,ha='center') fig.text(0.1, 0.7, '$K={:.3e}$'.format(10**(2.3)), fontsize=12) texts = {} texts['x']= fig.text(0.55, 0.65, '$x=$', fontsize=12,ha='center') texts['xeq']= fig.text(0.55, 0.45, '$x_{\\mathrm{eq}}=$', fontsize=12,ha='center') texts['Fe0']= fig.text(0.52, 0.7, '$n^0_\\mathrm{Fe^{3+}}=$', fontsize=12,ha='center') texts['Fe']= fig.text(0.52, 0.6, '$n_\\mathrm{Fe^{3+}}=$', fontsize=12,ha='center') texts['Feeq']= fig.text(0.52, 0.4, '$n^{\\mathrm{eq}}_\\mathrm{Fe^{3+}}=$', fontsize=12,ha='center') texts['SCN0']= fig.text(0.72, 0.7, '$n^0_\\mathrm{SCN^{-}}=$', fontsize=12,ha='center') texts['SCN']= fig.text(0.72, 0.6, '$n_\\mathrm{SCN^{-}}=$', fontsize=12,ha='center') texts['SCNeq']= fig.text(0.72, 0.4, '$n^{\\mathrm{eq}}_\\mathrm{SCN^{-}}=$', fontsize=12,ha='center') texts['FeSCN0']= fig.text(0.92, 0.7, '$n^0_\\mathrm{FeSCN^{2+}}=$', fontsize=12,ha='center') texts['FeSCN']= fig.text(0.92, 0.6, '$n_\\mathrm{FeSCN^{2+}}=$', fontsize=12,ha='center') texts['FeSCNeq']=fig.text(0.92, 0.4, '$n^{\\mathrm{eq}}_\\mathrm{FeSCN^{2+}}=$', fontsize=12,ha='center') texts['Q']=fig.text(0.1, 0.6, '$Q=$', fontsize=12) transp = 0.5 rects={} rects['Fe'] = ax.add_patch(patches.Rectangle((0.48,0.55),0.5*nmax/nmax*0.12,0.02, edgecolor = '#000000', facecolor = 'white', fill=True,transform=fig.transFigure)) rects['SCN'] = ax.add_patch(patches.Rectangle((0.68,0.55),nmax/nmax*0.12,0.02, edgecolor = '#000000', facecolor = 'white', fill=True,transform=fig.transFigure)) rects['FeSCN'] = ax.add_patch(patches.Rectangle((0.88,0.55),nmax/nmax*0.12,0.02, edgecolor = '#000000', facecolor = '#aa0000', fill=True,transform=fig.transFigure)) rects['Fe0'] = ax.add_patch(patches.Rectangle((0.48,0.55),nmax/nmax*0.12,0.02, edgecolor = '#000000', facecolor = '#dddddd', fill=True,transform=fig.transFigure,alpha=transp)) rects['SCN0'] = ax.add_patch(patches.Rectangle((0.68,0.55),nmax/nmax*0.12,0.02, edgecolor = '#000000', facecolor = '#dddddd', fill=True,transform=fig.transFigure,alpha=transp)) rects['FeSCN0'] =ax.add_patch(patches.Rectangle((0.88,0.55),nmax/nmax*0.12,0.02, edgecolor = '#000000', facecolor = '#dddddd', fill=True,transform=fig.transFigure,alpha=transp)) rects['Feeq'] = ax.add_patch(patches.Rectangle((0.48,0.545),0,0.03, edgecolor = '#000000', facecolor = 'black', fill=True,transform=fig.transFigure)) rects['SCNeq'] = ax.add_patch(patches.Rectangle((0.68,0.545),0,0.03, edgecolor = '#000000', facecolor = 'black', fill=True,transform=fig.transFigure)) rects['FeSCNeq']=ax.add_patch(patches.Rectangle((0.88,0.545),0,0.03, edgecolor = '#000000', facecolor = 'black', fill=True,transform=fig.transFigure)) lines = {} #Q for all values of xis lines['Q'], = ax1.plot([],[],color='#cccccc') #equilibrium value of Q lines['eq'], = ax1.plot([],[],marker='o',color='#cccccc') #value of Q for given xi lines['Qsingle'], = ax1.plot([],[],marker='o',color='red') param_widgets = widgets.make_param_widgets(parameters, plot_data,slider_box=[0.07, 0.85, 0.85, 0.15]) #print(dir(param_widgets)) #print(param_widgets['xi']) #choose_widget = widgets.make_choose_plot(lines, box=[0.015, 0.25, 0.2, 0.15]) #reset_button = widgets.make_reset_button(param_widgets) #plt.tight_layout() plt.show()