#!/usr/bin/perl

use strict;
#use warnings;

use lib 'c:/Programs/Perl/lib';
use lib 'd:/Programs/Perl/lib';
use lib '/home/tkamiya/bin/lib';
use lib '/home/tkamiya/bin';

use MyApplication;
use Crystal::CIF;
use Crystal::Crystal;
use Sci qw($torad);

#===============================================
# Applicationクラス作製
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

# 実行プログラム名（デフォルトでは$0から取得）
my $Program = $App->Program();
# アプリケーション名
$App->SetAppName($Program);
# バージョン
my $Version = $App->SetVersion("Ver 0.1");

#===============================================
# 環境設定
#===============================================
my $ProgramPath = $App->SpeculateProgramPath($0);
$App->print("Program Path: $ProgramPath\n");

#===============================================
# コマンドラインオプション読み込み
#===============================================
$App->AddArgument("--Mode",               "--Mode=[Atom|Site|All]: Default Atom");
$App->AddArgument("--MinimumDistance",    "--MinimumDistance=val: Default 0.1");
$App->AddArgument("--MaximumDistance",    "--MaximumDistance=val: Default 5.0");
$App->AddArgument("--nMesh",              "--nMesh=val: Default 201");
$App->AddArgument("--CoordinationNumber", "--CoordinationNumber=[val]: Default 4");
$App->AddArgument("--help",               "--help    : Show this help");
exit 1 if($App->ReadArgs() != 1);

