#!/usr/bin/perl

use lib 'D:/tkProg/tkProg.main/tklib/Perl/lib';

use strict;
use Sci qw($pi $todeg $torad asin);
use Crystal::CIF;
use Crystal::Crystal;
use Crystal::RietanFP;


my $XRDFile  = "XRD.csv";
my $CIFFile  = "d:/MyWebs/Research/CrystalStructure/Si5.cif"; #"Si5.cif";
   $CIFFile  = 'BaCeFe2As2.cif'; #"d:/LaCuOSe.cif";

my $lambda   = 0.154; # nm
my $thickness = 0; #150; # nm, set to zero for bulk

print("\n");
print("==============================================\n");
print("CIFFile: $CIFFile\n");
print("lambda : $lambda nm\n");
print("==============================================\n");
print("\n");

my $w = 0.1;
my $Q2Start =   5.0;
my $Q2Stop  = 150.0;
my $Q2Step  = 0.02;

my $Biso = 0.0;


$CIFFile = $ARGV[0] if(@ARGV >= 1);
$Q2Start = $ARGV[1] + 0.0 if(@ARGV >= 2);
$Q2Stop  = $ARGV[2] + 0.0 if(@ARGV >= 3);
$Q2Step  = $ARGV[3] + 0.0 if(@ARGV >= 4);
$w       = $ARGV[4] + 0.0 if(@ARGV >= 5);


print("\n");
print("========= simulate powder XRD =========");
print("CIF: $CIFFile\n");
print("X-ray wavelength: $lambda nm\n");
#print("Sample thickness: $thickness nm\n");
print("2Theta range: $Q2Start - $Q2Stop, $Q2Step step\n");

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

print "\n";
print "Crystal Information";
my $CrystalName = $Crystal->CrystalName();
print "CrystalName: $CrystalName\n";
my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "cell: $a $b $c  $alpha $beta $gamma\n";
my ($g11,$g12,$g13,$g21,$g22,$g23,$g31,$g32,$g33) = $Crystal->Metrics();
print "Metrics: $g11 $g22 $g33  $g12 $g23 $g31\n";
my ($Rg11,$Rg12,$Rg13,$Rg21,$Rg22,$Rg23,$Rg31,$Rg32,$Rg33) = $Crystal->RMetrics();
print "Metrics in reciprocal space: $Rg11 $Rg22 $Rg33  $Rg12 $Rg23 $Rg31\n";
my $vol = $Crystal->Volume();
print "Volume: $vol A^3\n";
my $Density = $Crystal->Density();
print "Density: $Density g/cm3\n";

my $gmax = 2.0 * sin($Q2Stop / 2.0 * 3.1416/180.0);
my $nx1 = int($gmax / $Rg11 + 2);
my $ny1 = int($gmax / $Rg22 + 2);
my $nz1 = int($gmax / $Rg33 + 2);
my $nx0 = -$nx1;
my $ny0 = -$ny1;
my $nz0 = -$nz1;

print("\n");
my $SPG = $Crystal->GetCSpaceGroup();
my $SPGName = $SPG->SPGName();
my $iSPG    = $SPG->iSPG();
my $LatticeSystem = $SPG->LatticeSystem();
print "Space Group: $SPGName ($iSPG)\n";
print "Lattice System: $LatticeSystem\n";
print "\n";

my @AtomTypeList                      = $Crystal->GetCAtomTypeList();
my $nAtomType                         = @AtomTypeList;
my @ExpandedAtomSiteList              = $Crystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite            = $Crystal->nTotalExpandedAtomSite();
my @nMultiplicityExpandedAtomSiteList = $Crystal->GetCnMultiplicityExpandedAtomSiteList();

