#===============================================
# ATOMS
#===============================================
package ATOMS;
use Exporter;
use Common;
@ISA = qw(Exporter Common);

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

use strict;
use File::Basename;

use JFile;
use Sci::Science;

use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::CIF;
use Crystal::RietanFP;

#============================================================
#　静的メンバー関数
#============================================================
sub GetSTRFile { return Deps::ReplaceExtension($_[0], ".str"); }
sub GetPDBFile { return Deps::ReplaceExtension($_[0], ".pdb"); }

#============================================================
#　変数取得関数など
#============================================================
sub SetFileName
{
	my ($this, $path) = @_;
	return $this->{'FileName'} = $path;
}
sub FileName
{
	my ($this) = @_;
	return $this->{'FileName'};
}
sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}
sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

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

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
#print "In ATOMS::CheckFileType for $path\n";
	if($ext =~ /^\.str$/i) {
		return "ATOMS str file";
	}
	elsif($ext =~ /^\.pdb$/i) {
		return "ATOMS pdb file";
	}
	return undef;
}

sub Read
{
	my ($this, $filename) = @_;
	my $FileType = $this->{'FileType'} = ATOMS::CheckFileType($filename);
	return undef unless($FileType);

	$this->SetFileName($filename);
	if($FileType eq "ATOMS str file") {
		return $this->ReadSTRFile($filename);
	}
	elsif($FileType eq "ATOMS pdb file") {
		return $this->ReadPDBFile($filename);
	}
	return undef;
}

sub ReadSTRFile
{
	my ($this, $filename, $IsPrint) = @_;
	$IsPrint = 1 unless(defined $IsPrint);

$this->print("Read ATOMS str file\n") if($IsPrint);

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($filename);
	my $SampleName = $filebody;

	my $crystal = new Crystal();
	my $in = new JFile($filename, "r");
	return undef unless($in);

	my $line = $in->SkipTo("TITL");
	$line =~ s/TITL\s+//;
	Utils::DelSpace($line);
	$SampleName = $line if($line !~ /^\d+-ICSD$/i);
	$SampleName = $line if($SampleName eq '');
	$crystal->SetCrystalName($SampleName);
	$crystal->SetSampleName($SampleName);
$this->print("SampleName: $SampleName\n") if($IsPrint);

	$line = $in->SkipTo("GNFP1");
	my ($head, $a, $b, $c, $alpha, $beta, $gamma)
		= Utils::Split("\\s+", $line);
	$a += 0;
	$b += 0;
	$c += 0;
	$alpha += 0;
	$beta  += 0;
	$gamma += 0;
	$crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
$this->print("Latt: $a $b $c   $alpha $beta $gamma\n") if($IsPrint);

	$line = $in->SkipTo("CSPGP");
	$line =~ s/^CSPGP\s+//;
	my $SPGName = Utils::DelSpace($line);
	$line = $in->SkipTo("CNSPG");
	$line =~ s/^CNSPG\s+//;
	my $iSPG    = Utils::DelSpace($line);
	$crystal->SetSpaceGroup($SPGName, $iSPG);
$this->print("SPG: $SPGName ($iSPG)\n") if($IsPrint);

	my $nAsymmetricAtom = 0;
#$this->print("nAsymmetricAtom: $nAsymmetricAtom\n") if($IsPrint);
	my %AtomCount;
	for(my $i = 0 ;  ; $i++) {
		$line = $in->SkipTo("ATBAS1");
		last unless(defined $line);
		my ($head, $AtomName, $atnum, $x, $y, $z)
			= Utils::Split("\\s+", $line);
		last unless(defined $z);

		$AtomName =~ s/\d+//g;
		$AtomCount{$AtomName}++;
		my $Label = $AtomName . $AtomCount{$AtomName};
		$crystal->AddAtomType($AtomName);
		$crystal->AddAtomSite($Label, $AtomName, $x, $y, $z, 1.0);
$this->print("$i: $Label: $AtomName ($x, $y, $z)\n") if($IsPrint);
	}

	my $Rietan = new Rietan();
#iCenterをSPGRAから読み取る
	my $SPG = $Rietan->ReadSpaceGroup($iSPG, 1);
	my $SPG = new SpaceGroup();
	$SPG->SetSPGName($SPGName);
	$SPG->SetiSPG($iSPG);
	my $iSet = 1;
	$SPG->SetiSet($iSet);
	$SPG->AnalyzeTranslation();
	$in->rewind();
	for(my $i = 0 ;  ; $i++) {
		my $SymOp = $in->SkipTo("SYMXYZ");
		last unless(defined $SymOp);
		$SymOp =~ s/SYMXYZ\s*//;
		Utils::DelSpace($SymOp);
		$SymOp = lc $SymOp;
		$SPG->AddSymmetryOperation($SymOp);
$this->print("symop[$i]: $SymOp<br>\n");
#		my @SymOpXYZ = split(/,\s*/, $SymOp);
#		my @XMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[0]);
#		my @YMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[1]);
#		my @ZMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[2]);
	}
$this->print("iCenter: ", $Rietan->iCenter(), "\n");
	$SPG->ExpandSymmetryOperation($Rietan->iCenter());

	$crystal->SetCSpaceGroup($SPG) if(defined $SPG);

	$in->Close();

	$crystal->ExpandCoordinates();

	return $crystal;
}