my $Mode = $App->GetGetArg('Mode');
$Mode = 'Atom' if($Mode =~ /^\s*$/);
my $MinimumDistance = $App->GetGetArg('MinimumDistance');
$MinimumDistance    = 0.1 if($MinimumDistance <= 0.0);
my $MaximumDistance = $App->GetGetArg('MaximumDistance');
$MaximumDistance    = 5.0 if($MaximumDistance <= 0.0);
my $nCN             = $App->GetGetArg("CoordinationNumber");
$nCN                = 4   if($nCN == 0);
my $nMesh           = $App->GetGetArg("nMesh");
$nMesh              = 201 if($nMesh <= 2);
my $CIFFile;
for(my $i = 0 ; ; $i++) {
	my $f = $App->GetGetArg($i);
	last if(!defined $f);
	next if($f =~ /^#/);
	
	$CIFFile = $f;
	last;
}

unless($CIFFile) {
	$App->Usage();
	exit 1;
}

print "Mode               : $Mode\n";
print "Maximum Distance   : $MaximumDistance\n";
print "nMesh              : $nMesh\n";
print "Coordination Number: $nCN\n";
print "CIF File: $CIFFile\n";
print "\n";

#CIFファイル名を引数に与えて読み込む
my $cif = $App->{'CIF'} = new CIF;
unless($cif->Read($CIFFile, 0)) {
	$App->print("Error: Can not read [$CIFFile]\n\n");
	exit 2;
}
$App->print("\n");

my $Crystal = $cif->GetCCrystal();
#$Crystal->SplitPartialSites();
$Crystal->ExpandCoordinates();
#pCrystal->m_ASFDCDBFile  = ASFDCPath;
#pCrystal->m_NonrelDBFile = NONRELPath;
#pCrystal->m_SPGDBFile    = SPGDBPath;
#pCrystal->SetFileName(InputFile);
#(pCrystal->SpaceGroup).SetNoSymmetry();
my $CrystalName = $Crystal->CrystalName();
$App->print("CrystalName: $CrystalName\n");
my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
$App->printf("cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a, $b, $c, $alpha, $beta, $gamma);
#my ($a11, $a12, $a13, $a21, $a22, $a23, $a31, $a32, $a33) = $Crystal->LatticeVectors();
#$App->printf("  Vectors: (%9.6f %9.6f %9.6f)\n", $a11, $a12, $a13);
#$App->printf("           (%9.6f %9.6f %9.6f)\n", $a21, $a22, $a23);
#$App->printf("           (%9.6f %9.6f %9.6f)\n", $a31, $a32, $a33);
#my ($g11, $g12, $g13, $g21, $g22, $g23, $g31, $g32, $g33) = $Crystal->Metrics();
#printf "  Metrics: |%12.6f %12.6f %12.6f|\n", $g11, $g12, $g13;
#printf "           |%12.6f %12.6f %12.6f|\n", $g21, $g22, $g23;
#printf "           |%12.6f %12.6f %12.6f|\n", $g31, $g32, $g33;
#$my $vol = $Crystal->Volume();
#$App->print("  Volume: $vol A^3\n");
#my $Density = $Crystal->Density();
#$App->print("  Density: $Density g/cm^3\n");
#my ($Ra,$Rb,$Rc,$Ralpha,$Rbeta,$Rgamma) = $Crystal->ReciprocalLatticeParameters();
#$App->printf("Reciprocal cell (w/o 2pi): %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $Ra, $Rb, $Rc, $Ralpha, $Rbeta, $Rgamma);
#my ($Ra11, $Ra12, $Ra13, $Ra21, $Ra22, $Ra23, $Ra31, $Ra32, $Ra33) = $Crystal->ReciprocalLatticeVectors();
#	printf "  Vectors: (%9.6f %9.6f %9.6f)\n", $Ra11, $Ra12, $Ra13;
#	printf "           (%9.6f %9.6f %9.6f)\n", $Ra21, $Ra22, $Ra23;
#	printf "           (%9.6f %9.6f %9.6f)\n", $Ra31, $Ra32, $Ra33;
#	my ($Rg11, $Rg12, $Rg13, $Rg21, $Rg22, $Rg23, $Rg31, $Rg32, $Rg33) = $Crystal->RMetrics();
#	printf "  Metrics: |%9.6f %9.6f %9.6f|\n", $Rg11, $Rg12, $Rg13;
#	printf "           |%9.6f %9.6f %9.6f|\n", $Rg21, $Rg22, $Rg23;
#	printf "           |%9.6f %9.6f %9.6f|\n", $Rg31 ,$Rg32, $Rg33;
#	my $Rvol = $Crystal->CalculateReciprocalVolume();
#	$App->print("  Volume: $Rvol A^-3\n");
#my $SPG = $Crystal->GetCSpaceGroup();
#my $SPGName = $SPG->SPGName();
#my $iSPG    = $SPG->iSPG();
#my $LatticeSystem = $SPG->LatticeSystem();
#$App->print("Space Group: $SPGName ($iSPG)\n");
#$App->print("Lattice System: $LatticeSystem\n");

my $rstep = $MaximumDistance / ($nMesh-1);
my $nLatticeX = 1;
my $nLatticeY = 1;
my $nLatticeZ = 1;
my $sinA = abs(sin($alpha * $torad));
my $sinB = abs(sin($beta  * $torad));
my $sinG = abs(sin($gamma * $torad));
$sinA = $sinB if($sinB < $sinA);
$sinA = $sinG if($sinG < $sinA);
my $nx = int($MaximumDistance / $a / $sinA);
my $ny = int($MaximumDistance / $b / $sinA);
my $nz = int($MaximumDistance / $c / $sinA);
$nLatticeX += $nx;
$nLatticeY += $ny;
$nLatticeZ += $nz;
printf("  Lattice Range: %d %d %d\n", $nLatticeX, $nLatticeY, $nLatticeZ);
print "\n";

my @AtomTypeList     = $Crystal->GetCAtomTypeList();
my $nAtomType        = @AtomTypeList;
my @AtomAsymSiteList = $Crystal->GetCAsymmetricAtomSiteList();
my $nAtomAsymSites   = @AtomAsymSiteList;
my @AtomSiteList     = $Crystal->GetCExpandedAtomSiteList();
my $nAtomSites       = $Crystal->nTotalExpandedAtomSite();

#$Mode = 'Atom';
#$Mode = 'Site';
#$Mode = 'All';
my (@site, @id, @mult);
if($Mode =~ /atom/i) {
	@site = @AtomTypeList;
	for(my $i = 0; $i < $nAtomSites ; $i++) {
		my $atom = $AtomSiteList[$i];
		my $name = $atom->AtomNameOnly();
		$id[$i] = $Crystal->FindiAtomType($name) - 1;
		$mult[$id[$i]]++;
#print "i=$i: id=$id[$i] $name m=$mult[$id[$i]]\n";
	}
}
elsif($Mode =~ /site/i) {
	@site = @AtomAsymSiteList;
	for(my $i = 0; $i < $nAtomSites ; $i++) {
		my $atom = $AtomSiteList[$i];
		my $name = $atom->AtomNameOnly();
		my $mult = $atom->Multiplicity();
		$id[$i] = $atom->Id() - 1;
		$mult[$id[$i]] = $mult;
	}
}
elsif($Mode =~ /all/i) {
	@site = @AtomSiteList;
	for(my $i = 0; $i < $nAtomSites ; $i++) {
		$id[$i] = $i;
		$mult[$id[$i]]++;
	}
}
else {
	print "\n";
	print "Error: Invalide Mode [$Mode], must be Atom, Site, or All\n";
	exit;
}

print "\n";
$App->print("nAtomType: $nAtomType\n");
for(my $i = 0 ; $i < $nAtomType ; $i++) {
	my $atom0  = $AtomTypeList[$i];
	my $name0          = $atom0->AtomNameOnly();
	if($Mode =~ /atom/i) {
		printf "  %3d: %2s  m=%2d\n", $i, $name0, $mult[$i];
	}
	else {
		printf "  %3d: %s\n", $i, $name0;
	}
}
print "\n";
$App->print("nAsymmetricAtomSite: $nAtomAsymSites\n");
for(my $i = 0 ; $i < $nAtomAsymSites ; $i++) {
	my $atom0  = $AtomAsymSiteList[$i];
	my $name0          = $atom0->AtomNameOnly();
	my ($x0, $y0, $z0) = $atom0->Position(1);
	my $occ0           = $atom0->Occupancy();
	my $mult0          = $atom0->Multiplicity();
	if($Mode =~ /site/i) {
		printf "  %3d: %2s (%8.4f, %8.4f, %8.4f)  occ=%6.3f  mult=%2d  m=%2d\n",
			$i, $name0, $x0, $y0, $z0, $occ0, $mult0, $mult[$i];
	}
	else {
		printf "  %3d: %2s (%8.4f, %8.4f, %8.4f)  occ=%6.3f  mult=%2d\n",
			$i, $name0, $x0, $y0, $z0, $occ0, $mult0;
	}
}
print "\n";
$App->print("nTotalSite: $nAtomSites\n");
for(my $i = 0 ; $i < $nAtomSites ; $i++) {
	my $atom0  = $AtomSiteList[$i];
	my $name0          = $atom0->AtomNameOnly();
	my ($x0, $y0, $z0) = $atom0->Position(1);
	my $occ0           = $atom0->Occupancy();
	my $mult0          = $atom0->Multiplicity();
	my $id0            = $atom0->Id() - 1;
	if($Mode =~ /all/i) {
		printf "  %3d: %2s (%8.4f, %8.4f, %8.4f)  occ=%6.3f  mult=%2d  SiteID=%3d  ModeID=%2d  \m=%2d\n",
			$i, $name0, $x0, $y0, $z0, $occ0, $mult0, $id0, $id[$i], $mult[$i];
	}
	else {
		printf "  %3d: %2s (%8.4f, %8.4f, %8.4f)  occ=%6.3f  mult=%2d  SiteID=%3d  ModeID=%2d\n",
			$i, $name0, $x0, $y0, $z0, $occ0, $mult0, $id0, $id[$i];
	}
}

&CalculateRDF($Crystal, \@site, \@AtomTypeList, \@AtomSiteList, \@id, \@mult);
exit;

my @nCN = Utils::Split("\\s*,\\s*", $nCN);
for(my $i = 0 ; $i < @nCN ; $i++) {
	next if($nCN[$i] <= 0);

	&CalculateCNRDF($Crystal, $nCN[$i], \@site, \@AtomSiteList);
}

exit 1;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# サブルーチン
#==========================================
sub CalculateRDF
{
	my ($Crystal, $pSite, $pAtomType, $pAllSite, $pId, $pMult) = @_;

	my $nAtomType = @$pAtomType;
	my $nSite     = @$pSite;
	my $nAllSite  = @$pAllSite;
	
	printf("\n");
	printf("Calculate RDF.\n");

	my (@r, @RDF, @RDFtot, @RCN, @RCNtot);
	my %iAtomType;
	for(my $k = 0 ; $k < $nMesh ; $k++) {
		$r[$k] = $rstep * $k;
		$RDFtot[$k] = 0.0;
	}
	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $atom = $pSite->[$i];
		my $name = $atom->AtomNameOnly();
		$iAtomType{$name} = $Crystal->FindiAtomType($name) - 1;
#print "$i: $name: $iAtomType{$name}\n";
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF[$i][$j][$k] = 0.0;
				$RCN[$i][$j][$k] = 0.0;
			}
		}
	}

	my @ndiv;
	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		my $id0 = $id[$i];
		$ndiv[$id0]++;
	}

	print "\n";
	print "ndiv:\n";
	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $atom0 = $pSite->[$i];
		my $name0 = $atom0->AtomNameOnly();
		print "  Site $name0$i: ndiv=$ndiv[$i]\n";
	}

	for(my $ia = 0 ; $ia < $nAllSite ; $ia++) {
		my $atom0  = $pAllSite->[$ia];
		my $name0          = $atom0->AtomNameOnly();
		my ($x0, $y0, $z0) = $atom0->Position(1);
		my $occ0           = $atom0->Occupancy();
		my $iAtomType0     = $iAtomType{$name0};
		my $id0            = $id[$ia];
#printf("ia=%d name0=%s  iAtomType0=%d: Count=%d\n", ia, name0, iAtomType0, nCountedAtoms[iAtomType0]);

		for(my $iz = -$nLatticeZ ; $iz <= $nLatticeZ ; $iz++) {
		for(my $iy = -$nLatticeY ; $iy <= $nLatticeY ; $iy++) {
		for(my $ix = -$nLatticeX ; $ix <= $nLatticeX ; $ix++) {
			for(my $i = 0 ; $i < $nAllSite ; $i++) {
				my $atom1  = $pAllSite->[$i];
				my $name1          = $atom1->AtomNameOnly();
				my ($x1, $y1, $z1) = $atom1->Position(1);
				my $occ1           = $atom1->Occupancy();
				my $iAtomType1     = $iAtomType{$name1};
				my $id1            = $id[$i];

				my $dis = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
				next if($dis < $MinimumDistance);
				next if($dis > $MaximumDistance);

				my $idx = int($dis / $rstep);
				$RDF[$id0][$iAtomType1][$idx] += $occ1;
				$RDFtot[$idx] += $occ1;
			}
		}
		}
		}
	}
	for(my $i = 0 ; $i < $nSite ; $i++) {
#print "ndiv=$ndiv[$i]\n";
		for(my $j = 0 ; $j < $nSite ; $j++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF[$i][$j][$k] /= $ndiv[$i];
			}
		}
	}

	for(my $i = 0 ; $i < $nSite ; $i++) {
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			for(my $k = 1 ; $k < $nMesh ; $k++) {
				$RCN[$i][$j][$k] = $RCN[$i][$j][$k-1] + $RDF[$i][$j][$k];
				$RDF[$i][$j][$k] /= $rstep;
#print "$k: $r[$k]\t$RDF[$i][$j][$k]\t$RCN[$i][$j][$k]\n";
				$RDFtot[$k] += $RDF[$i][$j][$k];
				$RCNtot[$k] += $RCN[$i][$j][$k];
			}
		}
	}

	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $atom0 = $pSite->[$i];
		my $name0 = $atom0->AtomNameOnly();
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			my $atom1 = $pAtomType->[$j];
			my $name1 = $atom1->AtomNameOnly();
