#!/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:  cut_cube_in_xyz [options] file.xyz

cuts a xyz file to include only atoms with a given cube
     and optionally a given time interval
options:
	-s fix side of the cube (default: $side)
	-lx fix x-side of the parallelepiped
	-ly fix y-side of the parallelepiped
	-lz fix z-side of the parallelepiped
	-x fix the x-coord of the center of the cube (default: $xcube)
	-y fix the y-coord of the center of the cube (default: $ycube)
	-z fix the z-coord of the center of the cube (default: $zcube)
	-i initial frame: keep atoms based on positions in first frame only
	-ti initial time	(assuming time is the second word in comment!)
	-tf final time  	(assuming time is the second word in comment!)
						version 29/08/2017 by Nick
EOF
    return;
}

$BIG=9.999e200;

# defaults:
$side=10;
$sidex=$side;
$sidey=$side;
$sidez=$side;
$xcube=0;
$ycube=0;
$zcube=0;
$initialframe=0;
$timein=-$BIG;
$timefi=$BIG;
$timeinActual=-$BIG;
$timefiActual=$BIG;
$filename="-";

while($p=shift(@ARGV)){
    if($p=~/^[^\-]/){ # all command line options start with "-", otherwise it a filename
	$filename=$p;
	next;
    }
    if($p eq "-s"){
        $side=shift(@ARGV)+0;
        $side>0 or die "Invalid cube side $side\n";
	$sidex=$side;
	$sidey=$side;
	$sidez=$side;
        next;
    }
    if($p eq "-lx"){
        $sidex=shift(@ARGV)+0;
        $sidex>0 or die "Invalid cube side $sidex\n";
        next;
    }
    if($p eq "-ly"){
        $sidey=shift(@ARGV)+0;
        $sidey>0 or die "Invalid cube side $sidey\n";
        next;
    }
    if($p eq "-lz"){
        $sidez=shift(@ARGV)+0;
        $sidez>0 or die "Invalid cube side $sidez\n";
        next;
    }
    if($p eq "-x"){
        $xcube=shift(@ARGV)+0;
        next;
    }
    if($p eq "-y"){
        $ycube=shift(@ARGV)+0;
        next;
    }
    if($p eq "-z"){
        $zcube=shift(@ARGV)+0;
        next;
    }
    if($p eq "-i"){
        $initialframe=1;
        next;
    }
    if($p eq "-ti"){
        $timein=shift(@ARGV)+0.;
        next;
    }
    if($p eq "-tf"){
        $timefi=shift(@ARGV)+0.;
        next;
    }
    if($p eq "-h"){ # help
	print_help;
	die "\n";
    }
    print_help;
    die "invalid option $p\n";
}

$sidex*=0.5;
$sidey*=0.5;
$sidez*=0.5;

$maxatoms=0;
$actualframes=0;

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]){
	$natms = $li[0];
	$_=<INPUT>;
	$commentline=$_;
	s/^\s*//; @li = split('\s+');
	$time=$li[1]+0.;
	if($timein <= $time && $time <= $timefi){
	    if($timeinActual==-$BIG){
		$timeinActual=$time;
	    }
	    $timefiActual=$time; # eventually it remains equal to the last processed time
	    $actualframes++;
#	    @nli=(); @xli=(); @yli=(); @zli=();
	    @eli=();
	    if($initialframe){
		for($i=0;$i<$natms;$i++){
		    $_=<INPUT>;
		    s/^\s*//; @li = split('\s+');
		    push @eli,$_;
		    if($actualframes==1 && 
		       abs($li[1]-$xcube)<$sidex &&
		       abs($li[2]-$ycube)<$sidey &&
		       abs($li[3]-$zcube)<$sidez ){
			push @indexes,$i
		    }
		}
	    }else{
		for($i=0;$i<$natms;$i++){
		    $_=<INPUT>;
		    s/^\s*//; @li = split('\s+');
		    if(abs($li[1]-$xcube)<$sidex &&
		       abs($li[2]-$ycube)<$sidey &&
		       abs($li[3]-$zcube)<$sidez ){
			push @eli,$_
		    }
		}
	    }
	    if($initialframe){
		$nat=$#indexes+1;
	    }else{
		$nat=$#eli+1;
	    }
	    if($nat>$maxatoms){
		$maxatoms=$nat;
	    }
	    $natms>=$nat or die "cut_cube_in_xyz: input file $filename incompatible with option -i from frame $actualframes\n";
	    print $nat,"\n";	# number of atoms in cube
	    print $commentline;
	    if($initialframe){
		for($j=0;$j<$nat;$j++){
		    $i=$indexes[$j];
		    print $eli[$i];
		}
	    }else{
		for($i=0;$i<=$#eli;$i++){
		    print $eli[$i];
		}
	    }
	}else{ # the following is not really needed, but makes skips quicker:
	    for($i=0;$i<$natms;$i++){<INPUT>;};
	}
    }
}
close(INPUT);
print STDERR "cut_cube_in_xyz: selected ",$actualframes," frames";
if($timeinActual!=-$BIG){
    print STDERR ", from t=", $timeinActual;
}
if($timefi!=$BIG){
    print STDERR ", until t=",$timefiActual;
}
print STDERR ", with ";
if(!$initialframe){
    print STDERR "at most ";
}
print STDERR $maxatoms," atoms per frame";
if($initialframe && $actualframes>1){
    print STDERR ", based on their positions in the first frame";
}
print STDERR ".\n";