#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';
#use lib 'd:/MyWebs/cgi-bin/lib';

use lib '/home/tkamiya/bin/lib';
use lib '/home/tkamiya/bin';

use strict;
#use warnings;

use JFile;

my $InputMode      = 'Cartesian';
#my $InputMode     = 'Reciprocal';

my $OutFile = "FreeElectronBand.csv";
my $nk      = 11;
my $krange  = 1;

#FCC Cartesian
my @KList;
if($InputMode eq 'Cartesian') {
	@KList = (
		[0.5,  0.0,  1.0, 'W'],
		[0.5,  0.5,  0.5, 'L'],
		[0.0,  0.0,  0.0, 'GAMMA'],
		[0.0,  0.0,  1.0, 'X'],
		[0.5,  0.0,  1.0, 'W'],
		[0.75, 0.75, 0.0, 'K'],
		);
}
else {
@KList = (
		[0.0,  0.0,  0.0,  'GAMMA'],
		[0.5,  0.5,  0.0,  'X'],
		[0.5,  0.75, 1.25, 'W'],
		[0.5,  0.5,  0.5,  'L'],
#		[0.375, 0.75, 0.375,  'W'],
		);
}

my @k;
my ($kxp, $kyp, $kzp);
my $ktot = 0.0;
for(my $i = 0 ; $i < @KList-1 ; $i++) {
	my ($kx0, $ky0, $kz0, $name0) = @{$KList[$i]};
	my ($kx1, $ky1, $kz1, $name1) = @{$KList[$i+1]};
	my $dkx = ($kx1 - $kx0) / ($nk - 1);
	my $dky = ($ky1 - $ky0) / ($nk - 1);
	my $dkz = ($kz1 - $kz0) / ($nk - 1);
#print "k0=$kx0,$ky0,$kz0\n";
#print "k1=$kx1,$ky1,$kz1\n";
	if($i == 0) {
		($kxp, $kyp, $kzp) = ($kx0, $ky0, $kz0);
	}
	for(my $j = 0 ; $j < $nk ; $j++) {
		my $kx = $kx0 + $j * $dkx;
		my $ky = $ky0 + $j * $dky;
		my $kz = $kz0 + $j * $dkz;
		my $k2 = ($kx - $kxp)**2 + ($ky - $kyp)**2 + ($kz - $kzp)**2;
		$ktot += $k2;
		push(@k, [$kx, $ky, $kz, $ktot, "$name0-$name1"]);
		($kxp, $kyp, $kzp) = ($kx, $ky, $kz);
	}
}

my $out = JFile->new($OutFile, 'w') or die "$!: Can not write to [$OutFile].\n";
$out->print("name,kx,ky,kz,k");
for(my $ikz = -$krange ; $ikz <= $krange ; $ikz++) {
for(my $iky = -$krange ; $iky <= $krange ; $iky++) {
for(my $ikx = -$krange ; $ikx <= $krange ; $ikx++) {
	$out->print(",E($ikx $iky $ikz)");
}
}
}
$out->print("\n");

my @E;
for(my $ikz = -$krange ; $ikz <= $krange ; $ikz++) {
for(my $iky = -$krange ; $iky <= $krange ; $iky++) {
for(my $ikx = -$krange ; $ikx <= $krange ; $ikx++) {
	for(my $i = 0 ; $i < @k ; $i++) {
		my $p = $k[$i];
		my ($kx, $ky, $kz, $ktot, $name) = @$p;
		my $k2 = ($kx + $ikx)**2 + ($ky + $iky)**2 + ($kz + $ikz)**2;
		$E[$i][$ikx+$krange][$iky+$krange][$ikz+$krange] = $k2;
	}
}
}
}

for(my $i = 0 ; $i < @k ; $i++) {
	my $p  = $k[$i];
	$out->print("$p->[4],$p->[0],$p->[1],$p->[2],$p->[3]");
	for(my $ikz = -$krange ; $ikz <= $krange ; $ikz++) {
	for(my $iky = -$krange ; $iky <= $krange ; $iky++) {
	for(my $ikx = -$krange ; $ikx <= $krange ; $ikx++) {
		my $E = $E[$i][$ikx+$krange][$iky+$krange][$ikz+$krange];
		$out->print(",$E");
	}
	}
	}
	$out->print("\n");
}

$out->Close();
exit;