#!/usr/bin/perl

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

use strict;
#use warnings;

use Sci qw($c $e $kB $NA $me $mn $hbar $h $pi $pi2);
use Sci::Algorism;

#===============================================
# デバッグ関係変数

my $Crystal = 'HT';
#my $Crystal = 'LT';

my ($a, $b, $c);
my $nAtom;

#potential in frac
# U in eV   an => an * latt [m]^n [in eV/
#         a0   a1   a2      a3   a4      a5   a6
my @Ua;

if($Crystal eq 'HT') {
# HT cubic Pm-3m
	($a, $b, $c) = (3.477, 3.477, 3.477);  # A
	$nAtom     = 1; # HT
	@Ua = (0.0, 0.0, 0.0173, 0.0, 0.0141, 0.0, 0.0); # HT
}
elsif($Crystal eq 'LT') {
# LT Pnma
	($a, $b, $c) = (7.285,  4.482,   5.506);
#Rietveld
#	($a, $b, $c) = (7.1117, 4.50672, 5.40546);
	$nAtom     = 4; # LT
	@Ua = (0.0, 0.0, 0.0373, 0.0, 0.0646, 0.0, 0.01292); # HT
}

my $atomname  = 'Cu';
my $M         =  63.546;
my $Q         = 0.0;
my $direction = '<110>';
#my $Kx        = sqrt($a*$a + $b*$b) * 1.0e-10;
my $Kx        = 1.0e-10;

my ($Tmin, $Tmax, $Tstep) = (0.0, 1200.0, 100.0);

# Integration
my $EPS         = 1.0e-10;
my $xmaxscale   = 10.0;
#my $IntegMethod = 'Rieman';
my $IntegMethod = 'IMT';
my $nx          = 501;

#==================================
# Main
#==================================

print "\n";
print "Einstein model F calculation for phase [$Crystal].\n";
print "\n";
print "Lattice parameters: $a $b $c A\n";
print "Displacement conversion factor to [m] $Kx\n";
#print "d$direction = $Kx A\n";
print "Oscillating atom: $atomname  M = $M\n";
printf "Uij(r) = %5.3f + %5.3f * r + %5.3f * r^2 + %5.3f * r^3 + %5.3f * r^4"
	        ." + %5.3f * r^5 + %5.3f * r^6\n", @Ua;
print "\n";

#F = m x'' ~ -2a2 * x
#x = x0 * exp(wt)
#  -mw2 * x ~ -2a2 * x
#  w = sqrt(2a2/m)
my $a2 = $Ua[2] * $e / $Kx**2; #J / m^2
my $omega = sqrt(2.0 * $a2 / ($M*$mn));
my $f     = $omega / $pi2;
my $lamda = $c / ($omega / $pi2); # m
my $k     = 1.0 / ($lamda * 100.0); # 1/cm
my $E     = $hbar * $omega / $e;
printf "Omega = %e /s\n",   $omega;
printf "f     = %g THz\n",  $f * 1.0e-12;
printf "lamda = %e m\n",    $lamda;
printf "k     = %g 1/cm\n", $k;
printf "E     = %g meV\n",  $E * 1.0e3;
print "\n";

print "F harmonic vibration from the analytical equation: Fvib = Ew/2 + ln[1 - exp(-Ew/kT)]\n";
print " T(K)     Fharm (eV)     Fharm (kJ/mol)\n";
for(my $i = 0 ; ; $i++) {
	my $T = $Tmin + $i * $Tstep;
	last if($T > $Tmax + 1.0e-6);
	next if($T == 0.0);
	
	my $kBTe = $kB * $T / $e;
	my $f = 1.0 - exp(-$E / $kBTe);
	my $lnf = log($f);
#print "  f=$f  lnf=$lnf ln(10)=", log(10.0), "\n";
	my $F = 0.5 * $E + $kBTe * $lnf * $nAtom; # eV
	my $F2 = $F * $e * $NA; # J/mol
	printf "%8.3f %12.3g %12.3g\n", $T, $F, $F2 * 1.0e-3;
}
print "\n";

print "Integration accuracy check\n";
my $Ttest = 1000.0; # K
my $xmax = $xmaxscale * (-log($EPS) * ($kB * $Ttest / $e))**(1.0/2.0);
#my $xmax = sqrt(-log($EPS) / $Ua[2] * ($kB * $Ttest / $e));

