#!/usr/local/bin/perl

if (@ARGV == 0) {
        die "usage: stats FILES [return]. The files should be in x,y format.\n";
}

### Subroutine 'numerically' for ordering the associative array
### by key order. If do simply a 'sort keys %array'
### will get a non-numerical, but alphabetical (even
### with numbers), ordering.

sub numerically { $a <=> $b }
# Note that $a and $b are "keywords" for this subroutine!!

#open(AVER,">average.out");
#open(STDDEV,">stddev.out");
open(STATS, ">result.out");

### First read in all data files, that should be
### in a x,y format, and put them in an associative
### array. The $key is the x and the field is the
### concatenation of all the y's.

print "reading... \n";

foreach $file (@ARGV) {
   open(FILE,"$file");
   while(<FILE>) {
	@tmp = split(' ',$_);
	$aa = shift @tmp;
	$bb = shift @tmp;
	$bb .= " ";
	$var{$aa} .= $bb;
   }
   close(FILE);
}
#foreach $key (sort numerically (keys(%var))) {
#   print $key, " ", $var{$key}, " \n";
#}  

### Calculate the averages of all the y's
### corresponding to a given x.

print "averages... \n";

foreach $key (keys(%var)) {
   @tmp2 = split(' ',$var{$key});
   $average{$key} = 0;
   foreach $myvar (@tmp2) {
	$average{$key} += $myvar;
   }
   $average{$key} /= ($#tmp2+1);
#  print $key," ", $average{$key},"\n";
}
#foreach $key (sort numerically (keys(%average)) ) {
##  print AVER $key, " ", $average{$key}, "\n";
#   print $key, " ", $average{$key}, "\n";
#}

### Calculate standard deviation.

print "standard dev... \n";

foreach $key (keys(%var)) {
   @tmp3 = split(' ',$var{$key});
   $stddev{$key} = 0;
   foreach $myvar (@tmp3) {
	$stddev{$key} += ($myvar - $average{$key}) * ($myvar - $average{$key});
   }
   $stddev{$key} /= ($#tmp3 + 1);
   $e = $stddev{$key};
   $stddev{$key} = sqrt($e);
#  print $key, " ", $stddev{$key}, "\n";
}	
#foreach $key (sort numerically (keys(%stddev)) ) {
##  print STDDEV $key, " ", $stddev{$key}, "\n";
#   print $key, " ", $stddev{$key}, "\n";
#}

### Print in single file in the x,y,dy format. The
### error bar then is in the -dy .. dy range.

foreach $key (sort numerically (keys(%stddev)) ) {
   print STATS $key, " ", $average{$key}, " ", $stddev{$key}, "\n";
}

#close(AVER);
#close(STDDEV);
close(STATS);