#!/usr/bin/perl

use lib 'd:/Programs/Perl/lib';
use lib 'c:/Programs/Perl/lib';
#use lib 'd:/MyWebs/cgi-bin/lib';

use lib '/home/tkamiya/bin/lib';
use lib '/home/tkamiya/bin';

use MyApplication;
use Crystal::CIF;
use Crystal::Crystal;

use strict;
#use warnings;

#===============================================
# Applicationクラス作製
#===============================================
sub terminate
{
	my $ret = 0;
	if(scalar @_ > 0) {
		$ret = $_[0];
	}
	print("\nPress ENTER to terminate>> ");
	<STDIN>;
	exit();
}

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("--MaximumDistance", "--MaximumDistance=val: Default 3.5");
$App->AddArgument("--MaximumNumber",   "--MaximumNumber=val  : Default 100");
$App->AddArgument("--DoSort",          "--DoSort=[0|1]       : Default 0");
$App->AddArgument("--help",            "--help    : Show this help");
terminate(1) if($App->ReadArgs() != 1);

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

&CalculateDistanceAngle();

terminate(0);

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

#==========================================
# サブルーチン
#==========================================

sub Sort
{
	my ($p) = @_;
	my @std = sort { $a->{'dis'} <=> $b->{'dis'} } @$p;
	return @std;
}

