function  e0_g = E0( T )
%function [ e0_g, e0_h ] = E0( T )
%E0 Standard electrode potential
% $Id: E0.m,v 1.2 2005/02/08 09:32:21 zenith Exp $

[ e0_g, e0_h ] = Gwater( T, 101.325 ) ;
e0_g = - e0_g / ( 2 * 96.485 ) ;
%e0_h = - e0_h / ( 2 * 96.485 ) 


function [ g, h ] = Gwater( T, p )
%Gwater Selects liquid or gas condition

if ( T < 0 )
    error( 'Absolute temperatures are positive!' ) ;
end

Pvapour = P0( T ) ;

if ( Pvapour < p )
    [ g, h ] = Gliq( T ) ;
else
    [ g, h ] = Ggas( T ) ;
end


function [ g, h ] = Gliq( T )
%Gliq Thermodynamical data for water in liquid state

if ( T < 298 | T > 500 )
    warning( 'Gliq function is defined only between 298K and 500K, expect bogus results' ) ;
end

t = T / 1000 ;
A = - 251.5433 ;
B = 1531.7293695 ;
C = - 3207.4290395 ;
D = 2477.25151 ;
E = 4.1248225 ;
F = - 241.6251335 ;
G = - 780.4783 ;

h = A * t + B * t^2 / 2 + C * t^3 / 3 + D * t^4 / 4 - E / t + F ;
s = A * log( t ) + B * t + C * t^2 / 2 + D * t^3 / 3 - E / ( 2 * t^2 ) + G ;
g = h - t * s ;


function [ g, h ] = Ggas( T )
%Ggas Thermodynamical data for water in gaseous state

if ( T < 298 | T > 1500 ) 
    warning( 'Ggas is defined only between 298K and 1500K, expect bogus results' );
end

t = T / 1000 ;
A = - 17.8445 ;
B = 15.2718835 ;
C = - 4.2226045 ;
D = 0.26203 ;
E = 0.3516355 ;
F = - 235.9583335 ;
G = - 68.3653 ;

h = A * t + B * t^2 / 2 + C * t^3 / 3 + D * t^4 / 4 - E / t + F ;
s = A * log( t ) + B * t + C * t^2 / 2 + D * t^3 / 3 - E / ( 2 * t^2 ) + G ;
g = h - t * s ;


function p = P0( T )
%P0 Vapour pressure of water
% This function calculates the vapour pressure of water
% between 0 and 170 °celsius.
% Temperature is expressed in kelvin.
% Pressure is expressed in kPa.
% Data from:	NIST (0-100°C)
%               http://www.s-ohe.com/water.html (100-170°C)
% $Id $

if ( T <= 303 )
	p = 10^2 * 10^( 5.40221 - ( 1838.675 / ( T - 31.737 ) ) ) ;

elseif ( T > 303 & T <= 333 ) 
	p = 10^2 * 10^( 5.20389 - ( 1733.926 / ( T - 39.485 ) ) ) ;

elseif ( T > 333 & T <= 363 )
	p = 10^2 * 10^( 5.07680 - ( 1659.793 / ( T - 45.854 ) ) ) ;

elseif ( T > 363 & T <= 373.15 )
	p = 10^2 * 10^( 5.08354 - ( 1663.125 / ( T - 45.622 ) ) ) ;

elseif ( T > 373.15 )
	p = exp( 0.0295815 * T - 6.42002 ) ;
end

if ( T < 273 | T > 443 )
    warning( 'The input temperature is outside the range of definition of function P0. Do not trust the results.' ) ;
end