package MoleculeObject;
use Exporter;
@ISA = qw(Exporter);

#公開したいサブルーチン
@EXPORT = qw();

use strict;
#use Math::Matrix;
use Math::MatrixReal;

use Sci qw($pi $a0 $torad $todeg);
use Sci::Algorism;
use Sci::Vector;

use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::Rietan;

#==================================================================
# 静的メンバー関数
#==================================================================

#==================================================================
# 変数取得関数
#==================================================================
sub Title           { return shift->{Title}; }
sub SetTitle        { my ($this,$t)=@_; return $this->{Title} = $t; }
sub SampleName      { return shift->{SampleName}; }
sub SetSampleName   { my ($this,$s)=@_; return $this->{SampleName} = $s; }
sub SetMoleculeName { my ($this,$name) = @_; return $this->{MoleculeName} = $name; }
sub MoleculeName    { my ($this) = @_; return $this->{MoleculeName}; }
sub SetCalculationCondition { my ($this,$c)=@_; return $this->{SetCalculationCondition} = $c; }
sub CalculationCondition    { my ($this)=@_; return $this->{SetCalculationCondition}; }
sub nAtomSite         { return shift->{nAtomSite}; }
sub SetnAtomSite      { my ($this,$n)=@_; return $this->{nAtomSite} = $n; }
sub nAtomType         { return shift->{nAtomType}; }
sub SetnAtomType      { my ($this,$n)=@_; return $this->{nAtomType} = $n; }
sub nZMatrix          { return shift->{nZMatrix}; }
sub SetnZMatrix       { my ($this,$n)=@_; return $this->{nZMatrix} = $n; }
sub pZMatrixList      { return shift->{RefZMatrixList}; }
sub ZMatrix           { my ($this,$i)=@_; return $this->{RefZMatrixList}->[$i]; }
sub SetnAtomType { my ($this,$n) = @_; return $this->{nAtomType} = $n; }
sub nAtomType    { my ($this) = @_; return $this->{nAtomType}; }

sub SetCAtomTypeList  
{
	my ($this, @atomlist) = @_;
	$this->{nAtomType} = @atomlist;
	return $this->{RefAtomTypeList} = \@atomlist;
}

sub GetCAtomTypeList
{
	my ($this) = @_;
	my $RefAtomList = $this->{RefAtomTypeList};
#print "Ref:$RefAtomList<br>\n";
	return () if($RefAtomList eq '');
	return @$RefAtomList;
}

sub GetCAtomType
{
	my ($this,$i) = @_;
	my $RefAtomList = $this->{RefAtomTypeList};
	my $atomtype =  $RefAtomList->[$i];
	if(!$atomtype) {
		return 0;
	}
	return $atomtype;
}

sub SetCAtomSiteList  
{
	my ($this, @atomlist) = @_;
	$this->{nAtomSite} = @atomlist;
	return $this->{RefAtomSiteList} = \@atomlist;
}

sub GetCAtomSiteList
{
	my ($this) = @_;
	my $RefAtomList = $this->{RefAtomSiteList};
#print "Ref:$RefAtomList<br>\n";
	return () if($RefAtomList eq '');
	return @$RefAtomList;
}

sub GetCAtomSite
{
	my ($this,$i) = @_;
	my $RefAtomList = $this->{RefAtomSiteList};
	my $atomsite =  $RefAtomList->[$i];
	if(!$atomsite) {
		return 0;
	}
	return $atomsite;
}

#How to build Molecule (Ex. from CIF.pm)
#	my $mol = new Molecule();
#	$mol->SetCAtomTypeList(
#		$this->GetCAtomTypeList() );
#	$mol->SetCAtomSiteList(
#		$this->GetCAtomSiteList() );

#==================================================================
# コンストラクタ、デストラクト
#==================================================================
sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	$this->Initialize();
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub Initialize
{
	my ($this) = @_;

	$this->{RefAtomTypeList} = [];
	$this->{nAtomType}       = 0;
	$this->{RefAtomSiteList} = [];
	$this->{nAtomSite}       = 0;
	$this->{RefZMatrixList}  = [];
	$this->{nZMatrix}        = 0;
	return $this;
}

