#!/usr/bin/perl

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
#print "\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "d:/Programs/Perl/lib", @INC);
}

use strict;

use MyHTMLApplication;
use Utils;
use Sci;
use JFile;
use Sci::Optics;
use Crystal::VASP;
use Crystal::XCrySDen;

my $App = new MyHTMLApplication;

my ($EStart, $EEnd) = (0.0, 10.0); # eV
my $OUTCAR        = "OUTCAR";
my $OpticsOut     = "Optics.csv";
my $VibrationOut  = "Vibration.txt";
my $XSFFile       = "a.xsf";
my $kDisplacement = 0.1;
$OUTCAR       = $ARGV[0] if(defined $ARGV[0]);
$OpticsOut    = $ARGV[1] if(defined $ARGV[1]);
$VibrationOut = $ARGV[2] if(defined $ARGV[2]);
$XSFFile      = $ARGV[3] if(defined $ARGV[3]);

&SaveOpticsCSV($App, $OUTCAR, $OpticsOut, $EStart, $EEnd);
&SaveVibrationFile($App, $OUTCAR, $VibrationOut);
exit;

sub SaveVibrationFile
{
	my ($App, $OUTCAR, $VibrationOut) = @_;

	$App->print("Save Vibration data from [$OUTCAR] to [$VibrationOut].\n");
	my $in = new JFile;
	$in->Open($OUTCAR, "r");
	if(!$in) {
		$App->print("Error: Can not read [$OUTCAR].\n");
		return;
	}
	my $line = $in->SkipTo("^\\s*Eigenvectors and eigenvalues of the dynamical");
	if(!defined $line) {
		$App->print("Error: can not find Vibration data\n");
		return;
	}

	my $out = new JFile();
	$out->Open($VibrationOut, "w");
	if(!$out) {
		$App->print("Error: Can not write to [$VibrationOut].\n");
		return;
	}
	$out->print($line);

	my @Vibration;
	my $count = 0;
	while(1) {
		$line = $in->ReadLine();
		last if($line =~ /^-----/);
		if($line =~ /\d+\s+f\s*=\s*([\+\-\d\.eE]+)/) {
			my $THz = $1;
			$in->ReadLine();
			my @lines;
			my $c = 0;
			while(1) {
				$line = $in->ReadLine();
				Utils::DelSpace($line);
				last if($line eq '');
				$lines[$c] = $line;
#$App->print("   $line\n");
				$c++;
			}
			$Vibration[$count] = [$THz, \@lines];
#			$Vibration[$count] = {THz => $THz, pLines => \@lines};
			$count++;
		}
	}
	@Vibration = sort { $a->[0] <=> $b->[0]; } @Vibration;
#	@Vibration = sort { $a->[THz} <=> $b->{THz}; } @Vibration;

	for(my $i = 0 ; $i < @Vibration ; $i++) {
		my $p = $Vibration[$i];
		my $THz = $p->[0];
		my $eV  = Optics::HzToeV($THz * 1.0e12);
		my $kiser = Optics::eVToKiser($eV);
		$App->print("$i: $THz THz   $eV eV   $kiser cm-1\n");
		$out->printf("%3d: %8.3f THz   %8.4g eV   %8.3f cm-1\n", $i, $THz, $eV, $kiser);
	}
	$out->print("\n");

	$App->print("             X         Y         Z           dx          dy          dz\n");
	$out->print("             X         Y         Z           dx          dy          dz\n");
	for(my $i = 0 ; $i < @Vibration ; $i++) {
		my $p = $Vibration[$i];
		my $THz = $p->[0];
		my $pl  = $p->[1];
		my $eV  = Optics::HzToeV($THz * 1.0e12);
		my $kiser = Optics::eVToKiser($eV);
		$App->print("$i: $THz THz   $eV eV   $kiser cm-1\n");
		$out->print("$i: $THz THz   $eV eV   $kiser cm-1\n");
		for(my $j = 0 ; $j < @$pl ; $j++) {
			$App->print("   $pl->[$j]\n");
			$out->print("   $pl->[$j]\n");
		}
	}

	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);
		last if($line =~ /^ writing wavefunctions/);
		$out->print($line);
