#!/usr/bin/perl
# $Id: activation_energy,v 1.1 2006/03/28 08:06:45 skoch Exp $
use strict;
use warnings;

my $debug = 0;
my $prob = 0.95;

for my $string (@ARGV)
{
    if ($string =~ /--debug/)
    {
	$debug = 1;
    }
    if ($string =~ /--prob\s+[\d\.]+/)
    {
	$string =~ /--prob\s+([\d\.]+)/;
	$prob = $1;
    }
    if ($string =~ /--help/)
    {
	&print_help;
	exit;
    }
}


my @data = ();
while (!eof(STDIN))
{
    my $st = <STDIN>;
    chomp $st;
    push @data, $st;
}
my $sumx = 0;
my $sumy = 0;
my $sumxy = 0;
my $sumxx = 0;
my $sumyy = 0;
my $n = 0;
for my $string (@data)
{
    print $string,"\n" if ($debug);
    if ($string =~ /^\s*\#/)
    {
    }
    else
    {
	unless ($string !~ /[\d\.\-e\+]+\s+[\d\.\-\+e]+/) 
	{
	    $string =~ /([\d\.\-e\+]+)\s+([\d\.\-\+e]+)/;
	    my $x = $1;
	    my $y = $2;
	    $sumx += $x;
	    $sumy += $y;
	    $sumxy += $x * $y;
	    $sumxx += $x * $x;
	    $sumyy += $y * $y;
	    $n++;
	}
    }
}
print "sumx $sumx\nsumy $sumy\nsumxy $sumxy\nsumxx $sumxx\nsumyy $sumyy\n" if ($debug);
if ($n < 2)
{
    print "Less than 2 data points were found, exiting!.\n";
    exit 0;
}
my $ssdx = $sumxx - $sumx * $sumx / $n;
my $ssdy = $sumyy - $sumy * $sumy / $n;
my $spd = $sumxy - $sumx * $sumy / $n;

my $alfa = $spd / $ssdx;
my $beta = ($sumy - $alfa * $sumx)/$n;
my $ssd02 = $ssdy - $spd * $spd / $ssdx;
my $s02_2 = 0;
unless ($n < 3)
{
    $s02_2 = $ssd02 / ($n - 2);
}
print "ssdx $ssdx\nssdy $ssdy\nspd $spd\nalfa $alfa\nbeta $beta\nssd02 $ssd02\ns02_2 $s02_2\n" if ($debug);

my $k = 1.38E-23;
my $e = 1.602E-19;

my $Ea = - $alfa * 1000 *$k / $e;
my $var = '';
if ($n < 3)
{
    $var = 'NA';
}
else
{
    my $df = $n - 2;
    my $p = 1 - $prob /2 ;
    my $command = '/usr/local/bin/probdist prob t ';
    $command = $command . $df . ' '.$p;
    print $command ,"\n" if ($debug);
    my $v = `$command` or die "Could not execute $command\n";
    chomp $v;
    my $factor = sqrt($s02_2 / $ssdx);
    $var = $v * $factor *1000 * $k / $e;
}
if ($var =~ /NA/)
{
    printf "Ea: %6.4f eV  n: $n.\n", $Ea;
}
else
{
    printf "Ea: %6.4f \+\/\- %6.4f eV  n: $n.\n", $Ea, $var;
}
exit;

sub print_help
{
    print "activation_energy calculates the activation energy based on linear regression on formatted data\n";
    print "Usage activation_energy [opt. --debug --prob \$prob].\n";
    print "The data must be formatted as following: 1000K/T  ln(1/R).\n";
    print "The data has to be typed (or piped) into <STDIN> one line at a time.\n";
}
