function [sigma]=sigma(T) 
  %surface tension [N/m] at T [celsius] 
 
 Tsat=[0:10:370,374.15];
  %saturation temperature in celsius

 surftens=1e-3*[75.6,74.24,72.78,71.23,69.61,67.93,66.19,64.40,62.57,60.69,58.78,56.83,54.85,52.83,50.79,48.70,46.59,44.44,42.26,40.05,37.81,35.53,33.23,30.90,28.56,26.19,23.82,21.44,19.07,16.71,14.39,12.11,9.89,7.75,5.71,3.79,2.03,0.47,0];

sigma=interp1(Tsat,surftens,T);