#!/usr/bin/perl  -I /home/manini/bin/perllib
# the above include is the location for math_module.pl
# couln't find a better way, e.g. using $HOME: does not work!

use POSIX;
use math_module;
# get it from http://www.lwebzem.com/cgi-bin/res/view.cgi?name=math_module.pm
# put math_module.pm in the folder named after -I in the first line above

sub print_help{
    print STDERR <<EOF;

Syntax:  xyz_in_cell [options] input_xyzfile

    generate an xyz file with particles translated back into one periodic cell

Options:
    -a1 a1x a1y a1z	the Cartesian components of the primitive vectors
    -a2 a2x a2y a2z
    -a3 a3x a3y a3z
    -ws			translate particles into the Wigner-Seitz cell
			(watch out, this is 9x slower!)
							by Nick 6 Sep 2017
EOF
    return;
}

$HUGE=1.2313e200;
# default: simple orthorombic cell:
$a=20.52;
$b=19.71;
$c=100.;
@a1=($a,0,0);
@a2=(0,$b,0);
@a3=(0,0,$c);
$ws=0; # by default report points in the CONVENTIONAL cell, not Wigner-Seitz

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 "-a1"){
	$x=shift(@ARGV)+0.; $y=shift(@ARGV)+0.; $z=shift(@ARGV)+0.;
	@a1=($x,$y,$z);
	next;
    }
    if($p eq "-a2"){
	$x=shift(@ARGV)+0.; $y=shift(@ARGV)+0.; $z=shift(@ARGV)+0.;
	@a2=($x,$y,$z);
	next;
    }
    if($p eq "-a3"){
	$x=shift(@ARGV)+0.; $y=shift(@ARGV)+0.; $z=shift(@ARGV)+0.;
	@a3=($x,$y,$z);
	next;
    }
    if($p eq "-ws"){
	$ws=1;
	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
}

# generate the transformation matrices:
for($i=0;$i<3;$i++){
    @mat[$i]=[$a1[$i],$a2[$i],$a3[$i]];
}    
@inv=math_module::get_inverse(3,@mat);

while ($filename=shift(@inputfiles)){
    $natom=0;
    $count=0;
    $totp=0;
    $frames=0;
    open(INPUT, $filename) || die "Cannot open $filename\n";
    
    while(<INPUT>){
	s/^\s*//;	# remove leading spaces
	@li = split('\s+');
	chomp();
	if($#li==0 && $count!=1){
	    $count=0;
	    print $_,"\n";
	}
	if($count==1){ # comment line
	    print $_,"refolded by xyz_in_cell\n";
	}
	if($count>1){
	    print $li[0];
	    for($i=0;$i<3;$i++){
		$r[$i]=$li[$i+1];
	    }
	    
	    for($i=0;$i<3;$i++){ # matrix multiplication by @inv
		$p[$i]=0;
		for($j=0;$j<3;$j++){
		    $p[$i]+=$inv[$i][$j]*$r[$j];
		}
	    }

	    for($i=0;$i<3;$i++){ # reduction to -0.5 to 0.5 interval
		$appo[$i]=$p[$i]-POSIX::floor($p[$i]+0.5);
	    }

	    @rk=(0.,0.,0.);
	    $mindist=$HUGE;
	    for($n1=-$ws;$n1<=$ws;$n1++){     # explore translations to nearest
		for($n2=-$ws;$n2<=$ws;$n2++){ # neighbors for the WS case
		    for($n3=-$ws;$n3<=$ws;$n3++){
			@nvec=($n1,$n2,$n3);
			for($i=0;$i<3;$i++){  # matrix multiplication by @mat:
			    $rr[$i]=0;
			    for($j=0;$j<3;$j++){
				$rr[$i]+=$mat[$i][$j]*($appo[$j]+$nvec[$j]);
			    }
			}
			$norm=$rr[0]*$rr[0]+$rr[1]*$rr[1]+$rr[2]*$rr[2];
			if($norm<$mindist){ # select shortest
			    $mindist=$norm;
			    @rk=@rr;
			}
		    }
		}
	    }
	    for($i=0;$i<3;$i++){
		print "\t",$rk[$i];
	    }
	    print "\n";
	}
	$count++;
    }
    close(INPUT);
}