package Diffraction;
use Common;
@ISA = qw(Common);

#公開したいサブルーチン
#@EXPORT   = qw();
@EXPORT_OK = qw();

use strict;
use Math::MatrixReal;
use Math::Complex;

use Sci qw($h $hbar $me $e $c $todeg $torad);# asin tan);
use Sci::Algorism;
use Sci::Optimize;
use Crystal::Crystal;
use Crystal::CIF;
#use Crystal::Rietan;
use Crystal::RietanFP;

#================================================
# 静的メンバー関数
#================================================
sub ElectroneVTonm
{
	my ($ElectronEnergy, $Relativistic) = @_;
	$Relativistic = 1 if(!defined $Relativistic);

	my $l0 = 2.0 * $me * $ElectronEnergy * $e;
	if($Relativistic) {
		my $K =1.0 + $ElectronEnergy * $e / 2.0 / $me / $c / $c;
		return $h / sqrt($l0 * $K) * 1.0e9; # in nm
	}
	return $h / sqrt($l0) * 1.0e9; # in nm
#	return sqrt($h * $h / 2.0 / $me / $ElectronEnergy / $e) * 1.0e9; # in nm
}

#==================================================================
# 変数取得関数
#==================================================================
sub TransmissionMode   { return shift->{TransmissionMode};                  }
sub Wavelength         { return shift->{Wavelength};                        }
sub SetWavelength      { my ($this,$wl)=@_; return $this->{Wavelength}=$wl; }
sub SourceSpectrum     { return shift->{pSourceSpectrum};                        }
sub SetSourceSpectrum  { my ($this,$p)=@_; return $this->{pSourceSpectrum}=$p; }
sub DiffractionMode    { return shift->{DiffractionMode};                   }
sub SetDiffractionMode 
{ 
	my ($this, $mode) = @_; 
	$this->{DiffractionMode} = $mode;
	if($mode =~ /Back/i) {
		$this->{TransmissionMode} = 0 
	}
	else {
		$this->{TransmissionMode} = 1;
	}
	return $this->{DiffractionMode} = $mode;
}
sub Source    { return shift->{Source};                        }
sub SetSource { my ($this,$s)=@_; return $this->{Source} = $s; }
sub Crystal   { return shift->{pCrystal};                      }
sub SetCrystal
{
	my ($this, $pCrystal, $Source, $StopByError) = @_;
	$StopByError = 1 if(!defined $StopByError);
	$Source = $this->{Source} if(!defined $Source);
	$this->{pCrystal} = $pCrystal;
	$pCrystal->ReadASF($Source, $StopByError);
	return $this->{pCrystal};
}
sub SetOrienation
{
	my ($this, $h, $k, $l, $shape, $p1, $p2) = @_;
#print "SO: $this, $h, $k, $l, $shape, $p1, $p2\n";
#exit;
	$this->{Orientationh} = $h;
	$this->{Orientationk} = $k;
	$this->{Orientationl} = $l;
	$this->{GrainShape}   = $shape;
	$this->{Orientationp1} = $p1;
	$this->{Orientationp2} = $p2;
}

sub hklRange { my ($this)=@_; return ($this->{hmax}, $this->{kmax}, $this->{lmax}); }
sub SethklRange
{
	my $this = shift;
	return ($this->{hmax}, $this->{kmax}, $this->{lmax}) = @_;
}
sub Q2Range { my ($this)=@_; return ($this->{Q2Start}, $this->{Q2Stop}, $this->{Q2Step}, $this->{nQ2}); }
sub SetQ2Range {
	my ($this, $Q2Start, $Q2Stop, $Q2Step) = @_;
	$this->{Q2Start} = $Q2Start;
	$this->{Q2Stop}  = $Q2Stop;
	$this->{Q2Step}  = $Q2Step;
	$this->{nQ2}     = int( ($Q2Stop - $Q2Start) / $Q2Step + 1.001);
	return ($Q2Start, $Q2Stop, $Q2Step, $this->{nQ2});
}

