#!/usr/bin/perl

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

use strict;
use Sci qw($pi $todeg $torad asin);
use Utils;
use IniFile;
use Template;
use Crystal::Crystal;
use Crystal::RietanFP;
use Sci::Optimize;

#=============================================
# 大域変数
#=============================================
my $TemplateIns = 'template.ins';
my $RietanDir   = 'D:\Programs\XRD\RIETAN-FP\RIETAN_VENUS';

my $Rietan = new RietanFP;

#=============================================
# 計算条件の初期値 (入力ファイルで書き換え可能)
#=============================================
my $ini = new IniFile;
$ini->{ScanBiso} = 0;
($ini->{Biso0}, $ini->{Biso1}, $ini->{BisoStep}) = (0.0, 1.0, 0.01);
$ini->{ScanzAs}  = 1;
($ini->{zAs0},  $ini->{zAs1},  $ini->{zAsStep})  = (0.3, 0.4, 0.01);

#$ini->{CEngine} = 'RIETAN|BUILTIN';
#$ini->{CEngine} = 'RIETAN';
$ini->{CEngine} = 'BUILTIN';

#=============================================
# ファイルパス設定
#=============================================
my $RietanHeader = $ARGV[0];
if($RietanHeader eq '') {
	print "\nusage: perl EpiXRD.pl File_Header\n\n";
	exit 0;
}

my ($drive, $directory, $filename, $ext1, $lastdir, $filebody) = Deps::SplitFilePath($RietanHeader);
$RietanHeader = $filebody;
my $IniFile   = "${RietanHeader}.ini";

if(!-f $IniFile) {
	print "\nError: File [$IniFile] does not exist.\n\n";
	exit 0;
}

my $InsFile = "${RietanHeader}.ins";
my $LstFile = "${RietanHeader}.lst";
my $cmd     = "RIETAN-FP.bat ${RietanHeader}";

$ENV{RIETAN}  = "${RietanDir}\\programs";

#=============================================
# 設定ファイル読み込み
#=============================================
print "Read [$IniFile]\n";
$ini->ReadAll($IniFile, '=', 0) or die "$!: Can not read [$IniFile].\n";

if($ini->{HKLM1} eq '') {
	print "\nHermann-Mauguin symbol is not specified with the key \"HKLM1\".\n";
	my ($iSPG, $iSet) = ($ini->{SPGName} =~ /(\d+)-(\d+)/);
	print "  Try to find the space group #$iSPG and the set #$iSet.\n";
	$ini->{HKLM1} = RietanFP->ReadHermannMauguin($iSPG, $iSet);
#$ini->{HKLM1} = 'I 4/m m m';
	print "  Use \"$ini->{HKLM1}\" for HKLM1.\n\n";
}

#foreach my $key (sort keys %$ini) {
#	print "$key: $ini->{$key}\n";
#}

# zAs, Bisoの初期値の読み込み
my @a = Utils::Split("\\s+", $ini->{Site4});
$ini->{zAs} = $a[4];
$ini->{Biso} = $a[5];

my $Crystal = &BuildCrystal($ini);

#=============================================
# 入力のリピート
#=============================================
print "Input data\n";
print "  CEngine: $ini->{CEngine}\n";
print "  zAs =$ini->{zAs}\n";
print "  Biso=$ini->{Biso}\n";
print "  Optimize=$ini->{Optimize} by $ini->{Method} EPS=$ini->{EPS} nMaxIter=$ini->{nMaxIter}iPrintLevel=$ini->{iPrintLevel}\n";
print "  ScanBiso=$ini->{ScanBiso} (Biso from $ini->{Biso0} to $ini->{Biso1} at step $ini->{BisoStep})\n";
print "  ScanzAs =$ini->{ScanzAs} (zAs from $ini->{zAs0} to $ini->{zAs1} at step $ini->{zAsStep})\n";
print "  Observed Intensities\n";
my $I0obs = $ini->{$ini->{NormalizeKey}};
foreach my $key (keys %$ini) {
	next if($key !~ /^I[\-\d]+$/);

	printf "    $key: %8.2f: ", $ini->{$key};
	$ini->{$key} /= $I0obs;
	printf "    %8.6f\n", $ini->{$key};
}
print "\n";

