#!/usr/bin/perl

#Start: use LauncherLibraries
use lib "/home/tkamiya/bin/Perl/lib";
use lib 'd:/Programs/Perl/lib';
#End: use LauncherLibraries

use strict;
use MultiColumnData;
use Crystal::VASP;

my $UseLog = 1;
my $Emin = -22;
my $Emax =  12;
my $Estep = 0.2;
my $OutFile = "DOS-up.out.csv";
my $EnergyFile = "Energy.out.csv";
#my $DOSFileTemplate = "DOS-a{}.csv";
#my $OUTCARTemplate  = "{}/SCF/OUTCAR";

my @files = glob("*");
my @FileHeaders = ();
for(my $i = 0 ; $i < @files ; $i++) {
	next if(!$files[$i] =~ /[\d+-\.eE]+$/);
	next if(!-d $files[$i]);
	next if(!-d "$files[$i]/DOS");
	push(@FileHeaders, $files[$i]);
}
@FileHeaders = sort { $a + 0.0 <=> $b + 0.0; } @FileHeaders;
print join(', ', @FileHeaders);
my @Values      = @FileHeaders;
my $nFiles  = @FileHeaders;
my $nE    = int( ($Emax - $Emin) / $Estep + 1.00001);
print "E: $Emin - $Emax, $Estep ($nE)\n";

my $VASP = new VASP();
my $out = new JFile($EnergyFile, "w") or die $!;
my $nIon;
for(my $i = 0 ; $i < @FileHeaders ; $i++) {
	my ($a) = ($FileHeaders[$i] =~ /([\d+\-\.eE]+)$/);
	my $f = "$FileHeaders[$i]/SCF/OUTCAR";
	$VASP->ReadOUTCARValues($f);
print "a=$a   E=$VASP->{FreeEnergy} eV\n";
	if($i == 0) {
		$nIon = $VASP->{"nIon"};
		$out->print("a,Total energy");
		for(my $iI = 1 ; $iI <= $nIon ; $iI++) {
			$out->print(",Ion$iI");
		}
		$out->print("\n");
	}

	$out->print("$a,$VASP->{FreeEnergy}");
	for(my $iI = 1 ; $iI <= $nIon ; $iI++) {
		$out->print(",$VASP->{Ion1_tot}");
	}
	$out->print("\n");
}
$out->Close();

my @in;
for(my $i = 0 ; $i < $nFiles ; $i++) {
	my $h = $FileHeaders[$i];
	my $f = "$h/DOS/DOS-up.csv"; #$InFile;
	$in[$i] = new MultiColumnData();
	my ($nData, $pLabelArray, @DataArray) = $in[$i]->Read($f, 1) or die $!;
	my $nDataArray = @$pLabelArray;
print "nData [$i]: $nData, $nDataArray, $pLabelArray\n";
}

my $out = new JFile($OutFile, "w");
die $! if(!$out);
my @K;
for(my $id = 0 ; $id < $nFiles ; $id++) {
	my $pX = $in[$id]->GetXData(0);
	my $pY = $in[$id]->GetYData(1);
	$K[$id] = 0;
	for(my $i = 0 ; $i < @$pY ; $i++) {
		$K[$id] = $pY->[$i] if($K[$id] < $pY->[$i]);
	}
	$K[$id] = 1.0 if($K[$id] < 1.0e-5);
}

for(my $i = 0 ; $i < $nE ; $i++) {
	my $E = $Emin + $Estep * $i;
	print ",$E";
	$out->print(",$E");
}
print("\n");
$out->print("\n");

my @Y;
for(my $id = 0 ; $id < $nFiles ; $id++) {
	my $pX = $in[$id]->GetXData(0);
	my $pY = $in[$id]->GetYData(1);

	$out->print("$Values[$id]");
	for(my $i = 0 ; $i < $nE ; $i++) {
		my $E = $Emin + $Estep * $i;
		$Y[$id] = $in[$id]->YVal($E) / $K[$id];
		$Y[$id] = 0.0 if($Y[$id] < 0.0);

		if($UseLog) {
			if($Y[$id] < 1.0e-5) {
				$Y[$id] = 0.0;
			}
			else {
				$Y[$id] = log($Y[$id]) / log(10.0) + 5.0;
			}
		}

		print ",$Y[$id]";
		$out->print(",$Y[$id]");
	}
	print("\n");
	$out->print("\n");
}
$out->Close();