print "Rieman sum\n";
for(my $nxtest  = 10 ; $nxtest <= 10001 ; $nxtest *= 2) {
	my $ZU = 2.0 * Algorism::IntegrateByRectangleWithFunction(
		sub { 
			my ($x) = @_;
			my $Uij = &Uij($x, \@Ua);
			my $y   = exp(-$Uij / ($kB * $Ttest / $e));
			return $y;
			},
		0.0, $xmax, $nxtest,
		);
	printf "  nx=%6d  ZU = %12.8g (xmax = %8.3f for T=$Ttest EPS=$EPS)\n", $nxtest, $ZU, $xmax;
}
print "\n";

print "IMT integral\n";
for(my $nxtest  = 10 ; $nxtest <= 10001 ; $nxtest *= 2) {
	my $ZU = Algorism::IntegrateByIMTWithFunction(
		sub { 
			my ($x) = @_;
			my $Uij = &Uij($x, \@Ua);
			my $y   = exp(-$Uij / ($kB * $Ttest / $e));
			return $y;
			},
		-$xmax, $xmax, $nxtest,
		);
	printf "  nx=%6d  ZU = %12.8g (xmax = %8.3f for T=$Ttest EPS=$EPS)\n", $nxtest, $ZU, $xmax;
}
print "\n";

print "F vibration from numerical integral: Fvib = ln Zvib, Zvib = int_x exp(U(x)/kT) / d$direction\n";
printf "  EPS=%8.3g  nx=%d\n", $EPS, $nx;
printf "  Integration by $IntegMethod\n";
#print " T(K)     Fvib (eV)     Fvib (kJ/mol)\n";
print " T(K)     Fkinetic (eV)   Fvib (kJ/mol)\n";
for(my $i = 0 ; ; $i++) {
	my $T = $Tmin + $i * $Tstep;
	last if($T > $Tmax + 1.0e-6);
	next if($T == 0.0);
	
	my $kBTe = $kB * $T / $e;

	my $xmax;
	my (@Z, @F, @F2);
	my $FK = -$kB * $T / $e * 0.5 * log($pi2 * $kB * $T / ($M * $mn)) * $nAtom; # eV
	printf "%8.2f ", $T;
	printf "  FK=%12.3g", $FK;
	for(my $i = 2 ; $i < @Ua ; $i += 2) {
		$xmax = $xmaxscale * (-log($EPS) * ($kB * $Ttest / $e))**(1.0/2.0);
#		$xmax = $xmaxscale * (-log($EPS) * ($kB * $Ttest / $e))**(1.0/$i);
#		$xmax = sqrt(-log($EPS) / $Ua[2] * ($kB * $T / $e));
#print "\nxmax[$i]=$xmax\n";

		if($IntegMethod eq 'IMT') {
			$Z[$i] = Algorism::IntegrateByIMTWithFunction(
			sub { 
				my ($x) = @_;
				my $Uij = $Ua[$i] * $x**$i;
				my $y   = exp(-$Uij / ($kB * $T / $e));
				return $y;
				},
			-$xmax, $xmax, $nx,
			);
		}
		elsif($IntegMethod eq 'Rieman') {
			$Z[$i] = 2.0 * Algorism::IntegrateByRectangleWithFunction(
			sub { 
				my ($x) = @_;
				my $Uij = $Ua[$i] * $x**$i;
				my $y   = exp(-$Uij / ($kB * $T / $e));
				return $y;
				},
			0.0, $xmax, $nx,
			);
		}

# convert to m, multiply ZK = sqrt(pi 2kBT / m)
		$Z[$i] *= $Kx;

		$F[$i]  = -$kB * $T / $e * log($Z[$i]) * $nAtom; # eV
		$F2[$i] = $F[$i] * $e * $NA; # J/mol
#		printf "  Z$i=%12.3g", $Z[$i];
		printf "  F$i=%12.3g", $F[$i];
#		printf "  F$i=%12.3g", $F2$i];
	}
	printf " (xmax=%6.3f)\n", $xmax;
}
print "\n";

exit;


#======================================
# Subroutine
#======================================
sub Uij
{
	my ($r, $pUa) = @_;
	
	my $rn = 1.0;
	my $Uij = 0.0;
	for(my $i = 0 ; $i < @$pUa ; $i++) {
		$Uij += $pUa->[$i] * $rn;
		$rn *= $r;
	}
	return $Uij;
}

sub Fij
{
	my ($r, $pUa) = @_;
	
	my $rn = 1.0;
	my $Fij = 0.0;
	for(my $i = 1 ; $i < @$pUa ; $i++) {
		$Fij += -$i * $pUa->[$i] * $rn;
		$rn *= $r;
	}
	return $Fij;
}