function er = PR(V,T,P,am,bm)
global Rgas
er = P - ( Rgas*T/(V-bm) - am / ( V^2 + 2*bm*V - bm^2));

