N=(6*(10^(23)));
m1=(2*(1.67*(10^(-27))));
m2=((2*16)*(1.67*(10^(-27))));
m3=(18*(1.67*(10^(-27))));
p1=0.09;
p2=1.428;
p3=998.21;
nu1=2;
nu2=1;
nu3=-2;
p0=101325;
epsilon=(p0*N)*(((nu1*m1)/p1)+((nu2*m2)/p2)+((nu3*m3)/p3));
c=5;
R=8.3;
F=96485;
K=-273;
KR=293;
f3=0.019;%liquid fugacity of water, solvent, deviation from Raoult's law, assume independent of pressure/temperature%
f1=0.019;%liquid fugacity of hydrogen, solute, deviation from Henry's law, assume independent of pressure/temperature%
f2=0.019;%liquid fugacity of oxygen, solute, deviation from Henry's law, assume independent of presure/temperature%
d1=1;%gas fugacity of hydrogen, assume independent of pressure/temperature%
d2=1;%gas fugacity of oxygen, assume independent of pressure/temperature%
k1=0.00078; %slope in Henry's law for hydrogen%
k2=0.00128; %slope in Henry's law for oxygen%
v0=(18*10^(-6)); %molar volume of water%
fixed=1; %the fixed dc power supply potential, you have to input the actual value;
for i=1:1000
for j=1:1000
measure(i,j)=0.1;
end
end
%Theoretically, you have to measure this as the pressure P' in the liquid mixture for which the%
%gas pressure of water vapour on its own in the mixture at equilibrium is atmospheric pressure divided by the%
%liquid fugacity of water. Using Raoult's law with fugacity,%
%(P_3)'*=(P_3)'/f3*x3=p0/f3, so that (P_3)'=p0*x3~po,%
%where (P_3)'* is the gas pressure of water vapour on its own,%
%(P_3)' is the partial pressure of water vapour, x3~1 is the%
%mole fraction of water in the liquid mixture. Assuming that P'~(P_3)', we have that v0*P'~0.1.%
for i=1:1000
for j=1:1000
delta3(i,j)=(v0*(p0+(i/10)))-measure(i,j);
end
end
for i=1:1000
for j=1:1000
kappa1(i,j)=((R*(293+(j/10)))*log(((d1*k1)/p0)));
end
end
for i=1:1000
for j=1:1000
kappa2(i,j)=((R*(293+(j/10)))*log(((d2*k2)/p0)));
end
end
for i=1:1000
for j=1:1000
potential2(i,j)=2+((i/1000)*(j/1000));
end
end
%you have to measure this, the potential in a given range of temperatures
%and pressures with the dc power supply on%
for i=1:1000
for j=1:1000
potential(i,j)=potential2(i,j)-fixed;
end
end
%the true potential of the reaction%
for j=1:1000
potential0(j)=potential(1,j);
end
%the potential in a range of temperatures at atmospheric pressure p0%
for i=1:1000
for j=1:1000
potential1(i,j)=potential(i,j)-potential0(j);
end
end
%the difference in potentials given by potential and potential1%
for i=1:1000
for j=1:1000
phi3(i,j)=(((R*(293+(j/10)))*log(f3))+delta3(i,j));
end
end
for i=1:1000
for j=1:1000
psi1(i,j)=(((R*(293+(j/10)))*log(f1))+kappa1(i,j));
end
end
for i=1:1000
for j=1:1000
psi2(i,j)=(((R*(293+(j/10)))*log(f2))+kappa2(i,j));
end
end
for i=1:1000
for j=1:1000
epsilon1(i,j)=nu3*phi3(i,j)+nu1*psi1(i,j)+nu2*psi2(i,j);
end
end
for i=1:1000
for j=1:1000
activity(i,j)=exp(((epsilon*log((p0+(i/10))))-epsilon1(i,j))/(R*(293+(j/10))));
end
end