#!/usr/bin/perl

use strict;

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';
use lib "$ENV{TkPerlDir}/lib";
use Deps;
use Utils;
use JFile;
use CSV;
use Sci qw($e $me $pi $h $hbar);
use Sci::Algorism;
use Crystal::VASP;

BEGIN {
#	$| = 1;
	open (STDERR, ">&STDOUT");
}

#===============================================
# ���ϐ�
#===============================================
my $ScriptCharCode = 'utf8';
my $FileSystemCharCode = 'sjis';

my $CARDir = ".";
$CARDir = $ARGV[0] if(defined $ARGV[0]);

my $nSmooth = 1;
my ($OutputEmin, $OutputEmax) = (-10.0, 10.0);

my $vasp = new VASP;
my $DOSDir = Utils::MakePath($CARDir, ["DOS"], '/', 0);
if(!-d $DOSDir) {
	$DOSDir = $CARDir;
}
my $POSCAR     = Utils::MakePath($DOSDir, ["POSCAR"], '/', 0);     #$vasp->GetPOSCARFileName($CARDir);
my $InputFile  = Utils::MakePath($DOSDir, ["DOS-up.csv"], '/', 0); #$vasp->GetPOSCARFileName($CARDir);
my $OutputFile = "DOS-me.csv";

print "nSmooth=$nSmooth\n";
print "POSCAR    : [$POSCAR]\n";
print "InputFile : [$InputFile]\n";
print "OutputFile: [$OutputFile]\n";

my $csv = new CSV;
if($csv->Read($InputFile)) {
	print("  Read [$InputFile]\n");
}
else {
	print("  Error!!: Can not read [$InputFile]\n");
	exit;
}

my $out = new JFile;
$out->Open($OutputFile, "w") or die "Error: Can not write to [$OutputFile].\n";

my $Crystal = $vasp->ReadStructureFromCARFiles($POSCAR, 1) or die "Error: Can not read [$POSCAR].\n";
my ($a, $b, $c, $alpha, $beta, $gamma) = $Crystal->LatticeParameters();
print "Latt: $a, $b, $c angstrom    $alpha, $beta, $gamma\n";
my $Volume = $Crystal->Volume();
print "Volume: $Volume A^3\n";

my $pLabelArray = $csv->LabelArray();
my $pDataArray  = $csv->DataArray();
my $nDataRow    = @$pDataArray;
my $pE          = $csv->GetXData("Energy.*");
my $pDOS        = $csv->GetXData("DOS.*");
print "pE  =$pE\n";
print "pDOS=$pDOS\n";
my $nData = @$pE;
print("  nDataRow: $nDataRow\n");
print("  nData   : $nData\n");

my (@Ne, @DOS2);
for(my $i= 0 ; $i < $nData ; $i++) {
	if($i == 0) {
		$Ne[0] = 0;
	}
	else {
		$Ne[$i] = $Ne[$i-1] + ($pDOS->[$i] + $pDOS->[$i-1]) / 2.0 * ($pE->[$i] - $pE->[$i-1]);
	}
	$DOS2[$i] = $pDOS->[$i] * $pDOS->[$i];
}
my $pdD2dE;
if($nSmooth > 1) {
	my $pDOS2Smoothed = Algorism::AdaptiveSmoothing(\@DOS2, $nSmooth);
	$pdD2dE = Algorism::DifferentiateByLinearAverage($pE, $pDOS2Smoothed); #\@DOS2);
}
else {
	$pdD2dE = [];
	$pdD2dE->[0] = ($DOS2[1]-$DOS2[0]) / ($pE->[1]-$pE->[0]);
	$pdD2dE->[$nData-1] = ($DOS2[$nData-1]-$DOS2[$nData-2]) / ($pE->[$nData-1]-$pE->[$nData-2]);
	for(my $i = 1 ; $i < $nData-1 ; $i++) {
		$pdD2dE->[$i] = ($DOS2[$i+1]-$DOS2[$i-1]) / ($pE->[$i+1]-$pE->[$i-1]);
	}
}

my $K = sqrt(2.0) / $pi / $pi / $hbar / $hbar / $hbar;
$out->print("E(eV),DOS(/uc),Ne,D0(cm-3eV-1.5),mde\n");
for(my $i = 0 ; $i < $nData ; $i++) {
	my $E = $pE->[$i];
	next if($E < $OutputEmin);
	next if($E > $OutputEmax);

	my $C = $pdD2dE->[$i] / ($Volume*$Volume*1.0e-48);
	my $sign = ($C > 0.0)? 1.0 : -1.0;
	my $D0 = $sign * sqrt(abs($C));
	my $mde = abs($D0) / $e**(3.0/2.0) * 1.0e6 / $K;
	$mde = $sign * $mde**(2.0/3.0) / $me;
	$out->print("$pE->[$i],$pDOS->[$i],$Ne[$i],$D0,$mde\n");
}

$out->Close();