#==================================================================
# 一般メンバー関数
#==================================================================
sub TranslateToCenter
{
	my ($this) = @_;
	my ($xc, $yc, $zc) = $this->CalculateCenterPosition();
	my $nSite = $this->nAtomSite();
	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $site = $this->GetCAtomSite($i);
#		last if(!$site);
		my $atomname = $site->AtomNameOnly();
		my ($AtmNum, $name, $charge) = AtomType::GetAtomInformation($atomname);
		my ($x, $y, $z) = $site->Position();
		$site->SetPosition($x-$xc, $y-$yc, $z-$zc);
	}
}

sub CalculateCenterPosition
{
	my ($this) = @_;
	
	my ($xc, $yc, $zc) = (0, 0, 0);
	my $nSite = $this->nAtomSite();
	for(my $i = 0 ; $i < $nSite ; $i++) {
		my $site = $this->GetCAtomSite($i);
		my $atomname = $site->AtomNameOnly();
		my ($AtmNum, $name, $charge) = AtomType::GetAtomInformation($atomname);
		my ($x, $y, $z) = $site->Position();
		$xc += $x;
		$yc += $y;
		$zc += $z;
	}
	return ($xc / $nSite, $yc / $nSite, $zc / $nSite);
}

sub AddZMatrix
{
	my ($this, $atom, $R, $Angle1, $Angle2, $iAtom2, $iAtom3, $iAtom4, 
		$idR, $idAngle1, $idAangle2) = @_;

	my %Z = (
		AtomName => $atom, 
		R        => $R, 
		Angle1   => $Angle1, 
		Angle2   => $Angle2, 
		iAtom2   => $iAtom2, 
		iAtom3   => $iAtom3, 
		iAtom4   => $iAtom4, 
		idR      => $idR, 
		idAngle1 => $idAngle1, 
		idAngle2 => $idAangle2
		);
	$this->AddAtomType($atom, 1);
	$this->{RefZMatrixList}->[$this->{nZMatrix}] = \%Z;
	return $this->{nZMatrix}++;
}

sub ZMatrixToXYZ
{
	my ($this, $TranslateToCenter, $IsPrint) = @_;
	$IsPrint = 0 if(!defined $IsPrint);

	my %iAtom;
	my (@x, @y, @z);
	for(my $i = 0 ; $i < $this->nZMatrix() ; $i++) {
		my $Z = $this->ZMatrix($i);
		
		my $atom  = $Z->{AtomName};
		$iAtom{$atom}++;
		my $label = $atom . $iAtom{$atom};
		my $R  = $Z->{R};
		my $angle1 = $Z->{Angle1};
		my $angle2 = $Z->{Angle2};

		my ($x, $y, $z);
		if($i == 0) {
			($x, $y, $z) = (0, 0, 0);
		}
		elsif($i == 1) {
			($x, $y, $z) = ($R, 0, 0);
		}
		elsif($i == 2) {
			$x = $R * cos($angle1 * $torad);
			$y = $R * sin($angle1 * $torad);
			$z = 0.0;
		}
		else {
			my $alpha = $torad * $angle1;
			my $beta  = $torad * $angle2;
			my $sina  = sin($alpha);
			my $cosa  = cos($alpha);
			my $sinb  = sin($beta);
			my $cosb  = cos($beta);
			my $j = $Z->{iAtom2}-1;
			my $k = $Z->{iAtom3}-1;
			my $l = $Z->{iAtom4}-1;
print "$i-$j-$k-$l: $R: $angle1: $angle2\n" if($IsPrint);

			my @a  = Vector::Difference($x[$k], $y[$k], $z[$k], $x[$j], $y[$j], $z[$j]);
			my @b  = Vector::Difference($x[$l], $y[$l], $z[$l], $x[$j], $y[$j], $z[$j]);
			my @c  = Vector::OuterProduct(@b, @a);
#			my @c  = Vector::OuterProduct(@a, @b);

			my $a  = sqrt(Vector::InnerProduct(@a, @a));
			my $b  = sqrt(Vector::InnerProduct(@b, @b));
			my $c  = sqrt(Vector::InnerProduct(@c, @c));

			my $ab  = Vector::InnerProduct(@a, @b);

printf "a=(%6.3f, %6.3f, %6.3f)\n", @a if($IsPrint);
printf "b=(%6.3f, %6.3f, %6.3f)\n", @b if($IsPrint);
printf "c=(%6.3f, %6.3f, %6.3f)\n", @c if($IsPrint);
#my $ac  = Vector::InnerProduct(@a, @c);
#my $bc  = Vector::InnerProduct(@b, @c);
#print "a*c=$ac  b*c=$bc\n" if($IsPrint);

#my $ab  = Vector::InnerProduct(@a, @b);
#my $angleab = $todeg * acos($ab / $a / $b);
#print "angle(a-b)=$angleab deg.\n" if($IsPrint);

			my $B = $R * $a * $sina * $cosb / $c;
			my $A = ($R * $a * $cosa - $B * $ab) / $a / $a;
			my $C  = sqrt($R*$R - $A*$A * $a*$a - $B*$B * $b*$b - 2.0 * $A * $B * $ab) / $c;
			$C = -$C if($angle2 < 0.0);

			($x, $y, $z) = ($x[$j] + ($A*$a[0] + $B*$b[0] + $C*$c[0]), 
							$y[$j] + ($A*$a[1] + $B*$b[1] + $C*$c[1]), 
							$z[$j] + ($A*$a[2] + $B*$b[2] + $C*$c[2]));
		}
		($x[$i], $y[$i], $z[$i]) = ($x, $y, $z);
		$this->AddAtomSite($label, $atom, $x, $y, $z, 1.0);
	}

	$this->TranslateToCenter() if($TranslateToCenter);
}

