function u=PiS(T,rho,omega,mu0)

u=(24.5*T-1400)*(4.44-0.78*omega)*rho/mu0;
% u=(24.5*T-1400)*(3.8-0.78*omega)*rho/mu0;