#print "$i,$j: $name0 $name1\n";

			my $fname = "$name0$i-$name1-RDF.csv";
			&SaveCSV($fname, "r(A),RDF($name0-$name1),RCN($name0-$name1)", \@r,$RDF[$i][$j], $RCN[$i][$j]);
		}
	}
#	&SaveCSV("Total-RDF.csv", "r(A),RDF(total),RCN(total)", \@r,\@RDFtot, \@RCNtot);
#	print "\n";

	return 1;
}

sub CalculateCNRDF
{
	my ($Crystal, $nCN, $pSite, $pAllSite) = @_;

	my $nSite = @$pSite;
	my $nAllSite = @$pAllSite;

	printf("\n");
	printf("Calculate CN RDF for coordination number $nCN.\n");

	my (@CNRDF, @RDF);
	my %iAtomType;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom = $AtomTypeList[$i];
		my $name = $atom->AtomNameOnly();
#print "i=$i: $name\n";
		$iAtomType{$name} = $i;
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF[$i][$j][$k] = 0.0;
				$CNRDF[$i][$j][$k] = 0.0;
			}
		}
	}

	for(my $ia = 0 ; $ia < $nAtomAsymSites ; $ia++) {
		my $atom0  = $AtomSiteList[$ia];
		my $name0          = $atom0->AtomNameOnly();
		my ($x0, $y0, $z0) = $atom0->Position(1);
		my $occ0           = $atom0->Occupancy();
		my $iAtomType0     = $iAtomType{$name0};
#printf("ia=%d name0=%s  iAtomType0=%d: Count=%d\n", ia, name0, iAtomType0, nCountedAtoms[iAtomType0]);

		for(my $iz = -$nLatticeZ ; $iz <= $nLatticeZ ; $iz++) {
		for(my $iy = -$nLatticeY ; $iy <= $nLatticeY ; $iy++) {
		for(my $ix = -$nLatticeX ; $ix <= $nLatticeX ; $ix++) {
			for(my $i = 0 ; $i < $nAtomSites ; $i++) {
				my $atom1  = $AtomSiteList[$i];
				my $name1          = $atom1->AtomNameOnly();
				my ($x1, $y1, $z1) = $atom1->Position(1);
				my $occ1           = $atom1->Occupancy();
				my $iAtomType1     = $iAtomType{$name1};

				my $dis = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
				next if($dis > $MaximumDistance);

				my $idx = int($dis / $rstep);
				$CNRDF[$iAtomType0][$iAtomType1][$idx] += $occ1;
			}
		}
		}
		}

		for(my $iAtomType1 = 0 ; $iAtomType1 < $nAtomType ; $iAtomType1++) {
#距離 r = rstep*k での全配位数
			my $ncn = 0;
			my $bClear = 0;
			for(my $k = 1 ; $k < $nMesh ; $k++) {
				$ncn += $CNRDF[$iAtomType1][$k];
				if(!$bClear and $ncn == $nCN) {
					$RDF[$iAtomType0][$iAtomType1][$k] += 1;
					$bClear = 1;
				}
				if($ncn >= $nCN+1) {
					if($bClear) {
						$RDF[$iAtomType0][$iAtomType1][$k] -= 1;
					}
					last;
				}
			}
		}
	}

	my (@r, @RCN);
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		for(my $j = $i ; $j < $nAtomType ; $j++) {
			$r[0] = 0.0;
			$RCN[$i][$j][0] = 0.0;
			for(my $k = 1 ; $k < $nMesh ; $k++) {
				$r[$k] = $rstep * $k;
				$RCN[$i][$j][$k] = $RCN[$i][$j][$k-1] + $RDF[$i][$j][$k];
				$RDF[$i][$j][$k] /= $rstep;
#print "$k: $r[$k]\t$RDF[$i][$j][$k]\t$RCN[$i][$j][$k]\n";
			}
		}
	}

	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom0  = $AtomTypeList[$i];
		my $name0  = $atom0->AtomNameOnly();

		for(my $j = $i ; $j < $nAtomType ; $j++) {
			my $atom1  = $AtomTypeList[$j];
			my $name1  = $atom1->AtomNameOnly();

			my $fname = "$name0-$name1-RDF-CN$nCN.csv";
			&SaveCSV($fname, "r(A),RDF($name0-$name1),RCN($name0-$name1)", \@r,$RDF[$i][$j], $RCN[$i][$j]);
		}
	}

	return 1;
}

sub SaveCSV
{
	my ($fname, $header, $px, $py1, $py2) = @_;

	printf("Save to [%s]...\n", $fname);
	my $out = JFile->new($fname, 'w') or die "Can not write to [$fname].\n";
	$out->printf("$header\n");
	for(my $k = 0 ; $k < $nMesh ; $k++) {
		$out->printf("%g,%g,%g\n", $px->[$k], $py1->[$k], $py2->[$k]);
	}
	$out->Close();
}