print "nAtomType: $nAtomType\n";
for(my $i = 0 ; $i < $nAtomType ; $i++) {
	my $label    = $AtomTypeList[$i]->Label();
	my $atomname = $AtomTypeList[$i]->AtomName();
	my $charge   = $AtomTypeList[$i]->Charge();
	my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
	my $AtomicMass   = $AtomTypeList[$i]->AtomicMass();
	my $i1 = $i+1;
	print "   #$i1: $label : $atomname [$AtomicNumber] ($charge) ($AtomicMass)\n";
}
print "\n";

print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";
#my $mu = 0.0;
for(my $i = 0 ; $i < $nTotalExpandedAtomSite ; $i++) {
	my $label     = $ExpandedAtomSiteList[$i]->Label();
	my $type      = $ExpandedAtomSiteList[$i]->AtomName();
	my ($x,$y,$z) = $ExpandedAtomSiteList[$i]->Position(1);
	my $occupancy = $ExpandedAtomSiteList[$i]->Occupancy();
	my $id = $ExpandedAtomSiteList[$i]->IdAsymmetricAtomSite();
	my $mult = $nMultiplicityExpandedAtomSiteList[$id-1];
	my $i1 = $i+1;
	print "   $i1: [id=$id]$label ($type): ($x, $y, $z) [occ=$occupancy][mult=$mult]\n";

#	if($MuRho{$type} eq '') {
#		print "Error: No MuRho data is provided for [$type].\n";
#		exit;
#	}
#	$mu += $MuRho{$type} * $Density / $nTotalExpandedAtomSite;
}
print "\n";

#print "mu = $mu cm-1\n";
#print "\n";


print "Thermal vibration parameters\n";
print("Biso=$Biso\n");

my @XRD;
my $nQ2     = int( ($Q2Stop - $Q2Start) / $Q2Step + 1.001);
my $nQ2w    = int( 20.0 * $w / $Q2Step);
for(my $i = 0 ; $i < $nQ2 ; $i++) {
	$XRD[$i] = 0.0;
}

for(my $h = $nx0 ; $h <= $nx1 ; $h++) {
for(my $k = $ny0 ; $k <= $ny1 ; $k++) {
for(my $l = $nz0 ; $l <= $nz1 ; $l++) {
	next if($h == 0 and $k == 0 and $l == 0);

	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);
#print " m=$m\n";
	my ($Fr, $Fi) = $Crystal->Fhkl($h, $k, $l);
	my $F         = sqrt($Fr*$Fr + $Fi*$Fi);
	next if($F < 1.0e-5);

	my $s = $sinQB / $lambda;
	my $P = CrystalObject::PolarizationFactor($QB);
	my $L = CrystalObject::LorentzFactor($QB);
	my $T = CrystalObject::TemperatureFactor($Biso, $s);
#	my $Abs = ($thickness == 0)? 1.0 : 1.0 - exp(-2.0 * $mu * $thickness*1.0e-7 / sin($QB));
	my $I = $F * $F * $P * $L * $T * $T;
#	my $I = $F * $F * $P * $L * $T * $T * $Abs;
#	my $I = $F * $F * $m * $P * $L * $T * $T; #等価な面も再計算しているので、多重度は考慮してはいけない
	my $dhkl = $Crystal->CalculateHKLInterplanarSpacing($h, $k, $l);
	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;

	my $iQ0 = int( ($QB2 - $Q2Start) / $Q2Step + 0.5);
	printf "   2QB=%7.3f P=%7.3f L=%7.3f T=%7.3f   I=%10.3f\n", 2.0*$QB, $P, $L, $T, $I;

	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);
		$XRD[$j] += $I * $Lorentz;
	}
}
}
}

print("\n");
print("Save CSV file to [$XRDFile]\n");
my $out = new JFile($XRDFile, "w");
$out->print("2Theta,I,P,L,T\n");
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");
}
$out->Close();

#my $r  = $Crystal->{ASFDCRed0};
#my $p  = $r->{pASFDCParameters};
#print "p=", join(', ', @$p), "\n";

print("\n");
print("Press ENTER to exit>> ");
<STDIN>;
exit;