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

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

use strict;
use File::Basename;

use Deps;
use Utils;
use Sci::Science;

use ProgVars;
#use Crystal::MyUtility;
use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::CIF;

#===============================================
# スクリプト大域変数
#===============================================
my $LF = "<br>\n";
my $DirectorySeparator = "\\";


#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $Program    = basename($0);
my $ProgramDir = ProgVars->ProgramDir();
my $LDDir      = ProgVars->LDDir();
my $MXDDir     = ProgVars->MXDDir();
my $MXDOrtoDir = ProgVars->MXDOrtoDir();
my $MXDExecDir = ProgVars->MXDExecDir();
my $SX1DBPath  = ProgVars->SX1DBPath();

BEGIN {
}

sub new
{
	my ($module) = @_;
	my $this = {};
	bless $this;
	return $this;
}

sub DESTROY
{
	my $this = shift;
}

sub SetSampleName
{
	my ($this, $name) = @_;
	return $this->{'SampleName'} = $name;
}

sub SampleName
{
	my ($this) = @_;
	return $this->{'SampleName'};
}

sub GetBusingParameter
{
	my ($this,$atomname) = @_;
	
	unless(open(IN, "<$SX1DBPath")) {
		print "Error in MXD.pm::GetBusingParameter: "
			."Can not read DB [$SX1DBPath].<br>\n";
		return ('');
	}
	<IN>;
	while(<IN>) {
		my ($name, $mass, $charge, $ai, $bi, $ci, $rad)
			= split(/\s+/, $_);
#print "name: $name, $atomname ($mass, $charge, $ai, $bi) <br>\n";
		if(uc $name eq uc $atomname) {
			return ($mass, $charge, $ai, $bi, $ci, $rad);
		}
	}
	
	close(IN);
	return ('');
}

sub MakeInputFile
{
	my ($this, $Crystal, $fname, $IonRadius, $Grange, $pModifiedParams) = @_;

	my $OutFile = "$fname.pme";
	my $out = JFile->new($OutFile, 'w');
	if(!$out) {
		print "Can't write to $OutFile.\n\n";
		return;
	}

print "IonRadius=$IonRadius\n";
print "Grange=$Grange\n";

	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $Crystal->SPGName();
#	$SPGName =~ s/[\s']//g;
#print OUT "$SPGName\n";
	my @SymmetryOperation  = $SPG->ExpandSymmetryOperation();
	my $nSymmetryOperation = @SymmetryOperation;
print "nSymmetry=$nSymmetryOperation\n";

	my @AtomTypeList    = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
print "nAtom=$nAtomType\n";
	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParametersByOutputMode(0);
	my @AsymmetricAtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite    = @AsymmetricAtomSiteList;
print "nAsymmetricAtomSite=$nAsymmetricAtomSite\n";
	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteList();
#	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nExpandedAtomSite    = @ExpandedAtomSiteList;
print "nExpandedAtomSite=$nExpandedAtomSite\n";

#O2 Si1                                                                  
	$out->print("$fname\n");
#  1.20  4.00  # Ion radius, max Ghkl in k space
	$out->printf(" %4.2f %5.2f\n", $IonRadius, $Grange);
#  6  2  0     # nSymmetry except inversion, nAsymmetric sites, Inversion (0: non-centrosym, 1: centrosym)
	$out->printf("%3d%3d  %d\n", 1, $nExpandedAtomSite, 0);
#   4.9210   4.9210   5.4160   90.000   90.000  120.000
	$out->printf(" %8.4f %8.4f %8.4f %8.4f %8.4f %8.4f\n", $a, $b, $c, $alpha, $beta, $gamma);
#   0.000000 1 0   0.000000 2 0   0.000000 3 0
	$out->print("   0.000000 1 0   0.000000 2 0   0.000000 3 0\n");
#Si1         4.0000   1.0000   0.5302   0.0000   0.3330
#O1         -2.0000   1.0000   0.4151   0.1476   0.1194
	for(my $i = 0 ; $i < $nExpandedAtomSite ; $i++) {
		my $atom = $ExpandedAtomSiteList[$i];
		my $atomname0   = $atom->AtomName();
		my $atomname    = $atom->AtomNameOnly();
		my $iAtom       = $Crystal->FindiAtomType($atomname);
		my $occ         = $atom->Occupancy();
		my ($x, $y, $z) = $atom->Position(1);
#		my $mult        = $atom->Multiplicity();
#		$mult = 1 if($IsExpandCoordinate or $IsChooseRandomly);
		my $charge0     = $atom->Charge();

		my ($Mass, $charge, $ai, $bi, $ci, $rad)
			= $this->GetBusingParameter($atomname);
#		$Mass = AtomType::GetAtomInformationByKey($atomname, "原子量") if($Mass eq '');
#		$Mass =~ s/\(.*?\)//g;
#		$Mass = eval($Mass);
#		$ai = 0.9  if($ai == 0);
#		$bi = 0.08 if($bi == 0);
#		$ci = 0.0  if($ci eq '');

		my ($zc, $sign) = ($charge0 =~ /([\d\.]+)([+\-])/);
		if($sign eq '-') {
			$charge0 = -$zc;
		}
		$charge = $charge0 if($charge0 ne '');

		my $iMP;
		if($pModifiedParams) {
			for(my $j = 0 ; $j < @$pModifiedParams ; $j++) {
				if($atomname eq $pModifiedParams->[$j]{Name}) {
					$iMP = $j;
					last;
				}
			}
#print "$atomname => index of Bader charge atom: $iMP\n";
		}
		if(defined $iMP) {
			my $p = $pModifiedParams->[$iMP];
			if(defined $p->{Charge}) {
				$charge = $p->{Charge};
			}
			if(defined $p->{ai}) {
				$ai = $p->{ai};
			}
			if(defined $p->{bi}) {
				$bi = $p->{bi};
			}
			if(defined $p->{ci}) {
				$ci = $p->{ci};
			}
			if(defined $p->{rad}) {
				$rad = $p->{rad};
			}
		}
#Si1         4.0000   1.0000   0.5302   0.0000   0.3330
		$out->printf("%-10s %7.4f %8.4f %8.4f %8.4f %8.4f\n", 
			$atomname0, $charge, $occ, $x, $y, $z);
	}
	close(OUT);

	return 1;
}

#	for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
#		my $SymOp = $SymmetryOperation[$i];
#		Utils::DelSpace($SymOp);
#
#		my @SymOpXYZ = split(/,\s*/, $SymOp);
#		my @XMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[0]);
#		my @YMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[1]);
#		my @ZMatrix = SpaceGroup::SymOpToMatrix($SymOpXYZ[2]);
#		printf OUT "%7.4f %7.4f %7.4f %7.4f %7.4f %7.4f "
#			."%7.4f %7.4f %7.4f\n",
#			$XMatrix[0], $XMatrix[1], $XMatrix[2],
#			$YMatrix[0], $YMatrix[1], $YMatrix[2],
#			$ZMatrix[0], $ZMatrix[1], $ZMatrix[2];
#		printf OUT "%7.4f %7.4f %7.4f\n",
#			eval($XMatrix[3]), eval($YMatrix[3]), 
#			eval($ZMatrix[3]);
#	}

1;