#		last if($line =~ /^  z  / or $line =~ /^ ZX /);
	}

	$out->Close();
	$in->Close();
	$App->print("Vibration data saved to [$VibrationOut].\n");

	my $vasp = new VASP;
	my $xcd = new XCrySDen();
	my $pCrystal = $vasp->ReadStructureFromCARFiles($OUTCAR, 1);

	my ($a, $b, $c, $alpha, $beta, $gamma) = $pCrystal->LatticeParameters();
	my ($a11,$a12,$a13,$a21,$a22,$a23,$a31,$a32,$a33) = $pCrystal->LatticeVectors();
	my @AtomType = $pCrystal->GetCAtomTypeList();
#	$pCrystal->ExpandCoordinates();
	my @AtomList = $pCrystal->GetCExpandedAtomSiteList();
	my $nAtom = @AtomList;
	my $nStep = scalar @Vibration;
	if(!$xcd->WriteXSFAnimationFileHeader($XSFFile, $nStep)) {
		$App->print("Error: can not write to [$XSFFile].\n");
		return;
	}

	unless(open(OUT,">>$XSFFile")) {
		return 0;
	}
	binmode(OUT);
#$App->print("n=$nAtom, $nStep\n");
	for(my $i = 0 ; $i < $nStep ; $i++) {
		my $p = $Vibration[$i];
		my $pl  = $p->[1];

		my $i1 = $i+1;
		print  OUT "PRIMVEC $i1\n";
		printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
		printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
		printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
		print  OUT "CONVVEC $i1\n";
		printf OUT " %16.10f %16.10f %16.10f\n", $a11,$a12,$a13;
		printf OUT " %16.10f %16.10f %16.10f\n", $a21,$a22,$a23;
		printf OUT " %16.10f %16.10f %16.10f\n", $a31,$a32,$a33;
		print  OUT "PRIMCOORD $i1\n";
		printf OUT " %d  1\n", $nAtom;
		for(my $j = 0 ; $j < @$pl ; $j++) {
			my ($xc, $yc, $zc, $vx, $vy, $vz) = Utils::Split("\\s+", $pl->[$j]);
#			my ($xc,$yc,$zc) = $pCrystal->FractionalToCartesian($x,$y,$z);
			my $atom         = $AtomList[$j];
			my $name         = $atom->AtomNameOnly();
#			my $AtomicNumber = $atom->AtomicNumber();
			my $AtomicNumber = AtomType::GetAtomicNumber($name);
			printf OUT "%3d   %12.7f %12.7f %12.7f   %12.7f %12.7f %12.7f\n",
				$AtomicNumber, $xc, $yc, $zc, $vx*$kDisplacement, $vy*$kDisplacement, $vz*$kDisplacement;
		}
	}
	close(OUT);
}

sub SaveOpticsCSV
{
	my ($App, $OUTCAR, $OpticsOut, $EStart, $EEnd) = @_;

	$App->print("Save OPTICS data from [$OUTCAR] to [$OpticsOut].\n");

	my @pData;
	for(my $i = 0 ; $i < 13 ; $i++) {
		$pData[$i] = [];
	}

	my $in = new JFile;
	$in->Open($OUTCAR, "r");
	if(!$in) {
		$App->print("Error: Can not read [$OUTCAR].\n");
		return;
	}
	my $line = $in->SkipTo("^\\s*frequency dependent IMAGINARY DIELECTRIC");
	if(!defined $line) {
		$App->print("Error: can not find IMAGINARY DIELECTRIC Function\n");
		return;
	}
#$App->print($line);
	$in->ReadLine();
	$in->ReadLine();
	my $count = 0;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);

