[t, p] = meshgrid(-40:40, 500:10:1000);
q = zeros(size(t));
[rh, rw, ew_salt] = get_rh(q, t, p);

figure
cont = [1 2:2:20 30 40 50];
% [ c, h] = contourf(t, p, 1000*rw, cont);
% colorbarf('vert', cont, cont)
[c, h] = contour(t, p, 1000*rw, cont);
clabel(c, h)
axis('ij')
grid
xlabel('temperature')
ylabel('pressure')
titstr{1} = 'saturation mixing ratio in gm/kg with Kraus correction';
titstr{2} = ['rw (g/kg) = ' num2str(1000*rw(end, 61)) ' at p = ' ...
	    num2str(p(end, 61)) ' & t = ' num2str(t(end, 61))];
title(titstr)

str = ['rw (g/kg) = ' num2str(1000*rw(end, 61)) ' at p = ' ...
	    num2str(p(end, 61)) ' & t = ' num2str(t(end, 61))];
disp(str)
disp(['ew_salt (mb) = ' num2str(ew_salt(end, 1)) ' at t = ' ...
      num2str(t(end, 1))]);
disp(['ew_salt (mb) = ' num2str(ew_salt(end, 41)) ' at t = ' ...
      num2str(t(end, 41))]);
disp(['ew_salt (mb) = ' num2str(ew_salt(end, 61)) ' at t = ' ...
      num2str(t(end, 61))]);
disp(['ew_salt (mb) = ' num2str(ew_salt(end, 81)) ' at t = ' ...
      num2str(t(end, 81))]);
