#!/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  = 0.01e3; # 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/Si5.cif"; #"Si5.cif";
my $OutFile = "LEED.csv";
my $Biso    = 0.5; #1.0;
my $W       = 10.0; #スクリーンサイズ

my $Q2Max = 45.0;

my $hmax = 15;
my $kmax = 15;
my $lmax = 15;

my @VIncident        = (0, 0, -1);
my @hklXPlane        = (1, 0,  0);
my @hklPerpendicular = (0, 0,  1);

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

my @hklYPlane = $Crystal->GetPerpendicularHKL(@hklXPlane, @hklPerpendicular);
print "Incident : (", join(',', @VIncident), ")\n";
print "hklXPlane: (", join(',', @hklPerpendicular), ")\n";
print "hklYPlane: (", join(',', @hklYPlane), ")\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(@VIncident);
my @hklXPlane  = $Crystal->GetVectorFromHKL(@hklXPlane);
my @hklYPlane  = $Crystal->GetVectorFromHKL(@hklYPlane);
my @hklPerpendicular     = $Crystal->GetVectorFromHKL(@hklPerpendicular);


my $peIncident      = Sci::NormalizeVector(\@VIncident);
my $peXPlane        = Sci::NormalizeVector(\@hklXPlane);
my $peYPlane        = Sci::NormalizeVector(\@hklYPlane);
my $pePerpendicular = Sci::NormalizeVector(\@hklPerpendicular);
print "Incident     : ", join(', ', @$peIncident), "\n";
print "XPlane       : ", join(', ', @$peXPlane),    "\n";
print "YPlane       : ", join(', ', @$peYPlane), "\n";
print "Perpendicular: ", join(', ', @$pePerpendicular),      "\n";
my $InnerIncidentPerpendicular = Sci::InnerProduct($peIncident, $pePerpendicular) / $lambda;

my @GhklPlaneX;
my @GhklPlaneY;
my $count = 0;
my $k2 = 1.0 / $lambda / $lambda;
for(my $h = -$hmax ; $h <= $hmax ; $h++) {
for(my $k = -$kmax ; $k <= $kmax ; $k++) {
for(my $l = -$lmax ; $l <= $lmax ; $l++) {
	my ($Fr, $Fi) = (1, 0);
	unless($h == 0 and $k == 0 and $l == 0) {
		($Fr, $Fi) = $Crystal->Fhkl($h, $k, $l, 1);
	}
	my $F         = sqrt($Fr*$Fr + $Fi*$Fi);
	next if($F < 1.0e-5);

	my @Vhkl = $Crystal->GetVectorFromHKL($h, $k, $l);
	my $c1   = Sci::InnerProduct(\@Vhkl, $pePerpendicular);
	my @GhklPlane;
	for(my $i = 0 ; $i < 3 ; $i++) {
		$GhklPlane[$i] = $VIncident[$i] + $Vhkl[$i] - $c1 * $pePerpendicular->[$i];
	}
	my $L2GhklPlane = Sci::InnerProduct(\@GhklPlane, \@GhklPlane);
#print "$h,$k,$l: Inc : ", join(', ', @$pVIncident), "\n";
#print "          Ghkl: ", join(', ', @Vhkl), "\n";
#print "               $k2 - $L2GhklPlane\n";
	my $PlaneX = Sci::InnerProduct(\@GhklPlane, $peXPlane);
	my $PlaneY = Sci::InnerProduct(\@GhklPlane, $peYPlane);
	my $IsNew = 1;
	for(my $i = 0 ; $i < $count ; $i++) {
		my $dx = $PlaneX - $GhklPlaneX[$i];
		my $dy = $PlaneY - $GhklPlaneY[$i];
		my $d2 = $dx*$dx + $dy*$dy;
#print "i=$i: $d2\n";
		if($d2 < 1.0e-10) {
			$IsNew = 0;
			last;
		}
	}
	if(!$IsNew) {
		next;
	}

	my $a2 = $InnerIncidentPerpendicular;
	my $a3 = $L2GhklPlane - $k2;
	my $check = $a2*$a2 - $a3;
	next if($check <= 0.0);

	my $c2 = -$a2 - sqrt($check);
	my @kdif;
	for(my $i = 0 ; $i < 3 ; $i++) {
		$kdif[$i] = $GhklPlane[$i] + $c2 * $pePerpendicular->[$i];
	}
	my $Q2 = Sci::CalVectorAngle(\@VIncident, \@kdif);
#print "          kdif: ", join(', ', @kdif), "\n";
#print "           Q=$Q2\tF=$F\n";
	next unless($Q2 > 90 and $Q2 > 90+$Q2Max);

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

	my $X = Sci::InnerProduct($pVPlane, $peXPlane);
	my $Y = Sci::InnerProduct($pVPlane, $peYPlane);

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

	printf "%2d%2d%2d: %6.2f, 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");

	$GhklPlaneX[$count] = $PlaneX;
	$GhklPlaneY[$count] = $PlaneY;
	$count++;
}
}
}
$out->Close();

exit;