#		my @a = Utils::Split("\\s+", $line);
		my @a = ($line =~ /^(.{12})(.{12})(.{12})(.{12})(.{12})(.{12})(.{12})/);
		last if(@a < 7);
		next if($a[0] < $EStart or $EEnd < $a[0]);
		for(my $i = 0 ; $i < 7 ; $i++) {
#Utils::DelSpace($a[$i]);
			$pData[$i]->[$count] = $a[$i] + 0.0;
		}
#$App->print("$count: $pData[0]->[$count], $pData[1]->[$count], $pData[3]->[$count]\n");
		$count++;
	}

	my $line = $in->SkipTo("^\\s*frequency dependent      REAL DIELECTRIC");
	if(!defined $line) {
		$App->print("Error: can not find REAL DIELECTRIC Function\n");
		return;
	}
#$App->print($line);
	$in->ReadLine();
	$in->ReadLine();
	$count = 0;
	while(1) {
		$line = $in->ReadLine();
		last if(!defined $line);

#		my @a = Utils::Split("\\s+", $line);
		my @a = ($line =~ /^(.{12})(.{12})(.{12})(.{12})(.{12})(.{12})(.{12})/);
		last if(@a < 7);
		next if($a[0] < $EStart or $EEnd < $a[0]);
		if($a[0]+0.0 ne $pData[0]->[$count]) {
			$App->print("Error: Energy values mismatch for E(real)=$a[0] and E(imag)=$pData[0]->[$count]\n");
			return;
		}
		for(my $i = 7 ; $i < 13 ; $i++) {
#Utils::DelSpace($a[$i-7]);
			$pData[$i]->[$count] = $a[$i-6] + 0.0;
		}
#$App->print("$count: $pData[0]->[$count], $pData[7]->[$count], $pData[9]->[$count]\n");
		$count++;
	}

	$in->Close();

	my $out = new JFile();
	$out->Open($OpticsOut, "w");
	if(!$out) {
		$App->print("Error: Can not write to [$OpticsOut].\n");
		return;
	}
	$out->print("E(eV),wl(nm),"
				."er(xx),ei(xx),n(xx),k(xx),a(xx),"
				."er(yy),ei(yy),n(yy),k(yy),a(yy),"
				."er(zz),ei(zz),n(zz),k(zz),a(zz),"
				."er(xy),ei(xy),er(yz),ei(yz),er(zx),ei(zx)\n");
	for(my $i = 0 ; ; $i++) {
		last if(!defined $pData[0]->[$i]);
		next if($pData[0]->[$i] <= 0.0);

		my $wl          = Optics::eVTonm($pData[0]->[$i]);
		my ($nxx, $kxx) = Optics::EpsToNK($pData[7]->[$i], $pData[1]->[$i]);
		my $axx         = Optics::KToAlpha($wl, $kxx);
		my ($nyy, $kyy) = Optics::EpsToNK($pData[8]->[$i], $pData[2]->[$i]);
		my $ayy         = Optics::KToAlpha($wl, $kyy);
		my ($nzz, $kzz) = Optics::EpsToNK($pData[9]->[$i], $pData[3]->[$i]);
		my $azz         = Optics::KToAlpha($wl, $kzz);
		$out->print("$pData[0]->[$i],$wl,"
					."$pData[7]->[$i],$pData[1]->[$i],$nxx,$kxx,$axx,"
					."$pData[8]->[$i],$pData[2]->[$i],$nyy,$kyy,$ayy,"
					."$pData[9]->[$i],$pData[3]->[$i],$nzz,$kzz,$azz,"
					."$pData[10]->[$i],$pData[4]->[$i],"
					."$pData[11]->[$i],$pData[5]->[$i],"
					."$pData[12]->[$i],$pData[6]->[$i]\n");
	}
	$out->Close();
	
	$App->print("OPTICS data saved to [$OpticsOut].\n");
}