#!/usr/bin/perl
# the above is the standard linux/unix path for perl: change if your computer is nonstandard

sub print_help{
    print STDERR <<EOF;

Syntax:  xyz_bond_length [options] xyzfile1 [xyzfile2, ...]

write out the length of one bond at all frames of the xyz file(s)
and report the second string in the label line

Options:
-a a1 a2	select the atoms to compute distance (default: 1 and 2)
							by Nick 21 Jun 2011
EOF
    return;
}

$BIG=99.e99;
$atom1=1;
$atom2=2;

while($p=shift(@ARGV)){
    if($p=~/^[^\-]/){ # all command line options start with "-", otherwise it a filename
        if(-e $p){
            push(@inputfiles,$p);
        }else{
            print STDERR "ERROR: file $p does not exist\n";
        }
        next;
    }
    if($p eq "-a"){ # help
	$atom1=shift(@ARGV)+0;
	$atom1>0 or die "Invalid atom index $atom1\n";
	$atom2=shift(@ARGV)+0;
	$atom2>0 or die "Invalid atom index $atom2\n";
	next;
    }
    if($p eq "-h"){ # help
	print_help;
	die "\n";
    }
    print_help;
    die "invalid option $p\n";
}

if($#inputfiles<0){
    push(@inputfiles,"-");      #       adds stdio for lack of arguments
}

@ave=(0.,0.,0.);
@ave2=(0.,0.,0.);
@max=(-$BIG,-$BIG,-$BIG);
@min=($BIG,$BIG,$BIG);

while ($filename=shift(@inputfiles)){
    for($i=0;$i<3;$i++){
	$di[$i]=0;
    }
    $count=-1;
    $frames=0;
    open(INPUT, $filename) || die "Cannot open $filename\n";
    while(<INPUT>){
	s/^\s*//;	# remove leading spaces
	@li = split('\s+');
	if($#li==0 && $count!=0){
	    $count=-1;
	    $frames++;
	    if($di[0]!=0 || $di[1]!=0 || $di[2]!=0){ # previous frame...
		$distance=0.;
		for($i=0;$i<3;$i++){
		    $distance+=$di[$i]*$di[$i];
		}
		print $distance."\n";
	    }
	    for($i=0;$i<3;$i++){
		$di[$i]=0;
	    }
	}
	if($count==0){
	    print $frames." ".$li[1]." ";
	}elsif($count>0){
	    if($count==$atom1){
		for($i=0;$i<3;$i++){
		    $di[$i]+=$li[$i+1];
		}
	    }
	    if($count==$atom2){
		for($i=0;$i<3;$i++){
		    $di[$i]-=$li[$i+1];
		}
	    }
	}
	$count++;
    }
    close(INPUT);
    $distance=0.;  # last frame...
    for($i=0;$i<3;$i++){
	$distance+=$di[$i]*$di[$i];
    }
    print $distance."\n";
}