#!/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_sphere_density [options] xyzfile1 [xyzfile2, ...]

Evaluate the average particles density in a sphere (circle) from a xyz file

Options:
	-r fix radius of the sphere (default: $radius)
	-x fix the x-coord of the center of the sphere (default: $xcube)
	-y fix the y-coord of the center of the sphere (default: $ycube)
	-z fix the z-coord of the center of the sphere (default: $zcube)
	-2 compute 2D density (based on a circle)
	-ti initial time	(assuming time is the second word in comment!)
	-tf final time  	(assuming time is the second word in comment!)
						version 26/06/2012 by Nick
EOF
    return;
}

$PI=3.14159265358979324;
$BIG=9.999e200;

# defaults:
$radius=10;
$xcenter=0;
$ycenter=0;
$zcenter=0;
$initialframe=0;
$timein=-$BIG;
$timefi=$BIG;
$timeinActual=-$BIG;
$timefiActual=$BIG;
$filename="-";
$twoD=0;

while($p=shift(@ARGV)){
    if($p=~/^[^\-]/){ # all command line options start with "-", otherwise it a filename
	$filename=$p;
	next;
    }
    if($p eq "-r"){
        $radius=shift(@ARGV)+0;
        $radius>0 or die "Invalid center radius $radius\n";
        next;
    }
    if($p eq "-x"){
        $xcenter=shift(@ARGV)+0;
        next;
    }
    if($p eq "-y"){
        $ycenter=shift(@ARGV)+0;
        next;
    }
    if($p eq "-z"){
        $zcenter=shift(@ARGV)+0;
        next;
    }
    if($p eq "-2"){
        $twoD=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";
}

$maxatoms=0;
$actualframes=0;
$deltaradius=0.1;
$nradii=19;

for($ir=0;$ir<$nradii;$ir++){
    $radius2[$ir]=($radius+($ir-$nradii/2)*$deltaradius)**2;
}
for($ir=0;$ir<$nradii;$ir++){
    $countin[$ir]=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>;
	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++;
	    for($i=0;$i<$natms;$i++){
		$_=<INPUT>;
		s/^\s*//; @li = split('\s+');
		$name=$li[0]; # store everything
		$x=$li[1];
		$y=$li[2];
		$z=$li[3];

		$d2=($li[1]-$xcenter)*($li[1]-$xcenter)+
		    ($li[2]-$ycenter)*($li[2]-$ycenter)+
		    ($li[3]-$zcenter)*($li[3]-$zcenter);

		for($ir=0;$ir<$nradii;$ir++){
		    if($d2<$radius2[$ir]){
			$countin[$ir]++;
                    }
		}
	    }
	}
    }
}
close(INPUT);

$avedens=0;
$avedens2=0;
for($ir=0;$ir<$nradii;$ir++){
    if($twoD==0){
	$volume=4.*$PI/3. * $radius2[$ir]**1.5; # sphere
    }else{
	$volume=$PI * $radius2[$ir]; # circle, 2D case
    }
    $dens=$countin[$ir]/($volume*$actualframes);
    $avedens+=$dens;
    $avedens2+=$dens*$dens;
}
$avedens/=$nradii;
$avedens2/=$nradii;

print "xyz_sphere_density: average density: ",$avedens," +/- ",sqrt(($avedens2-$avedens*$avedens)/$nradii),"\n";

if($twoD==0){
    $d=1.0/sqrt(2.)*(4.0/$avedens)**(1.0/3.); # fcc lattice
}else{
    $d=1.0/sqrt(sqrt(3.)/2. * $avedens); # triangular lattice
}
print "xyz_sphere_density: average atom-atom distance: ",$d,"\n";