#!/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クラス作製
#===============================================
sub terminate
{
	my $ret = 0;
	if(scalar @_ > 0) {
		$ret = $_[0];
	}
	print("\nPress ENTER to terminate>> ");
	<STDIN>;
	exi();
}

my $App = new MyApplication;
terminate() 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 null");
$App->AddArgument("--help",               "--help    : Show this help");
terminate(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                = '' 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();
	terminate(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");
	terminate(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, $nLatticeY, $nLatticeZ) = $Crystal->GetLatticeRange($MaximumDistance);
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";
	terminate();
}

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];
	}
}

my @nCN = Utils::Split("\\s*,\\s*", $nCN);
if(@nCN == 0) {
	&CalculateRDFs($Crystal, $rstep, $nMesh, \@site, \@AtomTypeList, \@AtomSiteList, \@id, \@mult);
}
else {
	&CalculateCNRDFs($Crystal, \@nCN, $rstep, $nMesh, \@site, \@AtomTypeList, \@AtomSiteList, \@id, \@mult);
}

terminate();

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

#==========================================
# サブルーチン
#==========================================
# なぜか、引数のポインタで与えられた$pRDFがクリアされるので使えない。CalRDFiを使う
sub AddRDFi
{
	my ($Crystal, $rstep, $nMesh, $MinimumDistance, $MaximumDistance, $pRDF, 
			$x0, $y0, $z0, $id0, $pAllSite, $pId, $pMult) = @_;

	my $nAllSite  = @$pAllSite;

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

	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		my $atom1          = $pAllSite->[$i];
		my $name1          = $atom1->AtomNameOnly();
		my $iAtomType1     = $Crystal->FindiAtomType($name1) - 1;
		my ($x1, $y1, $z1) = $atom1->Position(1);
		my $occ1           = $atom1->Occupancy();
#		my $id1            = $pId->[$i];
#print "id=$id0 - $name1 ia=$iAtomType1\n";
		for(my $iz = -$nLatticeZ ; $iz <= $nLatticeZ ; $iz++) {
		for(my $iy = -$nLatticeY ; $iy <= $nLatticeY ; $iy++) {
		for(my $ix = -$nLatticeX ; $ix <= $nLatticeX ; $ix++) {
			my $dis = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
#print "  d=$dis (r)=[$MinimumDistance, $MaximumDistance] (i)=($ix,$iy,$iz)\n";
			next if($dis < $MinimumDistance);
			next if($dis > $MaximumDistance);

			my $idx = int($dis / $rstep);
			$pRDF->[$iAtomType1][$idx] += $occ1 / $ndiv[$id0] / $rstep;
#print "  d=$dis idx=$idx: RDF=$pRDF->[$iAtomType1][$idx] (occ=$occ1)\n";
		}
		}
		}
	}

	return $pRDF;
}

sub CalRDFi
{
	my ($Crystal, $rstep, $nMesh, $MinimumDistance, $MaximumDistance, 
			$x0, $y0, $z0, $id0, $pAllSite, $pId, $pMult) = @_;

	my $nAllSite  = @$pAllSite;
	my $pRDF;
	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		my $atom1          = $pAllSite->[$i];
		my $name1          = $atom1->AtomNameOnly();
		my $iAtomType1     = $Crystal->FindiAtomType($name1) - 1;
		for(my $k = 0 ; $k < $nMesh ; $k++) {
			$pRDF->[$iAtomType1][$k] = 0.0;
		}
	}

	my @ndiv;
	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		my $idt = $pId->[$i];
		$ndiv[$idt]++;
	}

	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		my $atom1          = $pAllSite->[$i];
		my $name1          = $atom1->AtomNameOnly();
		my $iAtomType1     = $Crystal->FindiAtomType($name1) - 1;
		my ($x1, $y1, $z1) = $atom1->Position(1);
		my $occ1           = $atom1->Occupancy();
