    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
  