sub Incidenthkl { my ($this)=@_; return ($this->{Incidenth}, $this->{Incidentk}, $this->{Incidentl}); }
sub SetIncidenthkl
{
	my $this = shift;
	return ($this->{Incidenth}, $this->{Incidentk}, $this->{Incidentl}) = @_;
}
sub Perpendicularhkl { my ($this)=@_; return ($this->{Perpendicularh}, $this->{Perpendiculark}, $this->{Perpendicularl}); }
sub SetPerpendicularhkl
{
	my $this = shift;
	return ($this->{Perpendicularh}, $this->{Perpendiculark}, $this->{Perpendicularl}) = @_;
}
sub Horizontalhkl { my ($this)=@_; return ($this->{Horizontalh}, $this->{Horizontalk}, $this->{Horizontall}); }
sub SetHorizontalhkl
{
	my $this = shift;
	return ($this->{Horizontalh}, $this->{Horizontalk}, $this->{Horizontall}) = @_;
}
sub ScreenGeometry { my ($this)=@_; return ($this->{ScreenW}, $this->{ScreenH}, $this->{ScreenDistance}); }
sub SetScreenGeometry
{
	my ($this, $w, $h, $distance) = @_;
	$this->{ScreenW}        = $w;
	$this->{ScreenH}        = $h;
	$this->{ScreenDistance} = $distance;
}
sub PeakWidth              { return shift->{PeakWidth}; }
sub SetPeakWidth           { my ($this,$w)=@_; return $this->{PeakWidth} = $w; }
sub ThermalVibrationFactor { my ($this)=@_; return $this->{Biso}; }
sub SetThermalVibrationFactor
{
	my ($this, $Biso) = @_;
	return $this->{Biso} = $Biso;
}
sub IncidentAngle    { return shift->{IncidentAngle};                           }
sub SetIncidentAngle { my ($this, $Q) = @_; return $this->{IncidentAngle} = $Q; }
sub Q2Min            { return shift->{Q2Min};                                   }
sub SetQ2Min         { my ($this, $Q2Min) = @_; return $this->{Q2Min} = $Q2Min; }
sub Q2Max            { return shift->{Q2Max};                                   }
sub SetQ2Max         { my ($this, $Q2Max) = @_; return $this->{Q2Max} = $Q2Max; }

sub pXRDQ2           { return shift->{pXRDQ2}; }
sub SetpXRDQ2        { my ($this,$p)=@_; return $this->{pXRDQ2} = $p; }
sub pXRDIntensity    { return shift->{pXRDIntensity}; }
sub SetpXRDIntensity { my ($this,$p)=@_; return $this->{pXRDIntensity} = $p; }

sub pQ2                 { return shift->{pQ2}; }
sub SetpQ2              { my ($this,$p)=@_; return $this->{pQ2} = $p; }
sub ppFArray            { return shift->{ppF}; }
sub SetppFArray         { my ($this,$p)=@_; return $this->{ppF} = $p; }
sub pIntensityArray     { return shift->{pIntensity}; }
sub SetpIntensityArray  { my ($this,$p)=@_; return $this->{pIntensity} = $p; }
sub prIntensityArray    { return shift->{prI}; }
sub SetprIntensityArray { my ($this,$p)=@_; return $this->{prI} = $p; }
sub pScreenXArray       { return shift->{pScreenXArray};                        }
sub SetpScreenXArray    { my ($this,$p)=@_; return $this->{pScreenXArray} = $p; }
sub pScreenYArray       { return shift->{pScreenYArray};                        }
sub SetpScreenYArray    { my ($this,$p)=@_; return $this->{pScreenYArray} = $p; }
sub pTagArray           { return shift->{pTag}; }
sub SetpTagArray        { my ($this,$p)=@_; return $this->{pTag} = $p; }
sub phklArray           { return shift->{pphkl}; }
sub pWLArray            { return shift->{pWL}; }
sub pQ2Array            { return shift->{pQ2}; }
sub pISourceArray       { return shift->{pISource}; }

#==================================================================
# コンストラクタ、デストラクト
#==================================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

#==================================================================
# 一般メンバー関数
#==================================================================
sub Prepare
{
	my ($this, $DiffractionMode) = @_;
	my $Crystal     = $this->Crystal();
	my ($W, $H, $D) = $this->ScreenGeometry();

	$this->SetDiffractionMode($DiffractionMode) if(defined $DiffractionMode);
	my $ScreenL;
	if($this->DiffractionMode() =~ /Powder.*XRD/i) {
	}
	elsif($this->DiffractionMode() =~ /(Laue|TED)/i) {
		$ScreenL  = sqrt($W*$W + $H*$H);
		$this->SetQ2Max($todeg * atan($ScreenL / 2.0 / $D));
		$this->SetHorizontalhkl( $Crystal->GetPerpendicularHKL($this->Incidenthkl(), $this->Perpendicularhkl()) );
	}
	else {
		$ScreenL  = sqrt($W*$W + $H*$H);
		$this->SetQ2Max($todeg * atan($ScreenL / 2.0 / $D));
		$this->SetHorizontalhkl( $Crystal->GetPerpendicularHKL($this->Incidenthkl(), $this->Perpendicularhkl()) );
	}
}

sub ResetFhkl
{
	my ($this, $f) = @_;
	return if(!$f);

	delete $this->{pData};
	delete $this->{pphkl};
	delete $this->{pWL};
	delete $this->{pQ2};
	delete $this->{ppF};
	delete $this->{pIntensity};
	delete $this->{prI};
	delete $this->{pISource};
	delete $this->{pScreenXArray};
	delete $this->{pScreenYArray};
}

sub pData { return shift->{pData}; }

