#!/usr/athena/bin/perl

$gaussian = 1;

$|=1;
$res = 30;
$ymin = $xmin = 100000;
$xmax = $ymax = 0;
while(<>){
    split;
    push(@x, $_[0]);
    push(@y, $_[1]);
    push(@z, $_[2]);
    if($_[0] > $xmax){
	$xmax = $_[0];
    }
    elsif($_[0] < $xmin){
	$xmin = $_[0];
    }
    if($_[1] > $ymax){
	$ymax = $_[1];
    }
    elsif($_[1] < $ymin){
	$ymin = $_[1];
    }
    $ct++;
    if(($ct%500) == 0){
	print STDERR "Loaded $ct points...\n";
    }
}

open(X, ">/tmp/xmesh");
open(Y, ">/tmp/ymesh");
open(C, ">/tmp/conf");

print STDERR "Done reading ", $#x+1, " data points.\n";
for $exct (0..$res) {
    $ex = ((($xmax-$xmin)/$res)*$exct)+$xmin;
    for $eyct (0..$res) {
	$ey = ((($ymax-$ymin)/$res)*$eyct)+$ymin;

	$tw = 0;
	$wa = 0;
	print STDERR "Calculating hyperbolicly weighted average at $exct($ex), $eyct($ey)\n";
	for $p (0..$#x){
	    $dx = $x[$p] - $ex;
	    $dy = $y[$p] - $ey;
	    $dsq = ($dx*$dx)+($dy*$dy);
#	    print STDERR "Dsq($p) --> $dsq ($dx, $dy)\n";
	    if($gaussian){
		$g = exp(-10000*$dsq);
		$tw+=$g;
		$wa+=$z[$p]*$g;
	    }
	    else{
		$tw += 1/$dsq;
		$wa += $z[$p]/$dsq;
	    }
	}
	print C "$wa ";
	$wa = $wa/$tw;
	print X "$ex ";
	print Y "$ey ";
	print "$wa ";
    }
    print "\n";
    print X "\n";
    print C "\n";
    print Y "\n";
}

open(SCALED, ">/tmp/scaleddata");
for $p (0..$#x){
    $x = ($x[$p]-$xmin)*$res/($xmax-$xmin);
    $y = ($y[$p]-$ymin)*$res/($ymax-$ymin);
    print SCALED "$x $y $z[$p]\n";
}
close(SCALED);