#!/usr/bin/perl
# $Id: leak,v 1.3 2006/09/21 06:14:07 skoch Exp $
my $VERSION = 0.2;
my $PROGRAM = 'leak';
use constant F => 96485;
use constant R => 8.314;
use constant p0 => 1;

use strict;

if (@ARGV[0] =~ /--help/i)
{
    &help;
    exit;
}

for (my $x = 0; $x < scalar(@ARGV) ; $x++)
{
    if ($ARGV[$x] =~ /--version/)
    {
	print "$PROGRAM is: $VERSION\n";
	exit;
    }
}


my @arg = @ARGV;
my %input = &get_input(@arg);
#print %input;

our $atm = defined($input{'gas'}) ? $input{'gas'} : 1;
our $watertemp = defined($input{'water'}) ? $input{'water'}+273.15 : 273.15;
our $celltemp = defined($input{'cell'}) ? $input{'cell'}+273.15 : 1123.15;
# Flow is measured in L/h !!!
our $H2flow = defined($input{'H2flow'}) ? $input{'H2flow'} : 0.1;
our $O2flow = defined($input{'O2_add'}) ? $input{'O2_add'} : 0;
our $dry = defined($input{'dry'}) ? $input{'dry'} : 0;

# Coefficients for the water partial pressure are obtained by fitting a 6-
# power function to data from 'Databog Fysik Kemi', F og K forlaget 10.udgave 
my @a;
$a[6] = 6.3096E-11;
$a[5] = 2.5259E-6;
$a[4] = -3.2193E-3;
$a[3] = 1.6322;
$a[2] = -414.9;
$a[1] = 52920;
$a[0] = -2.7085E6;
# End ref
my $ph2o = 0;
for (my $x=0; $x<7; $x++)
{
    my $test = $watertemp;
    for (my $y=1;$y<$x;$y++)
    {
	$test = $test* $watertemp;
    }
    if ($x==0){$test = 1;};
    $ph2o += $a[$x]* $test;
} 

$ph2o = 1.0E-10 if ($dry);

our $command = 0;
our $showall = 0;
if (defined($input{'command'}))
{
    if ($input{'command'} =~ /all/) 
    {
	$showall = 1;
    }
    if ($input{'command'} =~ /ph2o/)
    {
	$command = 1;
    }
    if ($input{'command'} =~ /OCV/)
    {
	$command = 2;
    }
    if ($input{'command'} =~ /leak/)
    {
	$command = 3;
    }
    if ($input{'command'} =~ /po2/)
    {
	$command = 4;
    }
    if ($input{'command'} =~ /w\/h/)
    {
	$command = 5;
    }
    if ($input{'command'} =~ /h\/w/)
    {
	$command = 6;
    }
    if ($input{'command'} =~ /I/)
    {
	$command = 7;
    }
}


#####  Formulas  #####

# log(Kp) = 2.958-13022/T (from Hartung, see below)
# P_O2 = Kp^2 * (P_H2O / P_H2)^2  according to eqlibrium-function 
# Emf = (RT/4F)*ln(P_O2'/P_O2'')  versus air, according to Nernst
# P * V = n * R * T  Ideal gas law
# I = 4 * n * F  Faradays law for electrolysis for O2
#####  OCV calculation #######  
$ph2o = $ph2o/100000;
#  log(Kp) is according to R. Hartung and H.-H. Mobius:
#  "Potentiometrische Bestimmung des Wasserdampf
#  Dissoziationsgleichgewichtes zwischen 1000 und
#  1300K mit einer Festelektrolytzelle",
#  Chemie-ing.-techn. 40-12 pp592-600, 1968 
our $logKp = 2.958 - 13022/$celltemp; 
#  End ref
our $ph2 = $atm;

if( defined($input{'pH2O'})) 
{ 
    $ph2o = $input{'pH2O'}+0;
}
our $Kp = 10**$logKp;
our $nH2 = $H2flow * $atm / (60 * 60 * 0.0821 * 273.15);
our $nH2O = $H2flow * $ph2o / (60 * 60 * 0.0821 * 273.15 * (1-$ph2o));

our $nO2 = $O2flow / (60 * 60 * 0.0821 * 273.15);
$nH2O = $nH2O + 2 * $nO2;
$nH2 = $nH2 - 2 * $nO2; 
our $ph2_ = $nH2 / ($nH2 + $nH2O);
our $ph2o_ = $nH2O / ($nH2 + $nH2O);
our $pO2 = $Kp * $Kp * $nH2O * $nH2O / ($nH2 * $nH2 );
our $latm = log(0.21 / $pO2);
our $OCV = -$latm * R * $celltemp / (4 * F);