sub SetpData
{
	my ($this, $pData) = @_;

	my $IsLaue = 0;
	$IsLaue = 1 if($this->{DiffractionMode} =~ /Laue/i);
	my $n = @$pData;
	my $Imax = 0.0;
	for(my $i = 0 ; $i < $n ; $i++) {
		my $p = $pData->[$i];
		my ($h, $k, $l, $lambda, $Q2, $X, $Y, $Fr, $Fi, $F, $P, $T, $I, $ISource) = @$p;
		$Imax = $I if($I > $Imax);
	}
	my (@phkl, @WL, @Q2, @pF, @I, @rI, @ISource, @X, @Y, @Tag);
	for(my $i = 0 ; $i < $n ; $i++) {
		my $p = $pData->[$i];
		my ($h, $k, $l, $lambda, $Q2, $X, $Y, $Fr, $Fi, $F, $P, $T, $I, $ISource) = @$p;
		my $rI = $I / $Imax;
		my $tag;
		if($IsLaue) {
			$tag = sprintf("%d %d %d: l=%6.4f nm Q2=%6.4f\nF=%6.2f I=%6.2f Isrc=%6.4f", 
								$h, $k, $l, $lambda, $Q2, $F, $rI*1000.0, $ISource);
		}
		else {
			$tag = sprintf("%d %d %d: Q2=%6.4f\nF=%6.2f I=%6.2f", $h, $k, $l, $Q2, $F, $rI*1000.0);
		}
		$phkl[$i]    = [$h, $k, $l];
		$WL[$i]      = $lambda;
		$Q2[$i]      = $Q2;
		$pF[$i]      = [$Fr, $Fi, $F];
		$I[$i]       = $I;
		$rI[$i]      = $rI;
		$ISource[$i] = $ISource;
		$X[$i]       = $X;
		$Y[$i]       = $Y;
		$Tag[$i]     = $tag;
	}
	$this->{pData} = $pData;
	$this->{pphkl} = \@phkl;
	$this->{pWL}   = \@WL;
	$this->{pQ2}   = \@Q2;
	$this->SetppFArray(\@pF);
	$this->SetpIntensityArray(\@I);
	$this->SetprIntensityArray(\@rI);
	$this->{pISource} = \@ISource;
	$this->SetpScreenXArray(\@X);
	$this->SetpScreenYArray(\@Y);
	$this->SetpTagArray(\@Tag);
}