#		my $id1            = $pId->[$i];
#print "id=$id0 - $name1 ia=$iAtomType1\n";
		for(my $iz = -$nLatticeZ ; $iz <= $nLatticeZ ; $iz++) {
		for(my $iy = -$nLatticeY ; $iy <= $nLatticeY ; $iy++) {
		for(my $ix = -$nLatticeX ; $ix <= $nLatticeX ; $ix++) {
			my $dis = $Crystal->GetInterAtomicDistance($x0, $y0, $z0, $x1+$ix, $y1+$iy, $z1+$iz);
#print "  d=$dis (r)=[$MinimumDistance, $MaximumDistance] (i)=($ix,$iy,$iz)\n";
			next if($dis < $MinimumDistance);
			next if($dis > $MaximumDistance);

			my $idx = int($dis / $rstep);
			$pRDF->[$iAtomType1][$idx] += $occ1 / $ndiv[$id0] / $rstep;
#print "  d=$dis idx=$idx: RDF=$pRDF->[$iAtomType1][$idx] (occ=$occ1)\n";
		}
		}
		}
	}

	return $pRDF;
}

sub CalculateRDFs
{
	my ($Crystal, $rstep, $nMesh, $pSite, $pAtomType, $pAllSite, $pId, $pMult) = @_;

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

	my @RDF;
	for(my $ia = 0 ; $ia < $nSite ; $ia++) {
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF[$ia][$i][$k] = 0.0;
			}
		}
	}

	my ($pr, @RDF, @RCN);
	for(my $k = 0 ; $k < $nMesh ; $k++) {
		$pr->[$k] = $rstep * $k;
	}

	for(my $ia = 0 ; $ia < $nAllSite ; $ia++) {
		my $atom0       = $pAllSite->[$ia];
		my $name0       = $atom0->AtomNameOnly();
		my $iAtomType0  = $Crystal->FindiAtomType($name0) - 1;
		my ($x0, $y0, $z0) = $atom0->Position(1);
		my $id0         = $pId->[$ia];
# なぜか、引数のポインタで与えられた$pRDFがクリアされるので使えない。CalRDFiを使う
#		my $pRDF = &AddRDFi($Crystal, $rstep, $nMesh, $MinimumDistance, $MaximumDistance, $RDF[$id0], 
		my $pRDF = &CalRDFi($Crystal, $rstep, $nMesh, $MinimumDistance, $MaximumDistance,, 
				$x0, $y0, $z0, $id0, $pAllSite, $pId, $pMult);
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF[$id0][$i][$k] += $pRDF->[$i][$k];
			}
		}
	}

	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $atom = $pSite->[$i];
		my $name = $atom->AtomNameOnly();
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RCN[$i][$j][$k] = 0.0;
			}
		}
	}

	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] * $rstep;
#				$RCNtot[$k] += $RCN[$i][$j][$k];
			}
		}
	}
#&SaveCSV("c.csv", "r(A),RDF(Ga-O),RCN(Ga-O)", $pr,$RDF[2][3], $RCN[2][3]);
	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)", $pr,$RDF[$i][$j], $RCN[$i][$j]);
		}
	}
#	&SaveCSV("Total-RDF.csv", "r(A),RDF(total),RCN(total)", $pr, $pRDFtot, \@RCNtot);
#	print "\n";

	return 1;
}

