package BandStructure;
use Exporter;
use Common;
@ISA = qw(Exporter Common);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
#use Crystal::MyUtility;

BEGIN {
}

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub IsParallel
{
	my ($this, $vx0, $vy0, $vz0, $vx1, $vy1, $vz1) = @_;
	return 0 if(!defined $vz0 or !defined $vz1);
	
	my $r0  = sqrt($vx0*$vx0 + $vy0*$vy0 + $vz0*$vz0);
	my $r1  = sqrt($vx1*$vx1 + $vy1*$vy1 + $vz1*$vz1);
	return 1 if($r0 < 1.0e-3 or $r1 < 1.0e-3);
	my $k = 1.0 / $r0 / $r1;

	my $opx = $vy0 * $vz1 - $vy1 * $vz0;
	my $opy = $vz0 * $vx1 - $vz1 * $vx0;
	my $opz = $vx0 * $vy1 - $vx1 * $vy0;
#print "opy = $opy = $vz0 * $vx1 - $vz1 * $vx0\n";
#print "op=($opx,$opy,$opz)\n";
	return 1 if(abs($opx * $k) < 1.0e-3 and abs($opy * $k) < 1.0e-3 and abs($opz * $k) < 1.0e-3);
	return 0;
}

sub DivideKList
{
	my ($this, $Crystal, $pklist, $nKPointDiv) = @_;

	my @dk;
	my $dtot = 0.0;
	for(my $i = 0 ; $i < @$pklist-1 ; $i++) {
		my $p0 = $pklist->[$i];
		my $p1 = $pklist->[$i+1];
		my $d = $Crystal->GetReciprocalDistanceFromK($p0->{kx}, $p0->{ky}, $p0->{kz}, $p1->{kx}, $p1->{ky}, $p1->{kz});
#print "i=$i: dk=$d dtot=$dtot\n";
		$dtot += $d;
		$dk[$i] = $d;
	}
	
	my $delta = $dtot / $nKPointDiv;
#print "delta = $delta\n";
	my @nk;
	my @klist;
	for(my $i = 0 ; $i < @$pklist-1 ; $i++) {
		my $n = $dk[$i] / $delta;
		if($n < 4.5) {
			$nk[$i] = 4;
		}
		elsif($n <= 6.5) {
			$nk[$i] = 4;
		}
		elsif($n <= 9.0) {
			$nk[$i] = 8;
		}
		elsif($n <= 12.0) {
			$nk[$i] = 10;
		}
		elsif($n <= 14.0) {
			$nk[$i] = 12;
		}
		elsif($n <= 16.0) {
			$nk[$i] = 15;
		}
		elsif($n <= 18.0) {
			$nk[$i] = 16;
		}
		elsif($n <= 22.0) {
			$nk[$i] = 20;
		}
		elsif($n <= 28.0) {
			$nk[$i] = 25;
		}
		else {
			$nk[$i] = int($n / 5.0) * 5;
		}

		my $dkx = ($pklist->[$i+1]->{kx} - $pklist->[$i]->{kx}) / $nk[$i];
		my $dky = ($pklist->[$i+1]->{ky} - $pklist->[$i]->{ky}) / $nk[$i];
		my $dkz = ($pklist->[$i+1]->{kz} - $pklist->[$i]->{kz}) / $nk[$i];
		for(my $j = 0 ; $j < $nk[$i] ; $j++) {
			my $kx = $pklist->[$i]->{kx} + $j*$dkx;
			my $ky = $pklist->[$i]->{ky} + $j*$dky;
			my $kz = $pklist->[$i]->{kz} + $j*$dkz;
			my $name = ($j == 0)? $pklist->[$i]->{name} : '';
			push(@klist, {name => $name, kx => $kx, ky => $ky, kz => $kz});
		}
#print "i=$i: dk=$dk[$i] / $delta = $n => $nk[$i]\n";
	}
	push(@klist, $pklist->[@$pklist-1]);
#print "ntot=", scalar @klist, "\n";

	return (scalar @klist, \@klist);
}

