#!/usr/local/bin/perl

if (@ARGV == 0) {
        die "usage: $0 X_AXIS FILES [return]. \n\tThe FILES should be in x,y format.\n\tThe X-AXIS file is a single column of the x interpolation values.\n";
}

open(STATS, ">result.out");

# 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.
# NOTE THAT $A AND $B ARE "KEYWORDS" FOR THIS SUBROUTINE!!

sub numerically { $a <=> $b }

# FIRST READ THE X-AXIS TEMPLATE, OVER WHICH 
# THE INTERPOLATION WILL BE CARRIED

$first_x = shift @ARGV;
open(FILE,"$first_x");
$counter = 0;
while (<FILE>) {
	chop;
	$x_axis[$counter] = $_;
	$counter += 1;
}
close(FILE);

# SECOND READ IN ALL DATA FILES, THAT SHOULD BE IN A X,Y
# FORMAT, AND PUT THEM IN AN ASSOCIATIVE ARRAY AFTER
# INTERPOLATING THE READ X,Y INTO THE INPUT X-AXIS.  THE $key
# IS THE X AND THE FIELD IS THE CONCATENATION OF ALL THE Y'S.

print "reading... \n";

foreach $file (@ARGV) 
{
   print "... doing ",$file,"\n";
   open(FILE2,"$file");
   $counter = 0;
   while(<FILE2>) 
   {
     if (/^\#/)
     {
     }
     else
     {
	chop;
	@line = split(' ',$_);
 	$currx[$counter] = shift @line;
 	$curry[$counter] = shift @line;
	$counter += 1;
     }
   }
   close(FILE2);

   foreach $i (0 .. $#x_axis)
   {
      # Find the bounding array elements.
      $success = 0;
      for ($j=0; $j<$#currx; $j++)
      {
	if ( ($currx[$j]  <=$x_axis[$i]) 
	  && ($currx[$j+1]>=$x_axis[$i]) 
 	  && ($success != 1) )
	{
		$lw      = $j;
		$up      = $j+1;
	 	$success = 1;
	}
      }

      # INTERPOLATE
      if ($success==1)
      {
      	$interp_y  = ($curry[$up]-$curry[$lw])*($x_axis[$i]-$currx[$lw])/($currx[$up]-$currx[$lw]) + $curry[$lw];
      	$interp_y .= " ";
      	$var{$x_axis[$i]} .= $interp_y;
      }
   }
}

# 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);
}

# 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 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(STATS);