sub CalculateCNRDFs
{
	my ($Crystal, $pnCN, $rstep, $nMesh, $pSite, $pAtomType, $pAllSite, $pId, $pMult) = @_;

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

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

	printf("\n");
	printf("Calculate CN RDF.\n");

	my $pr;
	my (@RDF0, @RCN0);
	my (@RDF, @CNRDF);
	my (@RCN, @CNRCN);
	for(my $ia = 0 ; $ia < $nAllSite ; $ia++) {
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF0[$ia][$i][$k] = 0.0;
			}
		}
	}
	for(my $ia = 0 ; $ia < $nSite ; $ia++) {
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF[$ia][$i][$k] = 0.0;
				for(my $icn = 0 ; $icn < @$pnCN ; $icn++) {
					my $n = $pnCN->[$icn];
					$CNRDF[$n][$ia][$i][$k] = 0.0;
				}
			}
		}
	}
	for(my $k = 0 ; $k < $nMesh ; $k++) {
		$pr->[$k] = $rstep * $k;
	}
	my @ndiv;
	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		my $id0 = $id[$i];
		$ndiv[$id0]++;
	}

	my (@Id0, @Mult0);
	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		$Id0[$i] = $i;
		$Mult0[$i] = 1;
	}
	for(my $ia = 0 ; $ia < $nAllSite ; $ia++) {
		my $atom0       = $pAllSite->[$ia];
		my $name0       = $atom0->AtomNameOnly();
		my $iAtomType0  = $Crystal->FindiAtomType($name0) - 1;
		my ($x0, $y0, $z0) = $atom0->Position(1);
		my $id0         = $pId->[$ia];

# なぜか、引数のポインタで与えられた$pRDFがクリアされるので使えない。CalRDFiを使う
#		my $pRDF = &AddRDFi($Crystal, $rstep, $nMesh, $MinimumDistance, $MaximumDistance, $RDF0[$id0], 
		my $pRDF = &CalRDFi($Crystal, $rstep, $nMesh, $MinimumDistance, $MaximumDistance,, 
				$x0, $y0, $z0, $id0, $pAllSite, \@Id0, \@Mult0);
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			for(my $k = 0 ; $k < $nMesh ; $k++) {
				$RDF0[$ia][$i][$k] += $pRDF->[$i][$k];
				$RDF[$id0][$i][$k] += $pRDF->[$i][$k];
			}
		}
	}

	for(my $i = 0 ; $i < $nAllSite ; $i++) {
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			for(my $k = 1 ; $k < $nMesh ; $k++) {
				$RCN0[$i][$j][$k] = $RCN0[$i][$j][$k-1] + $RDF0[$i][$j][$k] * $rstep;
			}
		}
	}

	for(my $ia = 0 ; $ia < $nAllSite ; $ia++) {
		my $id0 = $pId->[$ia];
		for(my $i = 0 ; $i < $nAtomType ; $i++) {
			for(my $icn = 0 ; $icn < @$pnCN ; $icn++) {
				my $tnCN = $pnCN->[$icn];
				for(my $k = 0 ; $k < $nMesh ; $k++) {
					if($RCN0[$ia][$i][$k] > $tnCN - 1.0e-4) {
#my $r = $rstep * $k;
#print "hit: ia=$id0 (id0=$id0) iAtom=$i k=$k (r=$r) for nCN=$tnCN  RCN0=$RCN0[$ia][$i][$k]\n";
						$CNRDF[$tnCN][$id0][$i][$k]++;
						last;
#last if($tnCN >= 7);
					}
				}
			}
		}
	}

if(0) {
for(my $i = 0 ; $i < 8 ; $i++) {
&SaveCSV("Ga$i-O.csv", "r(A),RDF(Ga$i-O),RCN(Ga$i-O)", $pr,$RDF0[$i][1], $RCN0[$i][1]);
}
for(my $i = 0 ; $i < $nSite ; $i++) {
&SaveCSV("Ga$i-O-CNRDF4.csv", "r(A),RDF(Ga$i-O),RCN(Ga$i-O)", $pr,$CNRDF[4][$i][1], $RCN0[$i][1]);
}
#terminate();
}

	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] * $rstep;
				for(my $icn = 0 ; $icn < @$pnCN ; $icn++) {
					my $n = $pnCN->[$icn];
					$CNRCN[$n][$i][$j][$k] 
						= $CNRCN[$n][$i][$j][$k-1] + $CNRDF[$n][$i][$j][$k];
				}
			}
		}
	}

	for(my $i = 0 ; $i < $nSite ; $i++) {
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			for(my $k = 1 ; $k < $nMesh ; $k++) {
				$RDF[$i][$j][$k] /= $ndiv[$i];
				$RCN[$i][$j][$k] /= $ndiv[$i];
			}
		}
	}

	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 $fname1 = "$name0$i-$name1-RDF.csv";
			&SaveCSV($fname1, "r(A),RDF($name0-$name1),RCN($name0-$name1)", 
				$pr,$RDF[$i][$j], $RCN[$i][$j]);
			for(my $icn = 0 ; $icn < @$pnCN ; $icn++) {
				my $n = $pnCN->[$icn];
				next if(!defined $CNRCN[$n][$i]);
				my $fname2 = "$name0$i-$name1-CNRDF$n.csv";
				&SaveCSV($fname2, "r(A),CNRDF$n($name0-$name1),CNRCN$n($name0-$name1)", 
					$pr,$CNRDF[$n][$i][$j], $CNRCN[$n][$i][$j]);
			}
		}
	}
#	&SaveCSV("Total-RDF.csv", "r(A),RDF(total),RCN(total)", $pr, $pRDFtot, \@RCNtot);
#	print "\n";

	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();
}