#=============================================
# 計算実行
#=============================================
my $PMonochro = cos($ini->{Monochro2QB} * $torad);
if($ini->{MonochroType} eq 'SingleCrystal') {
}
else {
	$PMonochro = $PMonochro * $PMonochro;
}
my $mud = $ini->{mu} * $ini->{Thickness} * 1.0e-7; # mu in cm-1, thickness in nm

# Optimize
if($ini->{Optimize}) {
	print "=== Optimize by $ini->{Method} method\n";
	my $optimize = new Optimize;;
	if($ini->{Method} =~ /Simplex/i) {
		$optimize->AddParameters(
			"zAs",  \$ini->{zAs},   1, 0.05, 0.0, 1.0,
			sub {
				$ini->{zAs} = $_[0];
			},
			"Biso",  \$ini->{Biso}, 1, 0.1, 0.0, 2.0,
			sub {
				$ini->{Biso} = $_[0];
			},
			);
	}
	else {
		$optimize->{diffv}              = 0.001;
		$optimize->{DumpingFactorStart} = 10.0;
		$optimize->{DumpingFactorStep}  = -1.0;
		$optimize->{nDumpingFactor}     = 10;
		$optimize->AddParameters(
			"zAs",  \$ini->{zAs},   1, 0.01, 0.0, 1.0,
			sub {
				$ini->{zAs} = $_[0];
			},
			"Biso",  \$ini->{Biso}, 1, 0.01, 0.0, 2.0,
			sub {
				$ini->{Biso} = $_[0];
			},
			);
	}

#	$optimize->SetnS2Calculation(0);
	my ($OptVars, $MinVal) = $optimize->Optimize(
			$ini->{Method}, undef, undef, undef,
			$ini->{EPS}, $ini->{nMaxIter}, $ini->{iPrintLevel},
			sub { return &CalS2Opt(@_); },
			undef,
			sub { Optimize::BuildDifferentialMatrixes(@_); },
		);
	print "\nOptimized at:\n";

	$optimize->RecoverParameters($OptVars);
	$optimize->PrintParameters(1, $OptVars, $MinVal, 1);

	$ini->{zAs}  = $OptVars->[0];
	$ini->{Biso} = $OptVars->[1];

#=============================================
# 最小のS2をもつ結果の再計算と表示
#=============================================
	print "\n=== Final result candidate ===\n";
	print "Minmum S2 = $MinVal at zAs = $ini->{zAs}, Biso = $ini->{Biso}\n";
	my $S2 = &CalS2($ini->{zAs}, $ini->{Biso}, IsPrint => 1);
	printf "zAs=%8.6f, Biso=%8.6f: S2=%10.6g\n\n", $ini->{zAs}, $ini->{Biso}, $S2;
}

# Bisoのスキャン
if($ini->{ScanBiso}) {
	my $MinBiso;
	my $MinS2 = 1.0e10;
	print "=== Scanning for Biso from $ini->{Biso0} to $ini->{Biso1}\n";
	for($ini->{Biso} = $ini->{Biso0} ; $ini->{Biso} <= $ini->{Biso1} ; $ini->{Biso} += $ini->{BisoStep}) {
		my $S2 = &CalS2($ini->{zAs}, $ini->{Biso}, IsPrint => 0);
		printf "  Biso=%5.3f, zAs=%5.3f: S2=%10.6g\n", $ini->{Biso}, $ini->{zAs}, $S2;

		if($S2 == 0.0) {
			print "\nError: RIETAN-FP might be terminated incorrectly.\n";
			print "   The line including \"*** Summary of possible reflections\" could not be found\n";
			print "     in the output file \"$LstFile\"\n\n";
			exit;
		}

		if($MinS2 > $S2) {
			$MinS2 = $S2;
			$MinBiso = $ini->{Biso};
		}
	}

#=============================================
# 最小のS2をもつ結果の再計算と表示
#=============================================
	print "\n=== Final result candidate ===\n";
	print "Minmum S2 = $MinS2 at Biso = $MinBiso, zAs = $ini->{zAs}\n";
	$ini->{Biso} = $MinBiso;
	my $S2 = &CalS2($ini->{zAs}, $ini->{Biso}, IsPrint => 1);
	print "Biso=$ini->{Biso}: S2=$S2\n\n";
}