sub CalLauePattern
{
	my ($this, $callback) = @_;
	my $Crystal              = $this->Crystal();
	my ($Biso)               = $this->ThermalVibrationFactor();
	my $Q2Min                = $this->Q2Min();
	my $Q2Max                = $this->Q2Max();
	my ($hmax, $kmax, $lmax) = $this->hklRange();
	my ($W, $H, $D)          = $this->ScreenGeometry();
	my $TransmissionLaue     = $this->TransmissionMode();
	my $pSourceSpectrum      = $this->SourceSpectrum();
#print "h=($hmax, $kmax, $lmax)\n";
#print "w=($W, $H, $D)\n";
#print "t=$TransmissionLaue\n";
#print "Q=$Q2Min,$Q2Max\n";

	$this->Prepare("Laue");
	my @VIncident       = $Crystal->GetVectorFromHKL($this->Incidenthkl());
	my @VPerpendicular  = $Crystal->GetVectorFromHKL($this->Perpendicularhkl());
	my @VHorizontal     = $Crystal->GetVectorFromHKL($this->Horizontalhkl());
printf "Inc: %6.3f %6.3f %6.3f\n", @VIncident;
printf "Per: %6.3f %6.3f %6.3f\n", @VPerpendicular;
printf "Hor: %6.3f %6.3f %6.3f\n", @VHorizontal;
	my $peIncident      = Sci::NormalizeVector(\@VIncident);
	my $pePerpendicular = Sci::NormalizeVector(\@VPerpendicular);
	my $peHorizontal    = Sci::NormalizeVector(\@VHorizontal);

	my (@LaueX, @LaueY, @pData);
	my $c = 0;
	my $DispCount = 0;
	for(my $h = -$hmax ; $h <= $hmax ; $h++) {
	for(my $k = -$kmax ; $k <= $kmax ; $k++) {
	for(my $l = -$lmax ; $l <= $lmax ; $l++) {
		&$callback() if($callback);
		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 $ISource = 1.0;
		if($pSourceSpectrum) {
			$ISource = $pSourceSpectrum->YVal($lambda * 10.0) * 100.0;
		}
		my $I = $ISource * $F * $F * $P * $T * $T;
if($F < 1.0e-5) {
	$Fr = 0.0; 	$Fi = 0.0; 	$F  = 0.0; 	$I  = 0.0;
}

#		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];
			$pkdif->[$i] = $pkIncident->[$i] + $Vhkl[$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 = $D * Sci::InnerProduct($pkdif, $peHorizontal) / $InnerIncidentDiff;
#		my $Y = $D * Sci::InnerProduct($pkdif, $pePerpendicular) / $InnerIncidentDiff;
		my $X = $D * Sci::InnerProduct($pVPlane, $peHorizontal);
		my $Y = $D * Sci::InnerProduct($pVPlane, $pePerpendicular);
if($DispCount % 100 == 0) {
	printf "%2d %2d %2d: %5.2f,%7.4fnm (%6.3f,%6.3f)  (%6.3f,%6.3f) %6.3f (%6.3f)\n", 
			$h, $k, $l, $Q2, $lambda, $X, $Y, $Fr, $Fi, $F, $ISource;
}
$DispCount++;
		$LaueX[$c] = $X;
		$LaueY[$c] = $Y;
		$pData[$c] = [$h, $k, $l, $lambda, $Q2, $X, $Y, $Fr, $Fi, $F, $P, $T, $I, $ISource];
		$c++;
	}
	}
	}
	$this->SetpData(\@pData);
print "Diffraction::CalLauePattern: Calculation finished\n";
}

sub CalTEDPattern
{
	my ($this, $callback) = @_;
	my $Crystal              = $this->Crystal();
	my ($Biso)               = $this->ThermalVibrationFactor();
	my $Qcriteria            = $this->Q2Min();
	my $Q2Max                = $this->Q2Max();
	my ($hmax, $kmax, $lmax) = $this->hklRange();
	my ($W, $H, $D)          = $this->ScreenGeometry();
	my $lambda               = $this->Wavelength();

	$this->Prepare("TED");
	my @VIncident      = $Crystal->GetVectorFromHKL($this->Incidenthkl());
	my @VPerpendicular  = $Crystal->GetVectorFromHKL($this->Perpendicularhkl());
#	my @VHorizontal     = Sci::VectorProduct(\@VIncident, \@VPerpendicular);
	my @VHorizontal     = $Crystal->GetVectorFromHKL($this->Horizontalhkl());
printf "Inc: %6.3f %6.3f %6.3f\n", @VIncident;
printf "Per: %6.3f %6.3f %6.3f\n", @VPerpendicular;
printf "Hor: %6.3f %6.3f %6.3f\n", @VHorizontal;
	my $peIncident      = Sci::NormalizeVector(\@VIncident);
	my $pePerpendicular = Sci::NormalizeVector(\@VPerpendicular);
	my $peHorizontal    = Sci::NormalizeVector(\@VHorizontal);
	my (@TEDX, @TEDY, @pData);
	my $c = 0;
	for(my $h = -$hmax ; $h <= $hmax ; $h++) {
	for(my $k = -$kmax ; $k <= $kmax ; $k++) {
	for(my $l = -$lmax ; $l <= $lmax ; $l++) {
		&$callback() if($callback);
		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 $P = 1.0;
		my $T = CrystalObject::TemperatureFactor($Biso, $s);
		my $I = $F * $F * $P * $T * $T;
if($F < 1.0e-5) {
	$Fr = 0.0;
	$Fi = 0.0;
	$F  = 0.0;
	$I  = 0.0;
}
printf "%2d%2d%2d: Vhkl ((%7.3f,%7.3f,%7.3f)\n", $h, $k, $l, @Vhkl;
		my $X = Sci::InnerProduct(\@Vhkl, $peHorizontal);
		my $Y = Sci::InnerProduct(\@Vhkl, $pePerpendicular);
#		my $X = $D * Sci::InnerProduct(\@Vhkl, \@VHorizontal);
#		my $Y = $D * Sci::InnerProduct(\@Vhkl, \@VPerpendicular);
		my $K = $D / sqrt($X*$X + $Y*$Y) * tan($Q2 * $torad);
		$X *= $K;
		$Y *= $K;
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;
		$TEDX[$c] = $X;
		$TEDY[$c] = $Y;
		$pData[$c] = [$h, $k, $l, $lambda, $Q2, $X, $Y, $Fr, $Fi, $F, $P, $T, $I];
		$c++;
	}
	}
	}
	$this->SetpData(\@pData);
print "Diffraction::CalTEDPattern: Calculation finished\n";
}

sub CalRHEEDPattern
{
	my ($this, $callback, $StopByError) = @_;
	$StopByError = 1 if(!defined $StopByError);
	my $Crystal              = $this->Crystal();
	my ($Biso)               = $this->ThermalVibrationFactor();
	my $Qcriteria            = $this->Q2Min();
	my $Q2Max                = $this->Q2Max();
	my ($hmax, $kmax, $lmax) = $this->hklRange();
	my ($W, $H, $D)          = $this->ScreenGeometry();
	my $lambda               = $this->Wavelength() * 10.0; # in A
	my $IncidentAngle        = $this->IncidentAngle();

	$this->Prepare("RHEED");
	my @hklIncidentPlane = $this->Incidenthkl();
	my @hklPerpendicular = $this->Perpendicularhkl();
	my @hklHorizontal    = $this->Horizontalhkl();

	my @VIncidentPlane  = $Crystal->GetVectorFromHKL(@hklIncidentPlane); # in A-1
	my @VPerpendicular  = $Crystal->GetVectorFromHKL(@hklPerpendicular); # in A-1
	my @VHorizontal     = $Crystal->GetVectorFromHKL(@hklHorizontal);    # in A-1
	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 A-1
	}
	my $peIncident    = Sci::NormalizeVector($pVIncident);
printf "IncPlane     : %6.3f %6.3f %6.3f\n", @$peIncidentPlane;
printf "Horizontal   : %6.3f %6.3f %6.3f\n", @$peHorizontal;
printf "Perpendicular: %6.3f %6.3f %6.3f\n", @$pePerpendicular;
printf "Incident     : %6.3f %6.3f %6.3f\n", @$pVIncident;

	my $k2 = 1.0 / $lambda / $lambda; # in A-2
#	my (@GhklPlaneX, @GhklPlaneY);
	my (@RHEEDX, @RHEEDY, @pData);
	my $DiffractionCount = 0;
	my $count = 0;
	for(my $h = -$hmax ; $h <= $hmax ; $h++) {
print "h=$h\n";
	for(my $k = -$kmax ; $k <= $kmax ; $k++) {
	for(my $l = -$lmax ; $l <= $lmax ; $l++) {
		&$callback() if($callback);
		my @Vhkl = $Crystal->GetVectorFromHKL($h, $k, $l); # in A-1
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, $StopByError);
		}
		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);
#		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;
#			if($d2 < 1.0e-10) {
#				$IsNew = 0;
#				last;
#			}
#		}
#		if(!$IsNew) {
#			next;
#		}

		my $a2 = -$sinQinc / $lambda; # in A-1
		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($pVIncident, \@kdif);
		next if($Q2 > $Q2Max);

		my $sinQ = sin(0.5 * $Q2 * $torad);
		my $s = $sinQ / $lambda; # in A-1
		my $P = 1.0;
		my $T = CrystalObject::TemperatureFactor($Biso, $s);
		my $I = $F * $F * $P * $T * $T;
if($F < 1.0e-5) {
	$Fr = 0.0;
	$Fi = 0.0;
	$F  = 0.0;
	$I  = 0.0;
}

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

		my $X = $D * Sci::InnerProduct($pVPlane, $peHorizontal);
		my $Y = $D * Sci::InnerProduct($pVPlane, $pePerpendicular);
if($DiffractionCount % 100 == 0) {
	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;
}
$DiffractionCount++;
#		$GhklPlaneX[$count] = $PlaneX;
#		$GhklPlaneY[$count] = $PlaneY;
		$RHEEDX[$count]     = $X;
		$RHEEDY[$count]     = $Y;
		$pData[$count] = [$h, $k, $l, $lambda*0.1, $Q2, $X, $Y, $Fr, $Fi, $F, $P, $T, $I];
		$count++;
	}
	}
	}
	$this->SetpData(\@pData);
