#!/usr/bin/perl

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';

use strict;

use Utils;
use Sci::Algorism;
use JFile;
use CSV;

use Sci qw($pi $kB $R $NA $c $e $e0 $u0 $me $mp $mn $h $hbar $torad $todeg);

my $SourceCSV  = "TotalEnergy-BaSiO3-cubic-PBE96-hard.csv";
   $SourceCSV  = $ARGV[0] if(defined $ARGV[0]);
my $SourceCSV2 = "TotalEnergy-BaSiO3-P63mmc-PBE96-hard.csv";
   $SourceCSV2 = $ARGV[1] if(defined $ARGV[1]);

my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($SourceCSV2);
my $OutCSV    = Deps::MakePath("$drive$directory", "${filebody}-Out.csv", 0);
my $FitCSV    = Deps::MakePath("$drive$directory", "${filebody}-Fit.csv", 0);

my ($pInf,  $Pmin,  $Pmax,  $pP,  $pV,  $pE,  $pH)  = &ReadCSV($SourceCSV);
my ($pInf2, $Pmin2, $Pmax2, $pP2, $pV2, $pE2, $pH2) = &ReadCSV($SourceCSV2);


my $out = JFile->new($OutCSV, "w") or "$!: Can not write to [$OutCSV].\n";
$out->print("P(GPa),V(A^3),Emol(eV),Hmol(eV),V2(A^3),Emol2(eV),Hmol2(eV),dHmol(eV)\n");
#$out->print("P(GPa),a(A),b(A),c(A),V(A^3),Emol(eV),Hmol(eV),dHmol(eV)\n");
for(my $i = 0 ; $i < @$pP ; $i++) {
	my $P = $pP->[$i];
	next if($P < $Pmin2 or $Pmax2 < $P);
	my $V = $pV->[$i];
	my $E = $pE->[$i];
	my $H = $pH->[$i];

	my $V2 = Algorism::InterpolateBySimpson($pP2, $pV2, $P);
	my $E2 = Algorism::InterpolateBySimpson($pP2, $pE2, $P);
	my $H2 = Algorism::InterpolateBySimpson($pP2, $pH2, $P);
	my $dH = $H2 - $H;

	$out->print("$P,$V,$E,$H,$V2,$E2,$H2,$dH\n");
#	$out->print("$P,$pInf->[$i]{a},$pInf->[$i]{b},$pInf->[$i]{c},$V,$E,$H\n");
}
$out->Close();


exit;

#========================================================================
# Subroutines
#========================================================================
sub ReadCSV
{
	my ($SourceCSV) = @_;

	my $csv = new CSV;
	$csv->Read($SourceCSV) or die "$!: Can not read [$SourceCSV].\n";
	$csv->Close();
	my $pZ = $csv->GetXData("Z");
	my $pa = $csv->GetXData("a");
	my $pb = $csv->GetXData("b");
	my $pc = $csv->GetXData("c");
	my $pV = $csv->GetXData("V.*");
	my $pP = $csv->GetXData("PStress.*");
	my $pEtot = $csv->GetXData("Etot.*");
	my $pEmol = $csv->GetXData("Emol.*");
	my $nData = @$pV;
	my ($Pmin, $Pmax) = Utils::CalMinMax($pP, 0);

print "Source: [$SourceCSV]\n";
print "Out   : [$OutCSV]\n";
print "pa=$pa\n";
print "pb=$pb\n";
print "pc=$pc\n";
print "pV=$pV\n";
print "pP=$pP\n";
print "pE=$pEmol\n";
print "nData=$nData\n";
print "P = ($Pmin - $Pmax)\n";

	my @H;
	my @Inf;
	for(my $i = 0 ; $i < $nData ; $i++) {
		$H[$i] = ($pEtot->[$i] + $pP->[$i]*0.1 * 1.0e9 * $pV->[$i] * 1.0e-30 / $e) / $pZ->[$i];
print "$i: P=$pP->[$i] V=$pV->[$i]  E= $pEmol->[$i]  H= $H[$i]\n";
		$Inf[$i] = {
			Z    => $pZ->[$i],
			a    => $pa->[$i],
			b    => $pb->[$i],
			c    => $pc->[$i],
			V    => $pV->[$i],
			Vmol => $pV->[$i] / $pZ->[$i],
			P    => $pP->[$i] * 0.1, # kbar => GPa
			E    => $pEmol->[$i],
			Emol => $pEmol->[$i],
			Etot => $pEtot->[$i],
			Hmol => $H[$i],
			};
	}
	my @sorted = sort { return $a->{P} <=> $b->{P} if(defined $a->{P});
						return $b->{V} <=> $a->{V} } @Inf;

	for(my $i = 0 ; $i < @sorted ; $i++) {
#print "i=$i: P=$sorted[$i]{P}, V=$sorted[$i]{V}\n";
		$pP->[$i] = $sorted[$i]->{P};
		$pV->[$i] = $sorted[$i]->{V};
		$pE->[$i] = $sorted[$i]->{Emol};
		$H[$i]    = $sorted[$i]->{Hmol};
	}
#exit;

	return (\@sorted,$Pmin*0.1, $Pmax*0.1, $pP, $pV, $pE, \@H);
}