sub FindNearestEquivalentFractionCoordinate
{
	my ($this, $x2, $x) = @_;
	my $xret = $x2;
	my $mind = abs($x2 - $x);
	if($mind > abs($x2+1.0 - $x)) {
		$xret = $x2+1.0;
		$mind = abs($x2+1.0 - $x);
	}
	if($mind > abs($x2-1.0 - $x)) {
		$xret = $x2-1.0;
		$mind = abs($x2-1.0 - $x);
	}
	return $xret;
}

sub FractionalToCartesian
{
	my ($this, $x,$y,$z) = @_;
	
	my $a11 = $this->{a11};
	my $a12 = $this->{a12};
	my $a13 = $this->{a13};
	my $a21 = $this->{a21};
	my $a22 = $this->{a22};
	my $a23 = $this->{a23};
	my $a31 = $this->{a31};
	my $a32 = $this->{a32};
	my $a33 = $this->{a33};
	my $xc = $a11*$x + $a21*$y + $a31*$z;
	my $yc = $a12*$x + $a22*$y + $a32*$z;
	my $zc = $a13*$x + $a23*$y + $a33*$z;
	return ($xc,$yc,$zc);
}

sub FindAtomTypeByName
{
	my ($this, $atomname) = @_;
	$atomname =~ s/([\d+-]+)//;

	my $atomlistref = $this->{RefAtomTypeList};
	for(my $i = 0 ; $i < @$atomlistref ; $i++) {
		my $type = $atomlistref->[$i];
		my $name = $type->AtomNameOnly();
#print "name: $atomname <=> $name\n";
		return $type if(uc $name eq uc $atomname);
	}
	return undef;
}

sub AddAtomType
{
	my ($this, $atomname, $CheckRegistered) = @_;
	$CheckRegistered = 1 if(!defined $CheckRegistered);

	my $atomlistref = $this->{RefAtomTypeList};
	if($atomlistref) {
		if($CheckRegistered and $this->FindAtomTypeByName($atomname)) {
			return;
		}
		my $at = new AtomType();
		$at->SetAtomName($atomname);
		push(@$atomlistref, $at);

		$this->{"nAtomType"} = scalar @$atomlistref;
	}
	else {
		my @atomlist;
		my $at = new AtomType();
		$at->SetAtomName($atomname);
		push(@atomlist, $at);

		$this->{nAtomType}       = @atomlist;
		$this->{RefAtomTypeList} = \@atomlist;
	}

#my $ref = $this->{RefAtomTypeList} ;
#print "ref: ", $ref->[$this->{nAtomType}-1]->AtomName(), "<BR>";
	return $this->{nAtomType};
}