print "Diffraction::CalRHEEDPattern: Calculation finished\n";
}

sub CalLEEDPattern
{
	my ($this, $callback) = @_;
	my $Crystal              = $this->Crystal();
	my ($Biso)               = $this->ThermalVibrationFactor();
	my $Qcriteria            = $this->Q2Min();
	my $Q2Max                = $this->Q2Max();
	my ($hmax, $kmax, $lmax) = $this->hklRange();
	my ($W, $H, $D)          = $this->ScreenGeometry();
	my $lambda               = $this->Wavelength() * 10.0; # in A

	$this->Prepare("LEED");
	my @hklIncidentPlane = $this->Incidenthkl();
	my @hklPerpendicular = $this->Perpendicularhkl();
	my @hklHorizontal    = $this->Horizontalhkl();

	my @VIncidentPlane  = $Crystal->GetVectorFromHKL(@hklIncidentPlane); # in A-1
	my @VPerpendicular  = $Crystal->GetVectorFromHKL(@hklPerpendicular); # in A-1
	my @VHorizontal     = $Crystal->GetVectorFromHKL(@hklHorizontal);    # in A-1
	my $peIncidentPlane = Sci::NormalizeVector(\@VIncidentPlane);
	my $pePerpendicular = Sci::NormalizeVector(\@VPerpendicular);
	my $peHorizontal    = Sci::NormalizeVector(\@VHorizontal);
printf "IncPlane     : %6.3f %6.3f %6.3f\n", @$peIncidentPlane;
printf "Horizontal   : %6.3f %6.3f %6.3f\n", @$peHorizontal;
printf "Perpendicular: %6.3f %6.3f %6.3f\n", @$pePerpendicular;

	my $k2 = 1.0 / $lambda / $lambda; # in A-2
	my (@LEEDX, @LEEDY, @pData);
	my $count = 0;
	my $DiffractionCount = 0;
	for(my $h = -$hmax ; $h <= $hmax ; $h++) {
print "h=$h\n";
	for(my $k = -$kmax ; $k <= $kmax ; $k++) {
	for(my $l = -$lmax ; $l <= $lmax ; $l++) {
		next if($h == 0 and $k == 0 and $l == 0);
		&$callback() if($callback);

		my @Vhkl = $Crystal->GetVectorFromHKL($h, $k, $l); # in A-1
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) = $Crystal->Fhkl($h, $k, $l, 1);
		my $F = sqrt($Fr*$Fr + $Fi*$Fi);
#		next if($F < 1.0e-5);

		my $c1 = Sci::InnerProduct(\@Vhkl, $peIncidentPlane);
		my @GhklPlane;
		for(my $i = 0 ; $i < 3 ; $i++) {
			$GhklPlane[$i] = $VIncidentPlane[$i] + $Vhkl[$i] - $c1 * $peIncidentPlane->[$i];
		}
		my $L2GhklPlane = Sci::InnerProduct(\@GhklPlane, \@GhklPlane);
		my $PlaneX = Sci::InnerProduct(\@GhklPlane, $peHorizontal);
		my $PlaneY = Sci::InnerProduct(\@GhklPlane, $pePerpendicular);
		my $Theta = atan2($PlaneX, $PlaneY) * $todeg;
		my $PlaneR = sqrt($PlaneX*$PlaneX + $PlaneY*$PlaneY);
		my $sinPsi = $PlaneR * $lambda; # by A-1 * A
		next if(abs($sinPsi) >= 1.0);
		my $Psi   = asin($sinPsi) * $todeg;
	
		my $a2 = -1.0 / $lambda; # in A-1
		my $a3 =  $L2GhklPlane - $k2;
		my $check = $a2*$a2 - $a3;
		next if($check <= 0.0);

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

		my $Q2 = Sci::CalVectorAngle(\@VIncidentPlane, \@kdif);
		next if($Q2 < 90.0); #$Q2 > $Q2Max);

		my $sinQ = sin(0.5 * $Q2 * $torad);
		my $s = $sinQ / $lambda; # in A-1
		my $P = 1.0;
		my $T = CrystalObject::TemperatureFactor($Biso, $s);
		my $I = $F * $F * $P * $T * $T;
if($F < 1.0e-5) {
	$Fr = 0.0;
	$Fi = 0.0;
	$F  = 0.0;
	$I  = 0.0;
}
#my $Theta = atan2($PlaneX, $PlaneY) * $todeg;
#cos(Theta) = $PlaneX / $PlaneR;
		my $LScreen = $D * $Psi * $torad; #$sinPsi;
		my $X = $LScreen * $PlaneX / $PlaneR;
		my $Y = $LScreen * $PlaneY / $PlaneR;
if($DiffractionCount % 100 == 0) {
	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;
}
$DiffractionCount++;
		$LEEDX[$count] = $X;
		$LEEDY[$count] = $Y;
		$pData[$count] = [$h, $k, $l, $lambda*0.1, $Q2, $X, $Y, $Fr, $Fi, $F, $P, $T, $I];
		$count++;
	}
	}
	}
	$this->SetpData(\@pData);
