%%FY1005/TFY4165, Øving 10, oppgave d: Plotting av van der Waals
%%isotermer for kjølemediet R134a, 1,1,1,2-tetrafluoretan.
%%R = gasskonstanten = 8.314 J/mol*K
R=8.314;
%%T = absolutt temperatur (K)
Tmin=334;
Tmax=414;
DeltaT=20;
%%T = vektor med verdier mellom Tmin og Tmax, intervall DeltaT
T=Tmin:DeltaT:Tmax;
%%V = molart volum (m^3)
Vmin=0.12e-3;
Vmax=2.00e-3;
NV=1000;
%%V = vektor med verdier mellom Vmin og Vmax, i alt NV verdier
V=linspace(Vmin,Vmax,NV);
%%Verdier for a og b for R134a: 1.00466  9.5734e-5
%%Enheter: [a] = Pa*m^6 og [b] = m^3
a=1.00466;
b=9.5734e-5;
%%length(T) = antall elementer i vektoren T
%%Bruker for-løkke fra i=1 til i=length(T) til å regne ut van der Waals
%%isotermer p(V) for temperaturer T(1), T(2), ... , T(length(T))
for i = 1:length(T);
    %%p = trykket.
    p = R*T(i)./(V-b) - a./(V.*V);
    fig = plot(V*1e3,p*1e-6);  %gir enheter L og MPa
    %%p(V) for laveste valgte temperatur T(1): blå kurve
    %%p(V) for høyeste valgte temperatur T(length(T)): rød kurve
    %%Mellomliggende kurver: gradvis mellom blå og rød
    %%Tynne kurver, LineWidth = 1.0, for van der Waals tilstandsligning
    red=(i-1)/(length(T)-1);
    blue=1-red;
    green = 0.0;
    set(fig,'Color',[red green blue],'LineWidth',1.0);
    if i == 1;
        title('van der Waals isotermer for R134a','fontsize',18);
        xlabel('Molart volum (L)','fontsize',18);
        ylabel('Trykk (MPa)','fontsize',18);
        axis([0 1 0 10]);
        %%"hold on" sørger for at påfølgende kurver tegnes i samme figur
        hold on;
        %%Vi plotter p(V) for ideell gass for laveste
        %%verdi av temperaturen, dvs for T(1).
        %%Dette blir kurve nr 2 i listen for legend nedenfor.
        pideell = R*T(i)./V;
        fig = plot(V*1e3,pideell*1e-6);
        %%Tykk blå kurve for p(V) ved T(1), ideell gass
        set(fig,'LineWidth',1.5,'Color',[0 0 1]);
    end;
    if i == length(T);
        %%Vi plotter p(V) for ideell gass for høyeste
        %%verdi av temperaturen, dvs for T(length(T)).
        %%Dette blir siste kurve i listen for legend nedenfor.
        pideell = R*T(i)./V;
        fig = plot(V*1e3,pideell*1e-6);
        %%Tykk rød kurve for p(V) ved T(length(T)), ideell gass
        set(fig,'LineWidth',1.5,'Color',[1 0 0]);
    end;
    legend('334 K','334 K (ideell gass)','354 K','374 K (T_c)','394 K' ...
        ,'414 K','414 K (ideell gass)');
    grid on;
end;
hold off;