sub AddAtomSite
{
	my ($this, $label, $atomname, $x, $y, $z, $occ, $fx, $fy, $fz) = @_;
#print "AddAtomSite: F=$fx, $y, $fz\n";
#	return $this->{nAtomSite} if($atomname eq '');

	my $atomlistref = $this->{RefAtomSiteList};
	if(!$atomlistref) {
		$atomlistref = [];
	}

	my $at = new AtomSite();
	$at->SetLabel($label);
	$at->SetAtomName($atomname);
	$at->SetPosition($x, $y, $z);
	$at->SetForce($fx, $fy, $fz) if(defined $fz);
	$at->SetOccupancy($occ);
	$at->SetId(@$atomlistref + 1);

	push(@$atomlistref, $at);
	$this->{nAtomSite} = scalar @$atomlistref;

	$this->AddAtomType($atomname) if($atomname ne '');
	return $this->{nAtomSite};
}

sub GetInterAtomicDistance
{
	my ($this, $x0, $y0, $z0, $x1, $y1, $z1) = @_;;
	my $dx = $x1 - $x0;
	my $dy = $y1 - $y0;
	my $dz = $z1 - $z0;
	my $r2 = $dx*$dx + $dy*$dy + $dz*$dz;
	return sqrt($r2);
}

sub GetInterAtomicAngle
{
	my ($this, $x0, $y0, $z0, $x1, $y1, $z1,
		   $x2, $y2, $z2) = @_;;

	my $dis01 = $this->GetInterAtomicDistance($x0,$y0,$z0,$x1,$y1,$z1);
	return 0.0 if($dis01 == 0.0);
	my $dis02 = $this->GetInterAtomicDistance($x0,$y0,$z0,$x2,$y2,$z2);
	return 0.0 if($dis02 == 0.0);
#print "dis: $dis01, $dis02\n";
	my $dx01 = $x1 - $x0;
	my $dy01 = $y1 - $y0;
	my $dz01 = $z1 - $z0;
	my $dx02 = $x2 - $x0;
	my $dy02 = $y2 - $y0;
	my $dz02 = $z2 - $z0;
	my $ip  = $dx01*$dx02 + $dy01*$dy02 + $dz01*$dz02;

	my $cosa = $ip / $dis01 / $dis02;
#print "cos: $cosa<br>";
	my $angle = Sci::dacos($cosa);
#print "cos: $cosa $angle<br>";
	$angle = 360.0 - $angle if($angle > 180.0001);

	return $angle;
}

#iAtomTypeは1から始まる
sub FindiAtomType
{
	my ($this, $atomname) = @_;

	$atomname =~ s/^(\D+)(\d.*)$/$1/;
	$atomname = uc $atomname;
	my @AtomTypeList = $this->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name = $AtomTypeList[$i]->AtomName();
		$name =~ s/^(\D+)(\d.*)$/$1/;
		return $i+1 if($atomname eq uc $name);
	}
	return -1;
}

sub FillAtomTypeData
{
	my ($this) = @_;
	my $ref = $this->{"RefAtomTypeList"};
	for(my $i = 0 ; $i < @$ref ; $i++) {
		$ref->[$i]->FillAtomTypeData();
	}
}

# $s = sin Q / lambda (in nm)
sub AtomicScatteringFactor
{
	my ($this, $s) = @_;
	return $this->asf($s);
}

# $s = sin Q / lambda (in nm)
sub asfElectron
{
	my ($this, $iAtomType, $s) = @_;
	$s = 1.0e-5 if($s <= 0.0);
	my $s2 = $s * $s;
	my $k  = 1.0e-2 * 8.0 * $pi * $pi / ($a0 * 1.0e9);
	my $Z  = $this->{"ASFDCAtomNameZ$iAtomType"};

	my $asfXray = $this->asf($iAtomType, $s);
#print "s=$s, i=$iAtomType, Z=$Z, a=$asfXray\n";
#exit;
	my $val = $k * ($Z - $asfXray) / $s2;
#print "va [$iAtomType]=$val = $k * ($Z - $asfXray) / $s2\n";
	return $val;
}

