    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %    Calculate the Saturation      %
    %    Vapour Pressure using the     %
    %    Clausius-Clapeyron Equation   %
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%   Define the temperature
%%%%%%%%%   Range and increment.

NN = 61;
Tmin = -20;    %  Minimum temperature (Celcius)
Tmax = +40;    %  Maximum temperature (Celcius)
Tdiff = (Tmax-Tmin)/(NN-1);   %   Temperature increment.

%%%%%%%%  Define a few constants
Lv = 2.5*10^6;     %   Latent heat of vapourization
Rv = 461;          %   Gas constant for water vapour
T0 = 273;          %   Temperature offset for Kelvin scale.
es0 = 6.11;        %   Saturation vapour pressure at 0 deg C.

%%%%%%%%   Loop through the Tmperature values.
for N = 1:NN
  T = Tmin + (N-1)*Tdiff + T0;
  x = (Lv/Rv)*((1/T0)-1/T);
  es = es0*exp(x);
  TEMP(N) = T-T0;    %   Change temperature to Celcius
  ESAT(N) = es;      %   Store es values in array
end

%%%%%%%%   Plot the results
axis([-20,40,0,80])
plot(TEMP,ESAT)

title('Saturation Vapour Pressure')
xlabel('Temperature')
ylabel('e_s (g/kg)')

%%%%%%%%%%%%%%%  Fancy bit: shade below and above curve.
hold on
fill([TEMP,Tmax,Tmin],[ESAT,0,0],'r')  
fill([TEMP,Tmax,Tmin,TEMP(1)],[ESAT,80,80,ESAT(1)],'c')  

return
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