sub CalculateDistanceAngle
{
	my $MaximumDistance = $App->GetGetArg('MaximumDistance');
	$MaximumDistance    = 3.5 if($MaximumDistance <= 0.0);
	my $MaximumNumber   = $App->GetGetArg('MaximumNumber');
	$MaximumNumber      = 100 if($MaximumNumber <= 0);
	my $DoSort          = $App->GetGetArg('DoSort');
	$DoSort             = 0 if($DoSort eq '');

	$App->print("\n\n");
	$App->print("Calculate inter-atomic distances and angles\n");
	$App->print("CIF File: $CIFFile\n");
	$App->print("Maximum distance: $MaximumDistance A\n");
	$App->print("Maximum number to be shown: $MaximumNumber\n");
	$App->print("Sort by distance and ion names: $DoSort\n");

#CIFクラスの内容から、Crystalクラスを作成
#作成した時点で、Metricsの計算、イオン情報(原子量など)の読み込み、
#格子体積の計算は終わっている
my $Crystal = $cif->GetCCrystal();
$Crystal->SplitPartialSites();

	$App->print("\n");
	my $CrystalName = $Crystal->CrystalName();
	$App->print("CrystalName: $CrystalName\n");
	my ($a1,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
	$App->printf("cell: %9.6f %9.6f %9.6f   %9.6f %9.6f %9.6f \n", $a1, $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");
#対称要素を使い、非対称単位中の原子を展開する
#この際、Densityも計算される
	$Crystal->ExpandCoordinates();
	my $Density = $Crystal->Density();
	$App->print("  Density: $Density g/cm^3\n");

	$App->print("\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");

	$App->print("\n");
#Crystalクラスから、SpaceGroupクラスを取り出す
#この時点で、並進対称要素の抽出、LatticeSystemの抽出は完了している
	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");

	$App->print("\n");
	my $nTranslation = $SPG->nTranslation();
#	print "nTranslation: $nTranslation\n";
	for(my $i = 0 ; $i < $nTranslation ; $i++) {
		my ($x,$y,$z) = $SPG->TranslationVector($i+1);
#		print "   ($x, $y, $z)\n";
	}

	$App->print("\n");
	my $nSymmetryOperation = $SPG->nSymmetryOperation();
#	print "nSymmetryOperation: $nSymmetryOperation[n";
	for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
		my $symop = $SPG->SymmetryOperation($i+1);
#		print "   $symop\n";
	
		my ($x1,$y1,$z1,$t1,
			$x2,$y2,$z2,$t2,
			$x3,$y3,$z3,$t3)
			= $SPG->SymmetryOperationByMatrix($i+1);
#		print "      $x1 $y1 $z1  $t1\n";
#		print "      $x2 $y2 $z2  $t2\n";
#		print "      $x3 $y3 $z3  $t3\n";
	}

	$App->print("\n");
#Crystalクラス中の、原子の種類 AtomTypeクラスのリストをとる
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	$App->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;
		$App->print("   #$i1: $label : $atomname [$AtomicNumber] ($charge) ($AtomicMass)\n");
	}

	$App->print("\n");
#Crystalクラス中の、非対称単位中の原子 AsymmetricAtomSiteクラスのリストをとる
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
	$App->print("nAsymmetricAtomSite: $nAsymmetricAtomSite\n");
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
		my $label     = $AtomSiteList[$i]->Label();
		my $type      = $AtomSiteList[$i]->AtomName();
		my ($x,$y,$z) = $AtomSiteList[$i]->Position(1);
		my $occupancy = $AtomSiteList[$i]->Occupancy();
		my $iAtomType = $Crystal->FindiAtomType($type);

		$App->print("   $label ($type)[#$iAtomType]: ($x, $y, $z) [$occupancy]\n");
	}

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

	$App->print("\n");

	my $nLatticeX = int($MaximumDistance / $a1) + 1;
	my $nLatticeY = int($MaximumDistance / $a1) + 1;
	my $nLatticeZ = int($MaximumDistance / $a1) + 1;
	$App->print("Unit cell numbers to be examined: ($nLatticeX, $nLatticeY, $nLatticeZ)\n");

	$App->print("\n");
	$App->print("Inter-atomic distances\n");
	my $count = 0;
	my @iAtomListArray;
	my @CDistanceArray;
	for(my $ia = 0 ; $ia < $nAsymmetricAtomSite ; $ia++) {
		my $label0       = $AtomSiteList[$ia]->Label();
		my $type0        = $AtomSiteList[$ia]->AtomName();
		my ($x0,$y0,$z0) = $AtomSiteList[$ia]->Position(1);
		my $occupancy0   = $AtomSiteList[$ia]->Occupancy();
#		my $iAtomType0   = $Crystal->FindiAtomType($type);
		my $ia1 = $ia+1;
		my $k1;
		for(my $k = 0 ; $k < $nTotalExpandedAtomSite ; $k++) {
			my $id = $ExpandedAtomSiteList[$k]->IdAsymmetricAtomSite();
			if($id == $ia1) {
				$k1 = $k+1;
				last;
			}
		}
		$k1 = -1 if($k1 eq '');
		my @iAtomList;
		$iAtomListArray[$ia] = \@iAtomList;;

		for(my $iz = -$nLatticeZ ; $iz <= $nLatticeZ ; $iz++) {
		for(my $iy = -$nLatticeY ; $iy <= $nLatticeY ; $iy++) {
		for(my $ix = -$nLatticeX ; $ix <= $nLatticeX ; $ix++) {
			last if($count >= $MaximumNumber);
			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 $dis = $Crystal->GetInterAtomicDistance(
					$x0, $y0, $z0, $x+$ix, $y+$iy, $z+$iz );
				next if($dis > $MaximumDistance);
				
				my %Peak;
				$Peak{'sn'}     = $count+1;
				$Peak{'label0'} = $label0;
				$Peak{'ia1'}    = $ia1;
				$Peak{'k1'}     = $k1;
				$Peak{'label'}  = $label;
				$Peak{'i1'}     = $i1;
				$Peak{'ix'}     = $ix;
				$Peak{'iy'}     = $iy;
				$Peak{'iz'}     = $iz;
				$Peak{'dis'}    = $dis;
				$Peak{'x0'}     = $x0;
				$Peak{'y0'}     = $y0;
				$Peak{'z0'}     = $z0;
				$Peak{'x1'}     = $x+$ix;
				$Peak{'y1'}     = $y+$iy;
				$Peak{'z1'}     = $z+$iz;
				push(@CDistanceArray, \%Peak);
#				printf "   %04d: %4s(#%3d,\@%3d) - %4s(\@%3d)[%d,%d,%d]: %7.4f A  ---  "
#					."(%7.4f,%7.4f,%7.4f)-(%7.4f,%7.4f,%7.4f)\n",
#					$count+1, $label0, $ia1, $k1, $label, $i1, $ix, $iy, $iz, $dis,
#					$x0, $y0, $z0, $x+$ix, $y+$iy, $z+$iz;

				my $xx = $x+$ix;
				my $yy = $y+$iy;
				my $zz = $z+$iz;
				my $s = "$i $label $xx $yy $zz";
				push(@iAtomList, $s);

				$count++;
				if($count >= $MaximumNumber) {
					$App->print("Maximum number ($MaximumNumber) reached.\n");
					$App->print("Skip further calculation.\n");
					last;
				}
			}
		}
		}
		}
	}

	if($DoSort) {
#print("SORT\n");
		@CDistanceArray = &Sort(\@CDistanceArray);
#		@CDistanceArray = sort { $a->{'dis'} <=> $b->{'dis'} } @CDistanceArray;
#		@CDistanceArray = sort SortDistanceArrayFunc @CDistanceArray;
	}
	for(my $i = 0 ; $i < @CDistanceArray ; $i++) {
		my $peak = $CDistanceArray[$i];
		$App->printf("   %04d(%04d): %4s(#%3d,\@%3d) - %4s(\@%3d)[%2d,%2d,%2d]: %7.4f A  ---  "
			."(%7.4f,%7.4f,%7.4f)-(%7.4f,%7.4f,%7.4f)\n",
			$i+1, $peak->{'sn'}, $peak->{'label0'}, $peak->{'ia1'}, $peak->{'k1'},
			$peak->{'label'}, $peak->{'i1'}, $peak->{'ix'}, $peak->{'iy'}, $peak->{'iz'}, 
			$peak->{'dis'}, $peak->{'x0'}, $peak->{'y0'}, $peak->{'z0'}, 
			$peak->{'x1'}, $peak->{'y1'}, $peak->{'z1'});
	}

	$App->print("\n");
	$App->print("Inter-atomic Angles\n");
	my $count = 0;
	for(my $ia = 0 ; $ia < $nAsymmetricAtomSite ; $ia++) {
		next if($count >= $MaximumNumber);
		my $label0       = $AtomSiteList[$ia]->Label();
		my $type0        = $AtomSiteList[$ia]->AtomName();
		my ($x0,$y0,$z0) = $AtomSiteList[$ia]->Position(1);
		my $ia1 = $ia+1;
		my $iAtomList = $iAtomListArray[$ia];
		my $k1;
		for(my $k = 0 ; $k < $nTotalExpandedAtomSite ; $k++) {
			my $id = $ExpandedAtomSiteList[$k]->IdAsymmetricAtomSite();
			if($id == $ia1) {
				$k1 = $k+1;
				last;
			}
		}
		$k1 = -1 if($k1 eq '');

		for(my $i1 = 0 ; $i1 < @$iAtomList ; $i1++) {
			next if($count >= $MaximumNumber);
			my $s1 = $iAtomList->[$i1];
			my ($i1a, $label1, $x1, $y1, $z1) = (split(/\s/, $s1));
			my $i11 = $i1+1;
			my $dis01 = $Crystal->GetInterAtomicDistance(
				$x0, $y0, $z0, $x1, $y1, $z1);
#print "dis01($x0,$y0,$z0)-($x1,$y1,$z1): $dis01\n";
			next if($dis01 <= 0.0);
			for(my $i2 = 0 ; $i2 < @$iAtomList ; $i2++) {
				my $s2 = $iAtomList->[$i2];
				my ($i2a, $label2, $x2, $y2, $z2) = (split(/\s/, $s2));
				my $i21 = $i2+1;
				my $dis02 = $Crystal->GetInterAtomicDistance(
					$x0, $y0, $z0, $x2, $y2, $z2);
				next if($dis02 <= 0.0);

				my $dis12 = $Crystal->GetInterAtomicDistance(
					$x1, $y1, $z1, $x2, $y2, $z2);
				next if($dis12 <= 0.0);

				my $angle = $Crystal->GetInterAtomicAngle(
							$x0, $y0, $z0, 
							$x1, $y1, $z1,
							$x2, $y2, $z2 );
				$App->printf("   %04d: %4s(#%3d,\@%3d) - %4s(\@%3d) - %4s(\@%3d): "
					."%7.3f deg  ---  "
					."(%6.3f,%6.3f,%6.3f)-(%6.3f,%6.3f,%6.3f)-(%6.3f,%6.3f,%6.3f)\n",
					$count+1, 
					$label0, $ia1, $k1, $label1, $i11, $label2, $i21, 
					$angle,
					$x0, $y0, $z0, $x1, $y1, $z1, $x2, $y2, $z2);

				$count++;
				if($count >= $MaximumNumber) {
					$App->print("\nMaximum number ($MaximumNumber) reached.\n");
					$App->print("Skip further calculation.\n");
					last;
				}
			}
		}
	}

	$App->print("Calculation finished.\n")
}