our $Em = defined($input{'emf'}) ? $input{'emf'} : $OCV;
#### Leak calculation based on measured EMF ($Em) #####
our $po2m = 0.21* exp( 4* F * $Em / (R * $celltemp));
our $Z = sqrt($po2m)/$Kp;
our $t = ($Z * $nH2 - $nH2O) / (2 + 2 * $Z);
# $x is the leak in percent of the H2 flow! 
#### Leak current calculation, Current in A!!!   ####
our $deltanH = 2* $t;
our $x = 100 * $deltanH / $nH2; 
our $IO2 = 4 * F * $t;
our $pH2o = ($nH2-$deltanH)/($nH2+$nH2O);
our $pH2Oo = ($nH2O+$deltanH)/($nH2+$nH2O);

#### Output of calculations ########
if ($showall)
{
    print "OCV:\n";
    print "pH2O = ",$ph2o_,",  pH2 = ",$ph2_,",  pO2 = ",$pO2,"\n";
    print "OCV: ",$OCV," V\n";
    print "Real:\n";
    print "pH2O = ",$pH2Oo,",  pH2 = ",$pH2o,",  pO2 = ",$po2m,"\n";
    print "EMF: ",$Em," V\n";
    print "Oxygen leak = ",$x," % of H2 flow\n";
    print "Leak current = ",$IO2," A for a H2 flow of ",$H2flow," L/h\n";
    print "All calculations based on a cell temperature of ",$celltemp," K\n";
    print "and an OCV-water vapour pressure determined by eqlibrium\n";
    print "with liquid water at ",$watertemp," K\n";
};

if ($input{'command'} =~ /pH2o/)
{
    print "$pH2o\n";
    $command = -1;
}
if ($input{'command'} =~ /pH2Oo/ )
{
    print "$pH2Oo\n";
    $command = -1;
}

if ($command == 1)
{
    print $ph2o_,"\n";
    exit 1;
}
if ($command == 2)
{
    print $OCV,"\n";
    exit 1;
}
if ($command == 3)
{
    print $x,"\n";
    exit 1;
}
if ($command == 4)
{
    print $po2m,"\n";
    exit 1;
}
if ($command == 5)
{
    my $out = ($nH2O+2*$t)/($nH2-2*$t);
    print $out,"\n";
    exit 1;
}
if ($command == 6)
{
    my $out = ($nH2-2*$t)/($nH2O+2*$t);
    print $out,"\n";
    exit 1;
}
if ($command == 7)
{
    print $IO2,"\n";
    exit 1;
}

if ($command == 0)
{
    if (!$showall)
    {
	print "OCV = ",$OCV,",  Oxygen leak = ",$x," % of H2 flow \n";
	print "Leak current = ",$IO2," A\n";
    }
}
exit 1;

sub get_input
{
    my %out = ();
    my @input = @_;
    for (my $x = 0; $x < scalar(@input);$x++)
    {
	if ($input[$x] =~ /^*--/)
	{
	    $input[$x] =~ s/^*--//;
	    $input[$x+1] =~ s/^*\s//;
	    $out{$input[$x]} = $input[$x+1];
	    $x++;
	}
    }
    return %out;
}

sub help
{
    print "Leak  Usage:\n";
    print "leak --INPUT VALUE\n\n";
    print "  INPUT is the command for which the VALUE is the corresponding input value.\n\n";
    print "     --command : forces output into the following formats:\n";
    print "          all  : All values are displayed (verbose output).\n";
    print "          ph2o : Outputs only the water partial pressure on the anode side.\n";
    print "                 of an ideal cell at OCV (theoretical OCV).\n";
    print "          OCV  : Outputs only the open circuit voltage (theoretical).\n";
    print "          leak : Outputs only the oxygen leak (in \%) into the hydrogen.\n";
    print "          I    : Outputs leak current in A\n";
    print "          po2  : Outputs only the oxygen partia pressure on the anode side.\n";
    print "          w\/h  : Outputs only the ratio of the water and hydrogen partial\n";
    print "                 pressures on the anode side of the cell.\n";
    print "          h\/w  : Inverse of w\/h. \n\n";
    print "     --water   : inputs the water temperature, default is 25 C.\n\n";
    print "     --cell    : Inputs the cell temperature, default is 850 C.\n\n";
    print "     --emf     : Inputs the cell emf, default is OCV.\n\n";
    print "     --H2flow  : Inputs the hydrogen flow in L/h, default is 0.1.\n";
    print "     --O2_add  : Inputs the oxygen flow to the anode gas in L/h,\n";
    print "                 default is 0.\n";
    print "                 Note that this has to be at 0 degree C!!\n\n";
    print "     --gas     : Inputs the hydrogen gas composition.\n";
    print "                 Default is 1, coresponding to 100\% hydrogen\n";
    print "                 before equlibrium with liquid water.\n\n";
    print "     --pH2O    : Inputs the water partial pressure of the gas if this\n";
    print "                 this is not determined by equlibrium with liquid water.\n\n";
}
