%%FY1005/TFY4165, Oving 10, Oppgave 1b: van der Waals koeksistenstrykk
%%
%%Programmet tegner opp p(V) og ber om initialverdi for vaeskevolumet. 
%%Deretter beregnes gassvolumet som gir samme trykk.
%%Endelig beregnes integralet av p(V) numerisk og sammenlignes med
%%arealet under den rette linjen.
%%
%%Velg utskrift med flere gjeldende siffer enn det som er standard
format long;
%%Temperatur i Kelvin
T=123;
%%V = molart volum (L/mol)
Vmin=0.05;
Vmax=0.5;
NV=500;
%%V = vektor med verdier mellom Vmin og Vmax, i alt NV verdier
V=linspace(Vmin,Vmax,NV);
%%Verdier for a og b for luft: 1.368  0.0367
%%Enheter: [a] = bar*(L/mol)^2 og [b] = L/mol
R=8.134;
a=1.368;
b=0.0367;

%%van der Waals tilstandsligning:
p = (R*T./(V-b))/100 - a./(V.*V);
%%Plott av isotermen
fig = plot(V,p);
axis([0 Vmax 0 70]);
%%Tekst på figuren
title('Isoterm: van der Waals gass','fontsize',18);
xlabel('Molart volum (L/mol)','fontsize',18);
ylabel('Trykk (bar)','fontsize',18);
%%Med hold on; vil plottene "overleve", saa vi kan sammenligne.
hold on;
%%Les inn V1:
V1 = input('Gjett på V1:  ')
%%Riktig verdi er omtrent V1=0.06724129;
%%f.eks. 0.06724 gir 4*E-4 i feil for T=123K.
%%Bestemmer tilsvarende P1:
P1 = (R*T/(V1-b))/100 - a/(V1*V1);
%%Her har vi et volum som er mye storre enn V2:
V_stor=0.50;
%%Setter steglengden dV "liten nok"
dV=0.000001;
%%Bestemmer trykket for volumet V_stor
p_find=(R*T/(V_stor-b))/100 - a/(V_stor*V_stor);
%%Betingelsen V_stor > 0 gjor for at programmet stopper hvis noe gaar galt
while ( p_find < P1 && V_stor > 0)
	V_stor = V_stor-dV;
	p_find=(R*T/(V_stor-b))/100 - a/(V_stor*V_stor);
end;
%%Naa har vi funnet P2, V2, i den siste verdien:
V2=V_stor;
P2=p_find;
%%Skriver ut trykk- og volumverdier:
disp('V1: ')
disp(V1)
disp('V2: ')
disp(V2)
disp('P1: ')
disp(P1)
disp('P2: ')
disp(P2)
%%Estimat av usikkerhet i trykket P2:
dp = (- (R*T/((V2-b)*(V2-b)))/100 + 2.0*a/(V2*V2*V2))*dV;
disp('Estimert usikkerhet i trykket P2: ')
disp(abs(dp))
disp('Faktisk feil i P2: ')
disp(abs(P2-P1))
%%Tegner en linje paa plottet:
plot([V1,V2],[P1,P2]);
%%Saa skal vi integrere trykket fra V1 til V2, og sammenlikne med 
%%P*(V2-V1), dvs. arealet under linja. For aa integrere med quad() maa 
%%vi forst definere en symbolsk funksjon F:
%%"@(v)" sier at v er variabelen vaar, resten er likningen for P(v).
F = @(v) ((R*T)./(v-b))/100-a./(v.*v); 
%%Numerisk integrasjon av F fra V1 til V2 med saakalt adaptiv
%%Simpson-kvadratur.
I = quad(F, V1, V2); 
%%
disp('Feilen er:')
disp(I-P1*(V2-V1))
%%Tegner en linje paa plottet:
plot([V1,V2],[P1,P2]);
disp('Kjoer programmet paa nytt med ny gjetning for V1 dersom')
disp('feilen er for stor.')
