#!/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 $IsPrint = 1;
my $ElectronEnergy  = 0.03e3; # 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 = "RHEED.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 $IncidentAngle    = 3.0;
my @hklIncidentPlane = (1, 0, 0);
my @hklPerpendicular = (0, 1, 1);

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

my @hklHorizontal = $Crystal->GetPerpendicularHKL(@hklIncidentPlane, @hklPerpendicular);
print "Incident     : (", join(',', @hklIncidentPlane), ")\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 @VIncidentPlane  = $Crystal->GetVectorFromHKL(@hklIncidentPlane);
my @VPerpendicular  = $Crystal->GetVectorFromHKL(@hklPerpendicular);
my @VHorizontal     = $Crystal->GetVectorFromHKL(@hklHorizontal);
my $peIncidentPlane = Sci::NormalizeVector(\@VIncidentPlane);
my $pePerpendicular = Sci::NormalizeVector(\@VPerpendicular);
my $peHorizontal    = Sci::NormalizeVector(\@VHorizontal);
my $pVIncident = [];
my $cosQinc    = cos($IncidentAngle * $torad);
my $sinQinc    = sin($IncidentAngle * $torad);
for(my $i = 0 ; $i < 3 ; $i++) {
	$pVIncident->[$i] = $cosQinc * $peIncidentPlane->[$i] - $sinQinc * $pePerpendicular->[$i];
	$pVIncident->[$i] = $pVIncident->[$i] / $lambda; # in nm
}
my $peIncident    = Sci::NormalizeVector($pVIncident);
print "IncPlane     : ", join(', ', @$peIncidentPlane), "\n";
print "Horizontal   : ", join(', ', @$peHorizontal),    "\n";
print "Perpendicular: ", join(', ', @$pePerpendicular), "\n";
print "Incident     : ", join(', ', @$pVIncident),      "\n";

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 @Vhkl = $Crystal->GetVectorFromHKL($h, $k, $l);
my $c = Sci::InnerProduct(\@Vhkl, $peIncidentPlane);
next if($c > 0.0);
	my $c1 = Sci::InnerProduct(\@Vhkl, $pePerpendicular);
	next if($c1 < 0.0);

	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 @GhklPlane;
	for(my $i = 0 ; $i < 3 ; $i++) {
		$GhklPlane[$i] = $pVIncident->[$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, $peHorizontal);
	my $PlaneY = Sci::InnerProduct(\@GhklPlane, $peIncidentPlane);
	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 = -$sinQinc / $lambda;
	my $a3 =  $L2GhklPlane - $k2;
	my $check = $a2*$a2 - $a3;
	next if($check <= 0.0);

	my $c2 = -$a2 + sqrt($check);
#print "           c2=$c2\n";
	my @kdif;
	for(my $i = 0 ; $i < 3 ; $i++) {
		$kdif[$i] = $GhklPlane[$i] + $c2 * $pePerpendicular->[$i];
	}

	my $Q2 = Sci::CalVectorAngle($pVIncident, \@kdif);
#print "          kdif: ", join(', ', @kdif), "\n";
#print "           Q=$Q2\tF=$F\n";
#	next if($Q2 < $Q2Min or $Q2 > $Q2Max);
	next if($Q2 > $Q2Max);

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

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

	my $X = Sci::InnerProduct($pVPlane, $peHorizontal);
	my $Y = Sci::InnerProduct($pVPlane, $pePerpendicular);
	if($IsPrint) {
		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;