# zAsのスキャン
if($ini->{ScanzAs}) {
	my $MinzAs;
	my $MinS2 = 1.0e10;
	print "=== Scanning for zAs from $ini->{zAs0} to $ini->{zAs1}\n";
	for($ini->{zAs} = $ini->{zAs0} ; $ini->{zAs} <= $ini->{zAs1} ; $ini->{zAs} += $ini->{zAsStep}) {
		my $S2 = &CalS2($ini->{zAs}, $ini->{Biso}, IsPrint => 0);
		printf "zAs=%5.3f, Biso=%5.3f: S2=%10.6g\n", $ini->{zAs}, $ini->{Biso}, $S2;

		if($S2 == 0.0) {
			print "\nError: RIETAN-FP might be terminated incorrectly.\n";
			print "   The line including \"*** Summary of possible reflections\" could not be found\n";
			print "     in the output file \"$LstFile\"\n\n";
			exit;
		}

		if($MinS2 > $S2) {
			$MinS2 = $S2;
			$MinzAs = $ini->{zAs};
		}
	}

#=============================================
# 最小のS2をもつ結果の再計算と表示
#=============================================
	print "\n=== Final result candidate ===\n";
	print "Minmum S2 = $MinS2 at zAs = $MinzAs, Biso = $ini->{Biso}\n";
	$ini->{zAs} = $MinzAs;
	my $S2 = &CalS2($ini->{zAs}, $ini->{Biso}, IsPrint => 1);
	printf "zAs=%5.3f, Biso=%5.3f: S2=%10.6g\n\n", $ini->{zAs}, $ini->{Biso}, $S2;
}

exit;

#=============================================
# Subroutines
#=============================================
sub BuildCrystal
{
	my ($ini) = @_;

	my $Crystal = new Crystal;

print "Build Crystal\n";
	$Crystal->SetSampleName($ini->{Sample});
print "  Sample: $ini->{Sample}\n";
	$Crystal->SetCrystalName($ini->{Sample});
	my ($iSPG, $iSet) = ($ini->{SPGName} =~ /(\d+)-(\d+)/);
print "  SPG: $iSPG, $iSet\n";
	$Crystal->SetSpaceGroup(undef, $iSPG, $iSet);
	my $SPG = $Rietan->ReadSpaceGroup($iSPG, $iSet);
	$Crystal->SetCSpaceGroup($SPG);

	$Crystal->SetLatticeParameters($ini->{a}, $ini->{b}, $ini->{c}, $ini->{alpha}, $ini->{beta}, $ini->{gamma});
print "  Latt: $ini->{a}, $ini->{b}, $ini->{c}, $ini->{alpha}, $ini->{beta}, $ini->{gamma}\n";
	my @a = Utils::Split('\\s+', $ini->{Elements});
	for(my $i = 0 ; $i < @a ; $i++) {
		$a[$i] =~ s/\'//g;
print "  Atom$i: $a[$i]\n";
		$Crystal->AddAtomType($a[$i]);
	}
	for(my $i = 1 ; ; $i++) {
		my $key = "Site$i";
		last if(!defined $ini->{$key});
		
		my ($label, $atomname, $occ, $x, $y, $z) = ($ini->{$key} =~ /\s*(\S+)?\/(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)/);
print "  Site$i: $label $atomname  ($x, $y, $z) occ=$occ\n";
		$Crystal->AddAtomSite($label, $atomname, $x, $y, $z, $occ);
	}

	$Crystal->ExpandCoordinates();
my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
	my $site = $ExpandedAtomSiteList[$i];
	my ($x, $y, $z) = $site->Position();
	my $atomname    = $site->AtomName();
	print "  Expanded site$i: $atomname ($x, $y, $z)\n";
}
print "\n";

	return $Crystal;
}