sub ReduceKList
{
	my ($this, $plist) = @_;

	my @klist = ($plist->[0]);

	my $vx0 = $plist->[1]{kx} - $plist->[0]{kx};
	my $vy0 = $plist->[1]{ky} - $plist->[0]{ky};
	my $vz0 = $plist->[1]{kz} - $plist->[0]{kz};
#print "  : 1: ($vx0, $vy0, $vz0)\n";
	for(my $i = 2 ; $i < @$plist ; $i++) {
		my $vx = $plist->[$i]{kx} - $plist->[$i-1]{kx};
		my $vy = $plist->[$i]{ky} - $plist->[$i-1]{ky};
		my $vz = $plist->[$i]{kz} - $plist->[$i-1]{kz};

		if($this->IsParallel($vx0, $vy0, $vz0, $vx, $vy, $vz)) {
#print " p: $i: ($vx, $vy, $vz) // ($vx0, $vy0, $vz0)\n";
		}
		else {
#print "np: $i: ($vx, $vy, $vz) !// ($vx0, $vy0, $vz0)\n";
			push(@klist, $plist->[$i-1]);
		}

		$vx0 = $vx;
		$vy0 = $vy;
		$vz0 = $vz;
	}
	push(@klist, $plist->[@$plist-1]);

	return @klist;
}

sub CheckKPOINTSFileType
{
	my ($this, $klistFile) = @_;
	
	my $in = JFile->new($klistFile, 'r');
	return '' if(!$in);
	
	my $line = $in->ReadLine();
	return 'KPOINTS:Explicit' if($line =~ /Explicit k-points/i);

	$line = $in->ReadLine();
	$line = $in->ReadLine();
	$in->Close();

	return 'klist:XCrySDen'   if($line =~ /[+\-\d]+\s+[+\-\d]+\s+[+\-\d]+\s+[+\-\d]+\s+[+\-\d\.]/);
	return 'KPOINTS:LineMode' if($line =~ /Line-mode/i);
	return '';
}

sub ReadKListFromKListXCrySDen
{
	my ($this, $infile) = @_;

	my $in = JFile->new($infile, 'r');
	return '' if(!$in);
	
	my @klist;
	while(1) {
		my $line = $in->ReadLine();
#print "l: $line";
		last if(!defined $line);
		last if($line =~ /^\s*END/i);

		my $name = '';
		if($line =~ /^\s*([a-zA-Z]+)/) {
			$name = $1;
			$line =~ s/^\s*[a-zA-Z]+\s+//;
		}
		my ($kx, $ky, $kz, $ndiv, $w) = Utils::Split("\\s+", $line);
		if($ndiv ne '') {
			push(@klist, {
					name   => $name,
					kx     => $kx/$ndiv,
					ky     => $ky/$ndiv,
					kz     => $kz/$ndiv,
					weight => $w,
					ikx    => $kx,
					iky    => $ky,
					ikz    => $kz,
					ndiv   => $ndiv,
				});
		}
	}
	$in->Close();

	return @klist;
}

sub ReadKListFromKPOINTSExplicit
{
	my ($this, $infile) = @_;

	my $in = JFile->new($infile, 'r');
	return '' if(!$in);

	my $line = $in->ReadLine();
	return '' if($line !~ /Explicit k-points/i);

	my $nk = $in->ReadLine() + 0;
#print "nk=$nk\n";
#Reciprocal
	$in->ReadLine();

	my @klist;
	my ($kx0, $ky0, $kz0);
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);

		my ($kx, $ky, $kz, $w) = Utils::Split("\\s+", $line);
#print "(k)=($kx, $ky, $kz)\n";
		next if(!defined $kz);

		push(@klist, {
#				name   => $name,
				kx     => $kx,
				ky     => $ky,
				kz     => $kz,
				weight => $w,
#				ikx    => undef,
#				iky    => undef,
#				ikz    => undef,
#				ndiv   => undef,
			});
	}
	$in->Close();

	return @klist;
}

