#!/usr/bin/perl

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

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

my $CIFFile       = "d:/MyWebs/Research/CrystalStructure/Si5.cif"; #"Si5.cif";
my $OutFile       = "Laue.csv";

my $Biso          =  0.5; #1.0;

my $LTransmission =  5.0; # cm, “§‰ίLaue”z’u‚Ε‚ΜŽŽ—Ώ-ƒXƒNƒŠ[ƒ“‹——£
my $LBack         =  3.0; # cm, ”w–Κ”½ŽΛLaue”z’u‚Ε‚ΜŽŽ—Ώ-ƒXƒNƒŠ[ƒ“‹——£
my $W             =  8.9; # ƒtƒBƒ‹ƒ€•
my $H             = 11.4; # ƒtƒBƒ‹ƒ€‚‚³
my $Q2Min         =  0.1;
my $Q2Max         = 15.0;

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

my $hmax = 20;
my $kmax = 20;
my $lmax = 20;

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,Qinc,X,Y,Fr,Fi,F,I\n");

my @VIncident       = $Crystal->GetVectorFromHKL(@hklIncident);
my @VPerpendicular  = $Crystal->GetVectorFromHKL(@hklPerpendicular);
my @VHorizontal     = $Crystal->GetVectorFromHKL(@hklHorizontal);
my $peIncident      = Sci::NormalizeVector(\@VIncident);
my $pePerpendicular = Sci::NormalizeVector(\@VPerpendicular);
my $peHorizontal    = Sci::NormalizeVector(\@VHorizontal);

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);
	my $Q2      = -180.0 + 2.0 * $Qinchkl;
	my $Q       = 0.5 * $Q2;
	if($TransmissionLaue) {
		next if(abs($Q2) < $Q2Min or abs($Q2) > $Q2Max);
	}
	else {
		my $q2 = 180 - $Q2;
		next if(abs($q2) < $Q2Min or abs($q2) > $Q2Max);
	}

	my $dhkl   = $Crystal->CalculateHKLInterplanarSpacing($h, $k, $l); # in angstrom
	my $sinQ   = sin($Q * $torad);
	my $lambda = abs(2.0 * $dhkl * $sinQ * 0.1); # in nm

	my ($Fr, $Fi) = $Crystal->Fhkl($h, $k, $l);
	my $F         = sqrt($Fr*$Fr + $Fi*$Fi);
	next if($F < 1.0e-5);

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

	my $LVhkl  = sqrt(Sci::InnerProduct(\@Vhkl, \@Vhkl));
	my $peGhkl = Sci::NormalizeVector(\@Vhkl);

	my $Ghkl = 1.0 / $dhkl;            # in angstrom
	my $kinc = 1.0 / ($lambda * 10.0); # in angstrom
	my $pkIncident = Sci::NormalizeVector(\@VIncident, $kinc);
	my $pkhkl      = Sci::NormalizeVector(\@Vhkl,      $Ghkl);
	my $pkdif = [];
	for(my $i = 0 ; $i < 3 ; $i++) {
		$pkdif->[$i] = $pkIncident->[$i] + $pkhkl->[$i];
	}

	my $InnerIncidentDiff = Sci::InnerProduct($peIncident, $pkdif);
	my $pVPlane = [];
	for(my $i = 0 ; $i < 3 ; $i++) {
		$pVPlane->[$i] = $pkdif->[$i] / $InnerIncidentDiff - $peIncident->[$i];
	}

	my $X = Sci::InnerProduct($pVPlane, $peHorizontal);
	my $Y = Sci::InnerProduct($pVPlane, $pePerpendicular);
	printf "%2d%2d%2d: %6.2f, %7.4f nm XY:(%7.3f,%7.3f)  (%7.3f,%7.3f) %7.3f\n", 
			$h, $k, $l, $Q2, $lambda, $X, $Y, $Fr, $Fi, $F;
	$out->print("$h,$k,$l,$Q2,$X,$Y,$Fr,$Fi,$F,$I\n");

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