
getf("../equilibre/thermochim.sci") 
data=readData("../equilibre/thermodata.txt");
[DH,DS]=reaction(data,"CH_4(g)+2*H_2O(g)=CO_2(g)+4*H_2(g)")
            
plotA=scf();
plotK(DH,DS,300,1500,100)
Ti=DH*1e3/DS
K1=constanteK(DH,DS,1100)

function y=avancement(T,P)
 x=poly(0,"x");
 p=(4*x)^4*x-constanteK(DH,DS,T)/P^2*(2-2*x)^2*(1-x)*(3+2*x)^2;
 r=roots(p);
 a=0
 for i=1:length(r),
    if imag(r(i))==0 then 
        if ((real(r(i))<1)&(real(r(i))>0)) then 
            a=r(i)
        end,
    end,
 end,
 y=a
endfunction
            
 
n=100; Tmin=500; Tmax=1500;
temp=[Tmin:(Tmax-Tmin)/(n-1):Tmax];
a1=temp;
a2=temp;
for i=1:n,
    a1(i)=avancement(temp(i),1); 
    a2(i)=avancement(temp(i),10);            
end
            
plotB=scf();
plot2d(temp,a1,style=21)
plot2d(temp,a2,style=9)
legend(['1 bar';'10 bar'])
xtitle('Avancement','T (K)','x')
x1=avancement(1100,1)

function T=Tfinale(Ti,x)
    T=Ti-x*DH*1e3/(35*(1-x)+33.58*2*(1-x)+37.11*x+29*4*x)
endfunction
Ti=1800
n=10;
av=[0:1/(n-1):1];
Tf=av;
for i=1:n,
    Tf(i)=Tfinale(Ti,av(i));
end

            
plotC=scf();
plot2d(Tf,av,style=21)
plot2d(temp,a1)
legend(['Bilan H','Equilibre'])
xtitle('Avancement','T(K)','x')