print "Diffraction::CalLEEDPattern: Calculation finished\n";
}

sub CalPowderXRDPattern
{
	my ($this, $callback) = @_;
#return $this->CalPowderXRDFhkl($callback);

	if(!$this->pData()) { #!$this->ppFArray()) {
		$this->CalPowderXRDFhkl($callback);
	}

	my $Crystal              = $this->Crystal();
	my ($Biso)               = $this->ThermalVibrationFactor();
	my $lambda               = $this->Wavelength();

#print "p2: $this->{Orientationp2}\n";
	my $Orietantionh  = (!defined $this->{Orientationh})?  0   : $this->{Orientationh};
	my $Orietantionk  = (!defined $this->{Orientationk})?  0   : $this->{Orientationk};
	my $Orietantionl  = (!defined $this->{Orientationl})?  1   : $this->{Orientationl};
	my $IsSheet       = ($this->{GrainShape} eq 'sheet')?  1   : 0;
	my $Orietantionp1 = (!defined $this->{Orientationp1})? 1.0 : $this->{Orientationp1};
	my $Orietantionp2 = (!defined $this->{Orientationp2})? 0.0 : $this->{Orientationp2};
#print "p2: $Orietantionp1, $Orietantionp2\n";

	my $w                    = $this->PeakWidth();
	my ($Q2Start, $Q2Stop, $Q2Step, $nQ2) = $this->Q2Range();
	my $nQ2w    = int( 20.0 * $w / $Q2Step);
	my $pData = $this->pData();
	my $n     = @$pData;
	my (@XRDX, @XRDY);
	for(my $i = 0 ; $i < $nQ2 ; $i++) {
		$XRDX[$i] = $Q2Start + $i * $Q2Step;
		$XRDY[$i] = 0.0;
	}
	for(my $i = 0 ; $i < $n ; $i++) {
		my $pd = $pData->[$i];
		my ($h, $k, $l, $lambda, $QB2, $ScreenX, $ScreenY, $Fr, $Fi, $F, $P, $T, $I, $PO) = @$pd;
		my $QB    = $QB2 / 2.0;
		my $sinQB = sin($QB * $torad);
		my $s     = $sinQB / $lambda;
		my $L  = CrystalObject::LorentzFactor($QB);
		$T  = CrystalObject::TemperatureFactor($Biso, $s);
		$PO = $Crystal->PreferentialOrientation(
						$IsSheet, $Orietantionh, $Orietantionk, $Orietantionl,
						$Orietantionp1, $Orietantionp2, $h, $k, $l);
#print "PO=$PO ($h,$k,$l): $IsSheet $Orietantionh,$Orietantionk,$Orietantionl,  $Orietantionp1,$Orietantionp2\n";
		$I = $F * $F * $P * $L * $PO * $T * $T;
		&$callback() if($callback);

		my $iQ0 = int( ($QB2 - $Q2Start) / $Q2Step + 0.5);
		for(my $j = $iQ0 - $nQ2w ; $j <= $iQ0 + $nQ2w ; $j++) {
			next if($j < 0 or $j >= $nQ2);

			my $dQ = $Q2Step * ($iQ0 - $j) / $w;
			my $Lorentz = 1.0 / (1.0 + $dQ * $dQ);
			$XRDY[$j] += $I * $Lorentz;
		}
	}
	$this->SetpXRDQ2(\@XRDX);
	$this->SetpXRDIntensity(\@XRDY);
}