sub ReadASF
{
	my ($this, $Source) = @_;

	my @AtomTypeList = $this->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $Rietan = new Rietan;
		my $AtomName;
		my $AtomNameOnly = $AtomTypeList[$i]->AtomNameOnly();
		my $Charge   = $AtomTypeList[$i]->Charge();
		if($Charge == 1) {
			$AtomName = "$AtomNameOnly+";
		}
		elsif($Charge == -1) {
			$AtomName = "$AtomNameOnly-";
		}
		elsif($Charge > 1) {
			$AtomName = "$AtomNameOnly+";
		}
		elsif($Charge < -1) {
			$Charge = -$Charge;
			$AtomName = "$AtomNameOnly-";
		}
#print "a1: $AtomName\n";
		my @a = $Rietan->ReadASFParameters($AtomName);
		if(@a < 9) {
			$AtomName = $AtomTypeList[$i]->AtomNameOnly();
#print "a2: $AtomName\n";
			@a = $Rietan->ReadASFParameters($AtomName);
		}
		if(@a < 9) {
			print "Error: Can not read ASFDC for [$AtomName].\n";
			next;
		}
#print "array: ", join(',', @a), "\n";
#print "i=$i: $AtomName\n";
		$this->{"ASFDCRed$i"}              = $Rietan;
		$this->{"ASFDCRed$AtomName"}       = $Rietan;
		$this->{"ASFDCAtomName$i"}         = $Rietan->{ASFDCAtomName};
		$this->{"ASFDCAtomName$AtomName"}  = $Rietan->{ASFDCAtomName};
		my $AtomicNumber = $AtomTypeList[$i]->AtomicNumber();
		$this->{"ASFDCAtomNameZ$i"}            = $AtomicNumber;
		$this->{"ASFDCAtomNameZ$AtomNameOnly"} = $AtomicNumber;
print "ASFDC AtomName [$i]: $Rietan->{ASFDCAtomName} [Zi=$AtomicNumber]\n";
#print "   ao:$AtomNameOnly: $i\n";
	}
	$this->{ASFDCRed} = 1;
}

# $s = sin Q / lambda (in nm)
sub asf
{
	my ($this, $iAtomType, $s, $Source) = @_;
	$Source = 'Cu' if(!defined $Source);
	$this->ReadASF($Source) 	if(!defined $this->{ASFDCRed});

	my $Rietan = $this->{"ASFDCRed$iAtomType"};
	my $pP     = $Rietan->{pASFDCParameters};
	my ($A1,$B1,$A2,$B2,$A3,$B3,$A4,$B4,$C, $b,
		$Crdf1,$Crdf2, $Fedf1,$Fedf2, $Codf1,$Codf2, $Cudf1,$Cudf2, 
		$Modf1, $Modf2, $Agdf1,$Agdf2, $Kbdf1, $Kbdf2) = @$pP;
#	if($Source eq 'n') {
#		return $b;
#	}

	my $s2 = 0.01 * $s * $s;
	my $asf = $A1*exp(-$B1*$s2)+$A2*exp(-$B2*$s2)+$A3*exp(-$B3*$s2)+$A4*exp(-$B4*$s2)+$C;
	if($Source eq 'Cr' and defined $Crdf2) {
		$asf = cplx($asf + $Crdf1, $Crdf2);
	}
	elsif($Source eq 'Fe' and defined $Fedf2) {
		$asf = cplx($asf + $Fedf1, $Fedf2);
	}
	elsif($Source eq 'Co' and defined $Codf2) {
		$asf = cplx($asf + $Codf1, $Codf2);
	}
	elsif($Source eq 'Cu' and defined $Cudf2) {
#print "d=$Cudf1,$Cudf2\n";
		$asf = cplx($asf + $Cudf1, $Cudf2);
	}
	elsif($Source eq 'Mo' and defined $Modf2) {
		$asf = cplx($asf + $Modf1, $Modf2);
	}
	elsif($Source eq 'Ag' and defined $Agdf2) {
		$asf = cplx($asf + $Agdf1, $Agdf2);
	}
	elsif($Source eq 'Kb' and defined $Kbdf2) {
		$asf = cplx($asf + $Kbdf1, $Kbdf2);
	}

	return $asf;
}
1;