(*  Isothermen . Beispiel fuer (Subscript[H, 2])O . nach Kripfganz/Perlt,
   zur Physik siehe Landau,Achieser,Lifschitz S.172 *)

kb=UnitConvert[Quantity["BoltzmannConstant"]]
nA=UnitConvert[Quantity["AvogadroConstant"]]
RGsi=kb  nA
UnitConvert[Quantity["MolarGasConstant"]]
(*  Konversion von  kg m^2/s^2  in atm  geht mit:
9.86923 10^-3 , damit erhaelt man RG in atm l^3/(K mol) *)
RG = RGsi[[1]] * 9.86923*(10^(-3))
ct=273.15  (* Verschiebung Kelvin in Celsius *)

(* Van der Waals Gas, Zustandsgleichung *)
Clear[p,a,b,V,T]

p[a_,b_,V_,T_]:= RG (T+ct)/(V-b) -a/V^2

(*  Fuer (Subscript[H, 2])O Gas: a = 5.58 atm l^2 mol^-2
               b = 0.031 l mol^-1 und V>b ist vorausgesetzt *)

(* Umstellen der Gleichung zeigt fuer V eine Gl.3.Grades *)

Expand[p[a,b,V,T]*(V-b)*V^2- RG (T+ct) *V^2 +a*(V-b)]

(*  Eine oder drei Loesungen sind reell, die Stelle am
Uebergang ist der kritische Wert. *)

(* Suche kritische Temperatur *)
dpv=D[p[a,b,V,T],V]
dpVV=D[p[a,b,V,T],{V,2}] // Chop

eq=Eliminate[{dpv==0,dpVV==0},T] // FullSimplify

Solve[eq[[5]],V]

Solve[eq[[4]],V]//Chop // FullSimplify  (*  very crazy ... *)

Solve[eq[[3]],V]//Chop // FullSimplify (* etwas crazy ... *)

(*  verwende V = 3 b  aus der 1.Zeile weiter, da V>b ist, berechne kritische Temperatur: *)

loesKrit=Solve[(dpv /. V -> 3. b)==0,T]

Tkrit=loesKrit[[1,1]] /. {a-> 5.58, b-> 0.031}

Dann ist entsprechend   V_krit  -> 0.093  und  p_krit-> 218.5

a=5.58;  b=0.031;

(* Erstelle Diagramm *)
temps={200,300,376,400,500,600};
volStart={0.035,0.036,0.04,0.045,0.06,0.09};
Table[pp[i]=Plot[ p[a,b,V,temps[[i]] ],{V,volStart[[i]],0.5},PlotRange->{-240,620}],{i,1,6}]

poi=Graphics[{PointSize[Large],Point[{0.093,218}]}];

Show[pp[1],pp[2],pp[3],pp[4], pp[5],pp[6],poi,
PlotRange->{-250,620},AxesLabel->{"V(Liter)","p(atm)"},AxesOrigin->{0,0},PlotLabel->"Isothermen, T(Celsius)",Epilog->{Text["200",{0.095,-200}],
Text["600",{0.16,400}],Text["377",{0.092,180}] },
FormatType->TraditionalForm, TextStyle->{FontFamily->"Arial",FontSize->20}]

(* Globales Bild *)

Plot3D[ p[a,b,V,T],{V,0.04,0.3},{T,200,500},PlotRange->{-200,600},AxesLabel->{V,T,p},PlotPoints->50]

conts={-100,0,100,150,200,218.5,250,300,400};

ContourPlot[ p[a,b,V,T],{V,0.05,0.38},{T,200,475},Contours->conts,FrameLabel->{"Volumen","Temperatur"},ContourLabels->True]