sub CalS2Opt
{
	my ($pVars, $iPrintLevel) = @_;
	my ($zAs, $Biso) = @$pVars;

	my $S2 = &CalS2($zAs, $Biso, IsPrint => 0);
	if($iPrintLevel) {
		printf "zAs=%8.6f, Biso=%8.6f: S2=%10.6g\n", $zAs, $Biso, $S2;
	}
	return $S2;
}

sub CalS2
{
	my ($zAs, $Biso, %args) = @_;

	if($args{IsPrint}) {
		print "\n";
		print "Calculated Intensities for zAs=$zAs, Biso=$Biso\n";
	}

	my $mu = $ini->{mu};
	my $d  = $ini->{Thickness};

	my $pPeaks;
	my %Peaks;
	if($ini->{CEngine} =~ /RIETAN/i) {
			$pPeaks = &GetFbyRIETAN($zAs, $Biso);
#print "nPeaks=", scalar @$pPeaks, "\n";

		for(my $i = 0 ; $i < @$pPeaks ; $i++) {
			my $h    = $pPeaks->[$i]{h};
			my $k    = $pPeaks->[$i]{k};
			my $l    = $pPeaks->[$i]{l};
			my $Q2   = $pPeaks->[$i]{Q2};
			my $dhkl = $pPeaks->[$i]{dhkl};
			my $F    = $pPeaks->[$i]{F};

			$Peaks{"I$h$k$l"} = {
				h    => $h,
				k    => $k,
				l    => $l,
				Q2   => $Q2,
				dhkl => $dhkl,
				F    => $F,
				};
		}
	}

	my %Ical;
	my $c = 0;
	foreach my $key (sort keys %$ini) {
		next if($key !~ /^I\d+/);
		
		my ($h, $k, $l) = ($key =~ /I(\d)(\d)(\d)/);

		my ($Q2, $dhkl, $F);
		if($ini->{CEngine} =~ /RIETAN/i) {
			my $p = $Peaks{"I$h$k$l"};
			$Q2   = $p->{Q2};
			$dhkl = $p->{dhkl};
			$F    = $p->{F};
			my $s = 1.0 / 2.0 / ($dhkl * 0.1);
			my $T = CrystalObject::TemperatureFactor($Biso, $s);
			$F    /= $T;
printf "  RIETAN  : F$h$k$l @ %8.6f F = %8.6f (F*T=%8.6f) \n", $Q2, $F, $F*$T if($ini->{Debug});
		}
		if($ini->{CEngine} =~ /BUILTIN/i) {
			my $RefAtomList = $Crystal->{RefAsymmetricAtomSiteList};
			my $site =  $RefAtomList->[3];
			my ($x, $y, $z) = $site->Position();
			$site->SetPosition($x, $y, $zAs);
			$Crystal->ExpandCoordinates();

			$dhkl     = $Crystal->CalculateHKLInterplanarSpacing($h, $k, $l); # in angstrom
			my $sinQB = $ini->{lambda} * 10.0 / 2.0 / $dhkl;
			next if($sinQB >= 1.0);

			$Q2 = 2.0 * asin($sinQB) * $todeg;
			my ($Fr, $Fi) = $Crystal->Fhkl($h, $k, $l);
			$F  = sqrt($Fr*$Fr + $Fi*$Fi);
printf "  Built-in: F$h$k$l @ %8.6f F = %8.6f\n", $Q2, $F if($ini->{Debug});
		}
#print "c2=$c  $h $k $l $Q2 $dhkl $F\n";
		$c++;

# Ihklの計算
		my $sin2Q = sin($Q2 * $torad);
		my $sinQ  = sin($Q2 / 2.0 * $torad);
		my $cosQ  = cos($Q2 / 2.0 * $torad);
		my $s     = $sinQ / $ini->{lambda}; # nm-1
		my $A     = 1.0 - exp(-2.0 * $mud / $sinQ);
		my $P     = 1.0 + $cosQ * $cosQ * $PMonochro;
		my $L     = 1.0 / $sin2Q;
		my $T     = CrystalObject::TemperatureFactor($Biso, $s);
		my $I     = $F * $F * $P * $L * $T * $T * $A;
		$Ical{$key} = $I;

		if($args{IsPrint}) {
			printf "F(%d%d%d)=|%7.3f| d=%7.3f  I=%8.3f\n", $h, $k, $l, $F, $dhkl, $I;
			printf "  P=%7.3f  L=%7.3f  T=%7.3f  A=%7.3f\n", $P, $L, $T, $A;
#			print "$h $k $l   2Q=$Q2  F=$F\n" if($Code == 1 and $h == 0 and $k == 0 and $l <= 8);
		}
	}

	if($args{IsPrint}) {
		print "\n";
		print "Compare (hkl: Iobs Ical)\n";
	}

# 各IhklについてS2を計算
	my $S2 = 0.0;
	my $I0cal = $Ical{$ini->{NormalizeKey}};
	foreach my $key (sort keys %Ical) {
		$Ical{$key} /= $I0cal;
		my $d = ($Ical{$key} - $ini->{$key}) / $ini->{$key};
		$S2 += $d*$d;

		if($args{IsPrint}) {
			printf "$key: %10.6f %10.6f\n", $ini->{$key}, $Ical{$key};
		}
	}

	$S2 += 1e5 * abs($Biso) if($Biso <= 0.0);

	return $S2;
}


