#!/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_displacement [options] xyzfile1

write out another xyz file with
a) the displacements between subsequent frames of xyzfile1 (default behavior)
b) changed names of atoms moving more than a given threshold (when -t is set)

Options:
	-t threshold : set the displacement threshold to change the atom name
			by default threshold=0, so that output is displacements
							by Nick 10 May 2017
EOF
    return;
}

$BIG=99.e99;
$thresh=0;

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 xyz_displacement: file $p does not exist\n";
        }
        next;
    }
    if($p eq "-t"){ # help
	$thresh=shift(@ARGV)+0;
	$thresh>0 or die "Invalid threshold $thresh\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                     
}

$#inputfiles<2 or die "Can process either 1 file only!\n";

$filename=$inputfiles[0];
$count=0;
$totp=0;
$frames=0;
$nold=-99999;
$threshq=$thresh*$thresh;
open(INPUT, $filename) || die "Cannot open $filename\n";

while(<INPUT>){
    s/^\s*//;	# remove leading spaces
    @li = split('\s+');
    if($#li==0 && $li[0]+0==$li[0] && $count!=1){
	$n=$li[0]+0;
	if($n==$nold){
	    print;
        }elsif($nold!=-99999){
		print STDERR  "ERROR xyz_displacement: file $filename has frames with different numbers of particles: $n, $nold\n";
        }
	$frames++;
	$count=0;
    }elsif($count==1){
	$t=1. *$li[1];
	if($t==$li[1]){
	    $tnow=$t;
	    if($n==$nold){
		if($thresh>0){
		    print "#intermediate_time= ",0.5*($tnow+$told),"\n";
		}else{
		    print "#deltat= ",$tnow-$told,"\n";
		}
	    }
	}else{
	    if($n==$nold){print;}
	}
	$told=$tnow;
    }else{
	if($n==$nold){
	    $dx=$li[1]-$oldx[$count];
	    $dy=$li[2]-$oldy[$count];
	    $dz=$li[3]-$oldz[$count];
	    if($thresh>0){
		$delta=$dx*$dx+$dy*$dy+$dz*$dz;
		if($delta>$threshq){
		    print "Rb	$li[1]	$li[2]	$li[3]\n";
		}else{
		    print;
		}
	    }else{
		print "$li[0]	 $dx	$dy	$dz\n";
	    }
	}
	$oldx[$count]=$li[1];
	$oldy[$count]=$li[2];
	$oldz[$count]=$li[3];
	if($count==$n+1){
	    $nold=$n;
	}
    }
    $count++;
}
close(INPUT);