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

cuts a xyz file to include only atoms with a given sphere
     and optionally a given time interval
options:
	-r fix the radius of the sphere (default: $r)
	-x fix the x-coord of the center of the sphere (default: $xsphere)
	-y fix the y-coord of the center of the sphere (default: $ysphere)
	-z fix the z-coord of the center of the sphere (default: $zsphere)
	-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 25/08/2018 by Nick
EOF
    return;
}

$BIG=9.999e200;

# defaults:
$r=10;
$xsphere=0;
$ysphere=0;
$zsphere=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 "-r"){
        $r=shift(@ARGV)+0;
        $r>0 or die "Invalid cube side $side\n";
	next;
    }
    if($p eq "-x"){
        $xsphere=shift(@ARGV)+0;
        next;
    }
    if($p eq "-y"){
        $ysphere=shift(@ARGV)+0;
        next;
    }
    if($p eq "-z"){
        $zsphere=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";
}

$r2=$r*$r;
$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=();
	    if($initialframe){
		for($i=0;$i<$natms;$i++){
		    $_=<INPUT>;
		    s/^\s*//; @li = split('\s+');
		    push @nli,$li[0]; # store everything
		    push @xli,$li[1];
		    push @yli,$li[2];
		    push @zli,$li[3];
		    if($actualframes==1 && 
		        abs($li[1]-$xsphere)**2 +
			abs($li[2]-$ysphere)**2 +
			abs($li[3]-$zsphere)**2 <$r2 ){
			push @indexes,$i
		    }
		}
	    }else{
		for($i=0;$i<$natms;$i++){
		    $_=<INPUT>;
		    s/^\s*//; @li = split('\s+');
		    if( abs($li[1]-$xsphere)**2 +
			abs($li[2]-$ysphere)**2 +
			abs($li[3]-$zsphere)**2 <$r2 ){
			push @nli,$li[0]; # store only atoms in sphere
			push @xli,$li[1];
			push @yli,$li[2];
			push @zli,$li[3];
		    }
		}
	    }
	    if($initialframe){
		$nat=$#indexes+1;
	    }else{
		$nat=$#nli+1;
	    }
	    if($nat>$maxatoms){
		$maxatoms=$nat;
	    }
	    $natms>=$nat or die "cut_sphere_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 $nli[$i],"\t",$xli[$i],"\t",$yli[$i],"\t",$zli[$i],"\n";
		}
	    }else{
		for($i=0;$i<=$#nli;$i++){
		    print $nli[$i],"\t",$xli[$i],"\t",$yli[$i],"\t",$zli[$i],"\n";
		}
	    }
	}else{ # the following is not really needed, but makes skips quicker:
	    for($i=0;$i<$natms;$i++){<INPUT>;};
	}
    }
}
close(INPUT);
print STDERR "cut_sphere_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";