#!/usr/bin/perl
#
use POSIX;

sub print_help {
    print STDERR <<EOF;
usage:	intepolate filename [-i interval] [-m]
		generate a simple linear interpolation
		of a N-column data file x f1(x) f2(x) f3(x) ... 
	-i sets the interpolation interval
	-m outputs only the mesh points, suppressing original data
output: N columns with:                 x interpol._f1(x) interpol._f2(x) ...
			      			by Nick Manini, 13-10-2009
EOF
return;
}

$BIG=9.99e99;
$interval=1;
$originaldata=1;
$xmin=-$BIG;
$xmax=$BIG;

while($p=shift(@ARGV)){
    if($p=~/^[^\-]/){#filename unless starting with -
        if(-e $p){
            push(@inputfiles,$p);
        }else{
            print STDERR "File $p does not exist!!!\n";
        }
        next;
    }
    if($p eq "-h"){# [-h] print help
        print_help;
        exit;
    }
    if($p eq "-i"){# [-i int] forces interval int
        $interval = shift(@ARGV);
        next;
    }
    if($p eq "-m"){# [-m] forces to print out only mesh points
        $originaldata=0;
        next;
    }
    if($p eq "-r"){# [-r xmin xmax] forces range
        $xmin = 0+shift(@ARGV);
        $xmax = 0+shift(@ARGV);
	$xmin < $xmax || die "interpolate error: range $xmin - $xmax not valid\n";
        next;
    }
    print_help;
    print STDERR "invalid option $p!\n ";
    exit;
}

if($#inputfiles < 0){
   $filename="/tmp/interpolate_tmp".$$;
   `cat > $filename`;
   push(@inputfiles,$filename);
   $rmflag=1;
}

$mainfile="/tmp/interpolate".$$;
`cat "@inputfiles" > $mainfile`;# put everything in a single file

$hold[0]="notyet";
open(INPUT,$mainfile) || die "cannot open file $mainfile\n";
while(<INPUT>){
    s/^\s*//;	# remove leading spaces
    if(!/^#/){	# ignore commented lines
	my @fi =  split(/\s+/);
	if($#fi > 0){	# process lines with at least 2 fields
	    if($hold[0] eq "notyet"){		
		for($j=0;$j<=$#fi;$j++){
		    $hold[$j]=$fi[$j];
		}
		if($xmin!=-$BIG){
		    $xstart=POSIX::floor($xmin/$interval)*$interval;
		}else{
		    $xstart=POSIX::floor($hold[0]/$interval)*$interval;
		}
		for($i=0;;$i++){
		    $x=$xstart+$interval*$i;
		    if($x == $fi[0] && $originaldata){$i++;}
		    $x < $fi[0] || last;
		    if($x>=$xmin && $x<=$xmax){
			print $x;
			for($j=1;$j<=$#hold;$j++){
			    print "\t",$hold[$j];
			}
			print "\n";
		    }
		}
		if($originaldata && $fi[0]>=$xmin && $fi[0]<=$xmax){print;}
	    }else{
		for(;;$i++){
		    $x=$xstart+$interval*$i;
#		    print "quiii ", $x," ", $fi[0]," ",$i,"\n";
		    if($x == $fi[0] && $originaldata){$i++;}
		    $x < $fi[0] || last;
		    if($x>=$xmin && $x<=$xmax){
			print $x;
			for($j=1;$j<=$#fi;$j++){
			    print "\t",(($hold[0]*$fi[$j]-$fi[0]*$hold[$j]) + $x*($hold[$j]-$fi[$j]))/($hold[0]-$fi[0]);
			}
			print "\n";
		    }
		}
		if($originaldata && $fi[0]>=$xmin && $fi[0]<=$xmax){print;}
		for($j=0;$j<=$#fi;$j++){
		    $hold[$j]=$fi[$j];
		}
	    }
	}
    }
}
if($xmax!=$BIG){
    for(;;$i++){
	$x=$xstart+$interval*$i;
	$x <= $xmax || last;
	print $x;
	for($j=1;$j<=$#hold;$j++){
	    print "\t",$hold[$j];
	}
	print "\n";
    }
}

`rm $mainfile`;
if($rmflag){`rm $filename`;}