sub ReadPDBFile
{
	my ($this, $filename, $IsPrint) = @_;
	$IsPrint = 1 unless(defined $IsPrint);

$this->print("Read ATOMS pdb file\n") if($IsPrint);

	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($filename);
	my $SampleName = $filebody;

	my $crystal = new Crystal();
	my $in = new JFile($filename, "r");
	return undef unless($in);

	my $line = $in->SkipTo("TITLE");
	my ($title) = ($line =~ /TITLE\s+(\S.*)\s*$/);
	$SampleName = $line if($title !~ /^\d+-ICSD$/i);
	$SampleName = $line if($SampleName eq '');
	$crystal->SetCrystalName($SampleName);
	$crystal->SetSampleName($SampleName);
$this->print("SampleName: $SampleName\n") if($IsPrint);

	$line = $in->SkipTo("CRYST1");
	my ($a, $b, $c, $alpha, $beta, $gamma, $SPGName, $iSPG)
		= ($line =~ /CRYST1\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(\S+)\s+(.+)\s+(\d+)\s*$/);
	$a += 0;
	$b += 0;
	$c += 0;
	$alpha += 0;
	$beta  += 0;
	$gamma += 0;
$this->print("Latt: $a $b $c   $alpha $beta $gamma\n") if($IsPrint);
	$crystal->SetLatticeParameters($a, $b, $c, $alpha, $beta, $gamma);
$this->print("SPG: $SPGName ($iSPG)\n") if($IsPrint);
#	$crystal->SetSpaceGroup($SPGName, $iSPG);
	$crystal->SetSpaceGroup("P1", 1);

	$line = $in->SkipTo("SCALE1");
	my ($head1, $a11,$a12,$a13,$t1) = Utils::Split("\\s+", $line);
	$line = $in->SkipTo("SCALE1");
	my ($head2, $a21,$a22,$a23,$t2) = Utils::Split("\\s+", $line);
	$line = $in->SkipTo("SCALE1");
	my ($head3, $a31,$a32,$a33,$t3) = Utils::Split("\\s+", $line);

	my %AtomCount;
	for(my $i = 0 ;  ; $i++) {
		$line = $in->SkipTo("ATOM");
		last unless(defined $line);
		my ($head, $idx, $AtomName, $num, $xc, $yc, $zc)
			= Utils::Split("\\s+", $line);
		last unless(defined $zc);

		$AtomName =~ s/\d+//g;
		$AtomCount{$AtomName}++;
		my $Label = $AtomName . $AtomCount{$AtomName};
#		my $x = $rx;
#		my $y = $ry;
#		my $z = $rz;
#		my ($x, $y, $z) = $crystal->CartesianToFractional($xc, $yc, $zc);
		my $x = $a11*$xc + $a12*$yc + $a13*$zc + $t1;
		my $y = $a21*$xc + $a22*$yc + $a23*$zc + $t1;
		my $z = $a31*$xc + $a32*$yc + $a33*$zc + $t1;
		$x = Utils::Round($x+0.000099, 3);
		$y = Utils::Round($y+0.000099, 3);
		$z = Utils::Round($z+0.000099, 3);
		$crystal->AddAtomType($AtomName);
		$crystal->AddAtomSite($Label, $AtomName, $x, $y, $z, 1.0);
$this->print("$i: $Label: $AtomName ($x, $y, $z)\n") if($IsPrint);
	}

#	my $Rietan = new Rietan();
#	my $SPG = $Rietan->ReadSpaceGroup($iSPG, 1);
#	$crystal->SetCSpaceGroup($SPG) if(defined $SPG);

	$in->Close();

	$crystal->ExpandCoordinates();

	return $crystal;
}

1;