sub CalPowderXRDFhkl
{
	my ($this, $callback) = @_;
	my $Crystal              = $this->Crystal();
	my ($Biso)               = $this->ThermalVibrationFactor();
	my ($hmax, $kmax, $lmax) = $this->hklRange();
	my $lambda               = $this->Wavelength();

	my ($Q2Start, $Q2Stop, $Q2Step, $nQ2) = $this->Q2Range();
#	my $w                    = $this->PeakWidth();
#	my $nQ2w    = int( 20.0 * $w / $Q2Step);

	my $Orietantionh  = (!defined $this->{Orientationh})?  0   : $this->{Orientationh};
	my $Orietantionk  = (!defined $this->{Orientationk})?  0   : $this->{Orientationk};
	my $Orietantionl  = (!defined $this->{Orientationl})?  1   : $this->{Orientationl};
	my $IsSheet       = ($this->{GrainShape} eq 'sheet')?  1   : 0;
	my $Orietantionp1 = (!defined $this->{Orientationp1})? 1.0 : $this->{Orientationp1};
	my $Orietantionp2 = (!defined $this->{Orientationp2})? 0.0 : $this->{Orientationp2};
#print "p2: $Orietantionp1, $Orietantionp2\n";

	$this->Prepare("PowderXRD");
	my @pData;
	my (@XRDX, @XRDY);
	for(my $i = 0 ; $i < $nQ2 ; $i++) {
		$XRDX[$i] = $Q2Start + $i * $Q2Step;
		$XRDY[$i] = 0.0;
	}
	my $c = 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);
		&$callback() if($callback);
#print "$h,$k,$l\n";
		my $dhkl      = $Crystal->CalculateHKLInterplanarSpacing($h, $k, $l); # in angstrom
		my $sinQB     = $lambda * 10.0 / 2.0 / $dhkl;
		next if($sinQB >= 1.0);
		my $QB        = asin($sinQB) * $todeg;
		my $QB2       = $QB + $QB;
		next if($QB2 < $Q2Start or $QB2 > $Q2Stop);

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

#print "l=$lambda\n";
		my $s  = $sinQB / $lambda;
		my $P  = CrystalObject::PolarizationFactor($QB);
		my $L  = CrystalObject::LorentzFactor($QB);
		my $T  = CrystalObject::TemperatureFactor($Biso, $s);
		my $PO = $Crystal->PreferentialOrientation(
						$IsSheet, $Orietantionh, $Orietantionk, $Orietantionl,
						$Orietantionp1, $Orietantionp2, $h, $k, $l);
		my $I = $F * $F * $P * $L * $PO * $T * $T;
#		my $I = $F * $F * $m * $P * $PO * $L * $T * $T;
if($F < 1.0e-5) {
	$Fr = 0.0;
	$Fi = 0.0;
	$F  = 0.0;
	$I  = 0.0;
}
printf "F(%d %d %d)=(%7.3f,%7.3f) %7.3f [m=%d] d=%7.3f\n", $h, $k, $l, $Fr, $Fi, $F, $m, $dhkl;
printf "   2QB=%7.3f P=%7.3f L=%7.3f T=%7.3f PO=%6.3f  I=%10.3f\n", 2.0*$QB, $P, $L, $T, $PO, $I;

