#!/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_color_density [options] xyzfile

generate another xyz file with chemically changed atoms at higher local density
      NBB!! requires nearestneigbors.x to be installed somewhere in the path!!
      (contact the author for a copy of the source for nearestneigbors)

Options:
	-f1 %threshold  : set the first fractional threshold relative to the average d below which changes the atom kind (default=$fraction1)
	-f2 %threshold  : set the second fractional threshold relative to the average d above which changes the atom kind (default=$fraction2)
	-t1 threshold  : set the first threshold d below which changes the atom kind (default=$thresh1)
	-t2 threshold  : set the second threshold d above which changes the atom kind (default=$thresh2)
							by Nick v. 01 Aug 2013
EOF
    return;
}

$PI=3.14159265358979324;
$BIG=99.e99;
$a=1;
$fraction1=0.8;
$fraction2=1.2;
$thresh1="unset";
$thresh2="unset";

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 "-f1"){ # help
	$fraction1=shift(@ARGV)+0;
	$fraction1>0 or die "xyz_color_density: valid %threshold had better be positive\n";
	next;
    }
    if($p eq "-f2"){ # help
	$fraction2=shift(@ARGV)+0;
	$fraction2>0 or die "xyz_color_density: valid %threshold had better be positive\n";
	next;
     }
    if($p eq "-t1"){ # help
	$thresh1=shift(@ARGV)+0;
	$thresh1>0 or die "xyz_color_density: valid threshold had better be positive\n";
	next;
    }
    if($p eq "-t2"){ # help
	$thresh2=shift(@ARGV)+0;
	$thresh2>0 or die "xyz_color_density: valid threshold had better be positive\n";
	next;
    }
    if($p eq "-h"){ # help
	print_help;
	die "\n";
    }
    print_help;
    die "invalid option $p\n";
}

if($thresh1 ne "unset"){
    $fraction1="unset";
}
if($thresh2 ne "unset"){
    $fraction2="unset";
}

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

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

$filename=$inputfiles[0];
$count=0;
$totp=0;
$frames=0;

$tmpfile="/tmp/xyzcd.$$";

`nearestneigbors.x < $filename > $tmpfile`;

open(INPUT, $filename) || die "Cannot open $filename\n";
open(INNEI, $tmpfile) || die "Cannot open $tmpfile\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($#a>0){
	    &readandprintout;
	    @a=();
	    @x=();
	    @y=();
	    @z=();
	}
	$nold=$n;
	print $n,"\n";
	$frames++;
	$count=0;
    }elsif($count==1){
	print;
    }else{
	push @a,$li[0];
	push @x,$li[1];
	push @y,$li[2];
	push @z,$li[3];
    }
    $count++;
}
close(INPUT);
&readandprintout;
close(INNEI);
`rm $tmpfile`;


sub readandprintout{
    $totaldave=0;
    $totaln=0;
    $_=<INNEI>;
    $_=<INNEI>;
    for($i=1;$i<=$nold;$i++){
	$_=<INNEI>;
	s/^\s*//;	# remove leading spaces
	@ln = split('\s+');
	$i==$ln[0] or die "nonmatching in xyz_color_density $i $ln[0]\n";
	$dave=0;
	for($neig=1;$neig<=$ln[1];$neig++){
	    $dist=sqrt(($x[$i-1]-$x[$ln[$neig+1]-1])**2+
		    ($y[$i-1]-$y[$ln[$neig+1]-1])**2+
		    ($z[$i-1]-$z[$ln[$neig+1]-1])**2);
#		   print "QUII $i    $ln[1]  $neig ".$ln[$neig+1]." $d\n";
#		   print "QUAA $z[$i] ".$z[$ln[$neig+1]]."\n";
	    $dave+=$dist;
	}
	$totaldave+=$dave;
	$totaln+=$ln[1];
	if($ln[1]>0){
	    $dave/=$ln[1];
	}
	push @d,$dave;
    }
    $totaldave/=$totaln;

    if($fraction1 ne "unset"){
	$thresh1=$totaldave*$fraction1;
    }
    if($fraction2 ne "unset"){
	$thresh2=$totaldave*$fraction2;
    }

    print STDERR "frame ", $frames," average blength= ",$totaldave," over ", $totaln, " bonds\n";

    for($i=0;$i<$nold;$i++){
	$dave=$d[$i];
	if($dave<$thresh1){
	    $aname="Li";
	}elsif($dave<$thresh2){
	    $aname=$a[$i];
	}else{
	    $aname="Na";
	}
	print $aname,"\t",$x[$i]," ",$y[$i]," ",$z[$i]," ",$dave,"\n";
    }
    @d=();
}