#!/usr/bin/perl

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

use strict;
use Sci qw($h $hbar $me $e $todeg $torad asin tan);
use Crystal::CIF;
use Crystal::Crystal;
use Crystal::Rietan;

my $ElectronEnergy  = 100.0e3; # eV
my $lambda = sqrt($h * $h / 2.0 / $me / $ElectronEnergy / $e) * 1.0e9; # in nm
print "E=$ElectronEnergy eV: lambda = $lambda nm\n";

my $CIFFile = "d:/MyWebs/Research/CrystalStructure/ZnO.cif"; #"Si5.cif";
my $OutFile = "TED.csv";
my $Biso  = 0.5; #1.0;

my $Q2Max = 3.0;
my $Qcriteria = 0.1;

my $hmax = 10;
my $kmax = 10;
my $lmax = 10;

my @hklIncident      = (1,  0,  0);
my @hklPerpendicular = (0,  0,  1);

my $cif = new CIF($CIFFile);
$cif->Read($CIFFile);
my $Crystal = $cif->GetCCrystal();
$Crystal->CalcMetrics();

my @hklHorizontal = $Crystal->GetPerpendicularHKL(@hklIncident, @hklPerpendicular);
print "Incident     : (", join(',', @hklIncident), ")\n";
print "Perpendicular: (", join(',', @hklPerpendicular), ")\n";
print "Horizontal   : (", join(',', @hklHorizontal), ")\n";

my $out = new JFile($OutFile, "w");
if(!$out) {
	print "Error: Can not write to [$OutFile].\n";
	exit;
}
$out->print("h,k,l,2Theta,X,Y,Fr,Fi,F,I\n");
my @VIncident      = $Crystal->GetVectorFromHKL(@hklIncident);
my $L2VIncident    = Sci::InnerProduct(\@VIncident, \@VIncident);
my $LVIncident     = sqrt($L2VIncident);
my @eIncident;
for(my $i = 0 ; $i < 3 ; $i++) {
	$eIncident[$i] = $VIncident[$i] / $LVIncident;
}
my @VPerpendicular = $Crystal->GetVectorFromHKL(@hklPerpendicular);
my @VHorizontal    = $Crystal->GetVectorFromHKL(@hklHorizontal);
my $W = 10.0;
for(my $h = -$hmax ; $h <= $hmax ; $h++) {
for(my $k = -$kmax ; $k <= $kmax ; $k++) {
for(my $l = -$lmax ; $l <= $lmax ; $l++) {
	next if($h == 0 and $k == 0 and $l == 0);

	my @Vhkl    = $Crystal->GetVectorFromHKL($h, $k, $l);
	my $Qinchkl = Sci::CalVectorAngle(\@VIncident, \@Vhkl);
	next if(abs($Qinchkl - 90.0) > $Qcriteria);

	my $dhkl = $Crystal->CalculateHKLInterplanarSpacing($h, $k, $l); # in angstrom
	my $sinQ = $lambda * 10.0 / 2.0 / $dhkl;
	my $Q    = $todeg * asin($sinQ);
	my $Q2   = $Q * 2.0;
	next if($Q2 > $Q2Max);
	
	my ($Fr, $Fi) = $Crystal->Fhkl($h, $k, $l, 1);
	my $F         = sqrt($Fr*$Fr + $Fi*$Fi);
	next if($F < 1.0e-5);

	my $s = $sinQ / $lambda;
	my $T = CrystalObject::TemperatureFactor($Biso, $s);
	my $I = $F * $F * $T * $T;

	my $X = Sci::InnerProduct(\@Vhkl, \@VHorizontal);
	my $Y = Sci::InnerProduct(\@Vhkl, \@VPerpendicular);
	printf "%2d%2d%2d: %7.3f deg. XY:(%7.3f,%7.3f)  (%7.3f,%7.3f) %7.3f\n", 
			$h, $k, $l, $Q2, $X, $Y, $Fr, $Fi, $F;
	$out->print("$h,$k,$l,$Q2,$X,$Y,$Fr,$Fi,$F,$I\n");

}
}
}
$out->Close();
exit;