#		my $iQ0 = int( ($QB2 - $Q2Start) / $Q2Step + 0.5);
#		for(my $j = $iQ0 - $nQ2w ; $j <= $iQ0 + $nQ2w ; $j++) {
#			next if($j < 0 or $j >= $nQ2);
#
#			my $dQ = $Q2Step * ($iQ0 - $j) / $w;
#			my $Lorentz = 1.0 / (1.0 + $dQ * $dQ);
#			$XRDY[$j] += $I * $Lorentz;
#		}
		$pData[$c] = [$h, $k, $l, $lambda, $QB2, 0.0, 0.0, $Fr, $Fi, $F, $P, $T, $I, $PO];
		$c++;
	}
	}
	}
	@pData = sort { int($a->[4]*10000) <=> int($b->[4]*10000) or $a->[9] <=> $b->[9] } @pData;
#for(my $i = 0 ; $i < $nQ2 ; $i++) {
#	my $Q2 = $Q2Start + $i * $Q2Step;
#	my $Q  = $Q2 * 0.5;
#	my $s = sin($torad * $Q) / $lambda;
#	my $P = CrystalObject::PolarizationFactor($Q);
#	my $L = CrystalObject::LorentzFactor($Q);
#	my $T = CrystalObject::TemperatureFactor($Biso, $s);
#	$out->print("$Q2,$XRD[$i],$P,$L,$T\n");
#}

	$this->SetpData(\@pData);
#	$this->SetpXRDQ2(\@XRDX);
#	$this->SetpXRDIntensity(\@XRDY);
print "Diffraction::CalPowderXRDPattern: Calculation finished\n";
}

sub CalGmin {
	my ($this, $pCrystal, $pVars, $iPrintLevel) = @_;
#print "t: ($this, $pCrystal, $pVars, $iPrintLevel)\n";
	my ($h, $k, $l) = @$pVars;
	my $dhkl  = $pCrystal->CalculateHKLInterplanarSpacing($h, $k, $l); # in angstrom
	return 1.0 / $dhkl;
}

sub Guesshklmax
{
	my ($this) = @_;

	my $lambda   = $this->Wavelength() * 10; # in A
	my $pCrystal = $this->Crystal();
	my $Q2max    = $this->{Q2Stop};
	my $sinQmax  = sin($torad * $Q2max / 2.0);
	my ($hmax, $kmax, $lmax);

	my $Opt = new Optimize;
	my $Method = "Amoeba::Simplex";
	my $EPS = 0.01;
	my $nMaxIter = 100;
	my @Scale = (2, 2, 2);
	my $iPrintLevel = 0;

	my @OptId = (0, 1, 1);
	for($hmax = 1 ;  ; $hmax++) {
		my @Guess = ($hmax, 0, 0);
		my ($pOptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
					$EPS, $nMaxIter, $iPrintLevel,
					sub { $this->CalGmin($pCrystal, @_); }, sub {}, sub {},
				);
		my $sinQB = $lambda / 2.0 * $MinVal;
#print "h=$hmax: $sinQB ($sinQmax); (",join(',',@{$pOptVars}),")\n";
		if($sinQB >= $sinQmax) {
			--$hmax;
printf "hmax=%d, Gmin at (%d %3.1f %3.1f) with sinQ=%5.3f (lim=%5.3f)\n", $hmax, @{$pOptVars}, $sinQB, $sinQmax;
			last;
		}
	}
	@OptId = (1, 0, 1);
	for($kmax = 1 ;  ; $kmax++) {
		my @Guess = (0, $kmax, 0);
		my ($pOptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
					$EPS, $nMaxIter, $iPrintLevel,
					sub { $this->CalGmin($pCrystal, @_); }, sub {}, sub {},
				);
		my $sinQB = $lambda / 2.0 * $MinVal;
#print "k=$kmax: $sinQB ($sinQmax); (",join(',',@{$pOptVars}),")\n";
		if($sinQB >= $sinQmax) {
			--$kmax;
printf "kmax=%d, Gmin at (%3.1f %d %3.1f) with sinQ=%5.3f (lim=%5.3f)\n", $kmax, @{$pOptVars}, $sinQB, $sinQmax;
			last;
		}
	}
	@OptId = (1, 1, 0);
	for($lmax = 1 ; ; $lmax++) {
		my @Guess = (0, 0, $lmax);
		my ($pOptVars, $MinVal) = $Opt->Optimize($Method, \@Guess, \@OptId, \@Scale, 
					$EPS, $nMaxIter, $iPrintLevel,
					sub { $this->CalGmin($pCrystal, @_); }, sub {}, sub {},
				);
		my $sinQB = $lambda / 2.0 * $MinVal;
#print "l=$lmax: $sinQB ($sinQmax); (",join(',',@{$pOptVars}),")\n";
		if($sinQB >= $sinQmax) {
			--$lmax;
printf "lmax=%d, Gmin at (%3.1f %3.1f %d) with sinQ=%5.3f (lim=%5.3f)\n", $lmax, @{$pOptVars}, $sinQB, $sinQmax;
			last;
		}
	}

	$this->SethklRange($hmax, $kmax, $lmax);
	return ($hmax, $kmax, $lmax);
}

1;