#!/usr/bin/perl

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

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

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

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

my $CrystalName = $Crystal->CrystalName();
print "CrystalName: $CrystalName\n";
my $SPG = $Crystal->GetCSpaceGroup();
my $LatticeSystem = $SPG->LatticeSystem();
print "LatticeSystem: $LatticeSystem\n";
my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "cell: $a $b $c  $alpha $beta $gamma\n";

my @AtomTypeList = $Crystal->GetCAtomTypeList();
my $nAtomType = @AtomTypeList;
print "nAtomType: $nAtomType\n";
for(my $i = 0 ; $i < $nAtomType ; $i++) {
	my $label    = $AtomTypeList[$i]->Label();
	my $atomname = $AtomTypeList[$i]->AtomNameOnly();
	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";
}

$Crystal->ExpandCoordinates();
my @ExpandedAtomSiteList   = $Crystal->GetCExpandedAtomSiteList();
my $nTotalExpandedAtomSite = $Crystal->nTotalExpandedAtomSite();
my @nMultiplicityExpandedAtomSiteList = $Crystal->GetCnMultiplicityExpandedAtomSiteList();
print "nTotalExpandedAtomSite: $nTotalExpandedAtomSite\n";
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]$label ($type): ($x, $y, $z) [$occupancy][$mult]\n";
}

my $Rietan = new Rietan;
my $out = new JFile($OutFile, "w");
$out->print("s,asf(X-ray),asf(Electron),asf(H)\n");
for(my $s = 0.1 ; $s < 25.0 ; $s += 0.1) {
	my $asfX  = $Crystal->asf(0, $s);
	my $asfE  = $Crystal->asfElectron(0, $s);
	my $asfH  = $Rietan->HydrogenAtomicScatteringFactor($s);
#	print "$s,$asfX,$asfE\n";
	$out->print("$s,$asfX,$asfE,$asfH\n");
}
$out->Close();

exit;