#!/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 Utils;
use MyApplication;
use Crystal::VASP;
use Crystal::WIEN2k;
use Sci qw($RyToeV);

#===============================================
# スクリプト大域変数
#===============================================
my $dE = 2.7; # eV
my $DataDir = $ARGV[0];
$DataDir = "." if($DataDir eq '');
$DataDir = Utils::ReduceDirectory($DataDir);
print "DataDir: $DataDir\n";

my $INCARFile = "INCAR";

#===============================================
# メインルーチン
#===============================================
my $WIEN2k = new WIEN2k;

my $Title                 = &GetTitle($DataDir);
my $StructPath            = Deps::MakePath($DataDir, "$Title.struct");
my $InEnergyPath          = "SrGeO3-Wien2k.energy";
my $OutEnergyPath         = "SrGeO3-Wien2k.energy.scissors";

my $Crystal = $WIEN2k->ReadStructFile($StructPath, 0);#, tol => 1.0e-2);
my $LatticeSystem = $Crystal->LatticeSystem();
print "Lattice system in WIEN2k: $LatticeSystem\n";

my $EF = $WIEN2k->ReadFermiEnergy($StructPath, 1);
print "EF=$EF eV\n";

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

my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
my $nSite = @AtomSiteList;
for(my $i = 0 ; $i < $nSite ; $i++) {
	my $line = $in->ReadLine();
	$out->print($line);
	$line = $in->ReadLine();
	$out->print($line);
}

my $EFidx = 0;
my $pos = $in->tell();
my $line = $in->ReadLine();
last if(!defined $line);
#print "l: $line";
my ($nl, $w) = ($line =~ /(\d+)\s+([\d\.]+)\s*$/);
print "nl=$nl  w=$w\n";
#exit;
for(my $i = 0 ; $i < $nl ; $i++) {
	my $line = $in->ReadLine();
	last if(!defined $line);
	
	my ($idx, $e) = Utils::Split("\\s+", $line);
	my $eeV = $e*$RyToeV;
	if($eeV > $EF) {
print "E [$eeV eV / $e Ry] exceeds EF [$EF eV] for index=$idx\n";
		$EFidx = $idx;
		last;
	}
}
#exit;

$in->seek($pos, 0);
while(1) {
	my $line = $in->ReadLine();
	last if(!defined $line);
	$out->print($line);
#print "l: $line";
	my ($nl, $w) = ($line =~ /(\d+)\s+([\d\.]+)\s*$/);
print "nl=$nl  w=$w\n";
	for(my $i = 0 ; $i < $nl ; $i++) {
		my $line = $in->ReadLine();
		last if(!defined $line);

		my ($idx, $e) = Utils::Split("\\s+", $line);
		if($idx < $EFidx) {
			$out->print($line);
		}
		else {
			$e += $dE / $RyToeV;
			$out->printf("%12d %23.15e\n", $idx, $e);
		}
	}
}

$out->Close();
$in->Close();

exit;


sub GetTitle
{
	my ($path) = @_;

	my ($drive, $directory, $filename, $ext, $lastdir, $filebody) = Deps::SplitFilePath($path);
	my $Title   = $filebody;
	$Title =~ s/\s/_/g;
	return $Title;
}