sub GetFbyRIETAN
{
	my ($zAs, $Biso) = @_;

#=============================================
# テンプレートINSファイルからINSファイルの作成
#=============================================
	for(my $i = 1 ; ; $i++) {
		my $key = "Site$i";
		last if(!defined $ini->{$key});

		my @a = Utils::Split("\\s+", $ini->{$key});
#Site1=Ba1/Ba2+  1.0    0  0    0      0.7500  01101
		if($i == 4) {
			$ini->{$key} = "$a[0]    $a[1]    $a[2]  $a[3]  $zAs   $Biso  $a[6]";
		}
		else {
			$ini->{$key} = "$a[0]    $a[1]    $a[2]  $a[3]  $a[4]   $Biso  $a[6]";
		}
	}

	$ini->{StartAngle} =   5.0;
	$ini->{EndAngle}   = 150.0;
	$ini->{StepAngle}  =   0.01;
	my $text = Template->new()->GetTextByHash($TemplateIns, $ini, '{', '}', 0, undef, 1);
	open(OUT, ">$InsFile") or die "$!\n";
	print OUT $text;
	close(OUT);

#=============================================
# RIETAN-FPの実行
#=============================================
	if(system($cmd) >= 32) {
		die "Can not execute [$cmd]\n";
	}

#=============================================
# RIETAN-FPの出力ファイルの読み込み
#=============================================
#print "Read [$LstFile] (zAs=$zAs, Biso=$Biso)\n";
	my $in = new JFile($LstFile, 'r');
	if(!$in) {
		die "Can not read [$LstFile].\n";
	}

# |F|行の検索
	my $line = $in->SkipTo("Summary of possible reflections");
	$line = $in->SkipTo("No.\\s+Phase\\s+h\\s+k\\s+l");

	my @peaks;
	my $c = 0;
	while(1) {
		my $line = $in->ReadLine();
		my @a = Utils::Split("\\s+", $line);
		last if(@a <= 10);

# hkl, |F|の読み込み
#            No.   Phase   h    k    l    Code   2-theta        d       Ical  |F(nucl)| |F(magn)|    POF     FWHM    m     Dd/d
		my ($No, $Phase, $h, $k, $l, $Code, $Q2, $dhkl, $Ical, $F, $Fmagn, $POF, $FWHM, $m) = @a;
		last if(!defined $F);
		next if($Code > 1);
#print "c1=$c  $h $k $l $Q2 $dhkl $F\n";

		$peaks[$c] = {
			Phase => $Phase,
			h     => $h,
			k     => $k,
			l     => $l,
			Q2    => $Q2,
			dhkl  => $dhkl,
			Ical  => $Ical,
			F     => $F,
			Fmagn => $Fmagn,
			POF   => $POF,
			FWHM  => $FWHM,
			m     => $m
			};
		$c++;
	}
	$in->Close();

	return \@peaks;
}