sub ReadKListFromKPOINTSLineMode
{
	my ($this, $infile) = @_;

	my $in = JFile->new($infile, 'r');
	return '' if(!$in);

	my $line = $in->ReadLine();
	$line = $in->ReadLine();
	$line = $in->ReadLine();
	return '' if($line !~ /Line-mode/i);

	$line = $in->ReadLine();

	my @klist;
	my ($kx0, $ky0, $kz0);
	while(1) {
		$line = $in->ReadLine();
#print "l: $line";
		last if(!defined $line);

		Utils::DelSpace($line);
		next if($line eq '');

		my ($kx, $ky, $kz, $comment, $name) = Utils::Split("\\s+", $line);
		next if(!defined $kz);
		next if(abs($kx-$kx0) < 1.0e-4 and abs($ky-$ky0) < 1.0e-4 and abs($kz-$kz0) < 1.0e-4);

#print "GLH: ", $GLabelHash{a}, "\n";
		if(!defined $name) {
			if($comment) {
				$name = $comment;
			}
			else {
				$name = '';
			};
		}
#print "l:[$name, $kx, $ky, $kz]\n";

		push(@klist, {
				name   => $name,
				kx     => $kx,
				ky     => $ky,
				kz     => $kz,
				weight => undef,
				ikx    => undef,
				iky    => undef,
				ikz    => undef,
				ndiv   => undef,
			});
		$kx0 = $kx;
		$ky0 = $ky;
		$kz0 = $kz;
	}
	$in->Close();

	return @klist;
}

sub ReadKListFile
{
	my ($this, $klistFile, $ReduceKPoints) = @_;
	$ReduceKPoints = 1 if(!defined $ReduceKPoints);

	my $type = $this->CheckKPOINTSFileType($klistFile);
print "File type of [$klistFile]: $type\n";
	my @klist;
	if($type eq 'KPOINTS:LineMode') {
		@klist = $this->ReadKListFromKPOINTSLineMode($klistFile);
	}
	elsif($type eq 'KPOINTS:Explicit') {
		@klist = $this->ReadKListFromKPOINTSExplicit($klistFile);
	}
	elsif($type eq 'klist:XCrySDen') {
		@klist = $this->ReadKListFromKListXCrySDen($klistFile);
	}

#for(my $i = 0 ; $i < @klist ; $i++) {
#	my $p = $klist[$i];
#	print "$i: $p->{name}: k=($p->{kx}, $p->{ky}, $p->{kz}) w=$p->{weight}\n";
#}

	if($ReduceKPoints) {
		@klist = $this->ReduceKList(\@klist);
	}

for(my $i = 0 ; $i < @klist ; $i++) {
	my $p = $klist[$i];
	print "$i: $p->{name}: k=($p->{kx}, $p->{ky}, $p->{kz}) w=$p->{weight}\n";
}

	return @klist;
}

sub AnalyzeaKProduct
{
	my ($this, $Crystal, $aKProduct) = @_;

	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);

	my ($nx, $ny, $nz);
	my ($NKREDX, $NKREDY, $NKREDZ) = ();
	if($aKProduct =~ /^fix:(.*)$/i) {
		my ($a1, $a2, $a3) = Utils::Split(",", $1);
		($nx, $NKREDX) = Utils::Split("/", $a1);
		($ny, $NKREDY) = Utils::Split("/", $a2);
		($nz, $NKREDZ) = Utils::Split("/", $a3);
	}
	else {
		my $an = $aKProduct * 10.0;
		$nx = int($an / $a) + 1;
		$ny = int($an / $b) + 1;
		$nz = int($an / $c) + 1;
	}
	return ($nx, $ny, $nz, $NKREDX, $NKREDY, $NKREDZ);
}

1;