#!/usr/bin/perl -w
#
# k-point generator for energy band plot
# version 1.01a
#
# input file:
# distance
# b1x b2x b3x
# b1y b2y b3y
# b1z b2z b3z
# k1 k2 k3 denom # X
# k1 k2 k3 denom # Gamma
# k1 k2 k3 denom # L
# ....
#
# Copyright (C) 2003-2004, Takenori YAMAMOTO
# Institute of Industrial Science, The University of Tokyo
#

my @dummy;

if(@ARGV<1) {
print "Usage: band_kpoint.pl InputFile\n";
print <<end;
Format of input file:
distance
b1x b2x b3x
b1y b2y b3y
b1z b2z b3z
k1 k2 k3 denom # X
k1 k2 k3 denom # Gamma
k1 k2 k3 denom # L
....
end
exit;
}

open(IN,$ARGV[0]);
$distance = <IN>;
chomp($distance);
@{$rlvec[0]} = split(' ',<IN>);
@{$rlvec[1]} = split(' ',<IN>);
@{$rlvec[2]} = split(' ',<IN>);
$i=0;
while($line = <IN>) {
  ($line,@dummy) = split('#',$line);
  @{$skp[$i]} = split(' ',$line);
  die "bad input: i=$i skp=@{$skp[$i]}\n" if(@{$skp[$i]} != 4);
  for($j=0;$j<3;$j++) {
    $skv[$i][$j] = $skp[$i][$j]/$skp[$i][3];
  }
  $i++;
}
close(IN);
$nsk = @skv;

# Division numbers will be calculated by the inputed distance.
for($i=1;$i<$nsk;$i++) {
  @dkc=(0,0,0);
  for($j=0;$j<3;$j++) {
    for($k=0;$k<3;$k++) {
      $dkc[$j] += $rlvec[$j][$k]*($skv[$i][$k]-$skv[$i-1][$k]);
    }
  }
  $dst = sqrt($dkc[0]**2+$dkc[1]**2+$dkc[2]**2);
  print "Distance of $i = $dst\n";
  $ndiv[$i-1] = int($dst/$distance);
}

print "division numbers = @ndiv\n";

# orientation vectors
for($i=0;$i<$nsk-1;$i++) {
  for($j=0;$j<3;$j++) {
    $orivec[$i][$j] = ($skv[$i+1][$j]-$skv[$i][$j])/$ndiv[$i];
  }
  print "i=$i @{$orivec[$i]}\n";
}

# generate k-points
$nkv=0;
$deno[$nkv] = $skp[0][3]*$skp[1][3]*$ndiv[0];
for($k=0;$k<3;$k++) {
  $kpt[$k] = $skv[0][$k];
  if($kpt[$k] > 0) {
    $nume[$nkv][$k] = int($kpt[$k]*$deno[$nkv]+0.5);
  } else {
    $nume[$nkv][$k] = int($kpt[$k]*$deno[$nkv]-0.5);
  }
}
$n=$nkv+1;
print "$n : @kpt\n    => $nume[$nkv][0]/$deno[$nkv] $nume[$nkv][1]/$deno[$nkv] $nume[$nkv][2]/$deno[$nkv]\n";
$nkv++;
for($i=1;$i<$nsk;$i++) {
  @kpt = @{$skv[$i-1]};
  for($j=0;$j<$ndiv[$i-1];$j++) {
    $deno[$nkv] = $skp[$i-1][3]*$skp[$i][3]*$ndiv[$i-1];
    for($k=0;$k<3;$k++) {
      $kpt[$k] += $orivec[$i-1][$k];
      if($kpt[$k] > 0) {
        $nume[$nkv][$k] = int($kpt[$k]*$deno[$nkv]+0.5);
      } else {
        $nume[$nkv][$k] = int($kpt[$k]*$deno[$nkv]-0.5);
      }
    }
    $n=$nkv+1;
    print "$n : @kpt\n    => $nume[$nkv][0]/$deno[$nkv] $nume[$nkv][1]/$deno[$nkv] $nume[$nkv][2]/$deno[$nkv]\n";
    $nkv++;
  }
}

# write data in PHASE format
open(OUT,">kpoint.data");

print OUT "$nkv $nkv\n";
for($i=0;$i<$nkv;$i++) {
  print OUT "$nume[$i][0] $nume[$i][1] $nume[$i][2] $deno[$i] 1\n";
}

close(OUT);