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

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

use strict;

use JFile;
use Crystal::Crystal;

BEGIN {
}

#===============================================
# スクリプト大域変数
#===============================================
my $Anions = "(O|S|Se|F|Cl|I)";

#===============================================
# パス（読み込みDB）
# Web関係変数
# CGI の仮想パス
# プログラム名
#===============================================
my $LF = "\n";
my $ProgramDir = "d:\\Programs";
my $VESTADir   = Deps::MakePath($ProgramDir, "Perl",   0);
   $VESTADir   = Deps::MakePath($VESTADir, "VESTA", 0);

my $thisAtomInfFile = "All.vesta";
   $thisAtomInfFile = Deps::MakePath($VESTADir, $thisAtomInfFile, 0);

#============================================================
#　変数等取得、再定義関数
#============================================================
sub VESTAAtomInfFile { return $thisAtomInfFile; }

#============================================================
#　継承クラスで定義しなおす関数
#============================================================

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

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　メンバー関数
#============================================================
sub GetAtomInformation
{
	my ($this, $name) = @_;
	if(!defined $this->{pAtomInfo}) {
		if(!$this->ReadAtomInformation()) {
			return undef;
		}
	}
	my $inf =  $this->{pAtomInfo}->{$name};
	if(!$inf) {
		$inf =  $this->{pAtomInfo}->{'H'};
	}
	return $inf;
}

sub ReadAtomInformation
{
	my ($this) = @_;
	my %AtomInfo;
	my $in = new JFile;
	if(!$in->Open($thisAtomInfFile, "r")) {
		return 0;
	}
	$in->SkipTo("ATOMT");
	while(1) {
		my $line = $in->ReadLine();
		last if($line + 0 == 0);
		Utils::DelSpace($line);
		my ($number, $name) = Utils::Split("\\s+", $line);
		$line =~ s/^\s*\d+\s+\w+\s+//;
		$AtomInfo{$name} = $line;
#print "[$name]: $AtomInfo{$name}\n";
	}
	$in->Close();
	$this->{pAtomInfo} = \%AtomInfo;
	return 1;
}


#ISURF
#render_sign iso_type sec_type isolevel sec_max sec_min
#tex_color_sign tex_max tex_min
#
#render_sign: 正負何れの領域を描くか; 0=両方, 1=正のみ, 2=負のみ
#iso_type: Isosurfaceのスタイル; 0=fill, 1=wireframe, 2=dotted
#sec_type: 断面(section)のスタイル
#isolevel: Isosufrace level
#sec_max: Sectionのsaturation level (max)
#sec_min: Sectionのsaturation level (min)
#tex_color_sign: 3D texture(静電ポテンシャル等)の色 1又は-1
#tex_max: Textureのsaturation level(max)
#tex_min: Textureのsaturation level(min)
#
#Isosurface levelの初期値はデフォルトが1, データの全ての値が1以上または1
#以下のときは最大値と最小値の中間に設定しています。VASPデータをisosurface
#用のデータとして読み込んだ場合、格子体積で割り算しています(ファイル名が
#LOCCAR, ELFCAR以外の場合だけ例外)。
#
#sec_max, sec_minの値もisosurface levelの値と同様、実データに対応する値で
#指定します。デフォルトでは実データの最大値の最小値に設定されます。
#
#tex_max, tex_minの値は0 - 1に規格化された値で、isosurface上のtextureデー
#タの最大値と最小値がそれぞれ1, 0に対応します。
#
#sec_type: 上位3ビットは2D data displayウィンドウで使われるもの
#//--- Section type flags
##define TYPE_COLOR_ABS   1   // 絶対値に基づいて彩色
##define TYPE_RECURSIVE   2   // Assign colors recursively
##define TYPE_LOGARITHTIC 4   // not used
##define TYPE_C_M_Y       8   // C-M-Y color mode
##define TYPE_MONOCHROME  16  // gray scale mode
##define TYPE_POLYGON     32  // fill polygons of sections
##define TYPE_CONTOUR     64  // draw contour lines
##define TYPE_GRIDEDGE    128 // draw mesh lines
sub WriteVESTAFileFromCrystal
{
	my ($this, $Crystal, $DensityFile, $path, 
		$BondLength, $Boundary, $BeyondBoundary,
		$DensityMin, $DensityMax, $DensityIsoLevel,
		%hash
		) = @_;
	$BondLength     = 3.0 if($BondLength <= 0);
	$BeyondBoundary = 3.0 if(!defined $BeyondBoundary);
	my ($LBoundary, $UBoundary) = split(/:/, $Boundary);
	$LBoundary       = 0.0 if(!defined $LBoundary);
	$UBoundary       = 1.0 if(!defined $UBoundary);
	$DensityMin      = 0.0 if(!defined $DensityMin);
	$DensityMax      = 1.0 if(!defined $DensityMax);
	$DensityIsoLevel = 1.0 if(!defined $DensityIsoLevel);

	if(!$this->ReadAtomInformation()) {
		return -1;
	}
#my $l = $this->GetAtomInformation("H");
#print "H: $l\n";

	my $out = new JFile;
	if(!$out->Open($path, "w")) {
		return 0;
	}
	$out->print("#VESTA_FORMAT_VERSION 1\n");
	$out->print("\n");

	my $CrystalName = $Crystal->CrystalName();
print "CrystalName: $CrystalName\n";
	$out->print("TITLE\n");
	$out->print("$CrystalName\n");
	$out->print("\n");

	if($DensityFile) {
		$out->print("IMPORT_DENSITY\n");
		$out->print("+1.000000 $DensityFile\n");
		$out->print("\n");
	}

	my $SPG = $Crystal->GetCSpaceGroup();
	my $SPGName = $SPG->SPGName();
	my $iSPG    = $SPG->iSPG();
#my $LatticeSystem = $SPG->LatticeSystem();
print "Space Group: $SPGName ($iSPG)$LF";
#print "Lattice System: $LatticeSystem$LF";
	$out->print("GROUP\n");
	$out->print("$iSPG 1 $SPGName\n");

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

	$out->print("SYMOP\n");
	my $nSymmetryOperation = $SPG->nSymmetryOperation();
	for(my $i = 0 ; $i < $nSymmetryOperation ; $i++) {
		my $symop = $SPG->SymmetryOperation($i+1);
print "   $symop$LF";

		my ($x1,$y1,$z1,$t1,
			$x2,$y2,$z2,$t2,
			$x3,$y3,$z3,$t3) = $SPG->SymmetryOperationByMatrix($i+1);
		$out->printf("%9.6f %9.6f %9.6f  %2d %2d %2d  %2d %2d %2d  %2d %2d %2d\n",
			 $t1, $t2, $t3,
			 $x1, $y1, $z1,
			 $x2, $y2, $z2,
			 $x3, $y3, $z3);
	}
	$out->print(" -1.0 -1.0 -1.0  0 0 0  0 0 0  0 0 0\n");
	$out->print("TRANM\n");
	$out->print(" 0.000000  0.000000  0.000000  1  0  0   0  1  0   0  0  1\n");

	my ($a,$b,$c,$alpha,$beta,$gamma) = $Crystal->LatticeParameters();
print "cell: $a $b $c  $alpha $beta $gamma$LF";
	$out->print("CELLP\n");
	$out->printf(" %9.6f  %9.6f  %9.6f  %9.6f  %9.6f  %9.6f\n",
			$a, $b, $c, $alpha, $beta, $gamma);
	$out->printf("  0.000000   0.000000   0.000000   0.000000   0.000000   0.000000\n");

	$out->print("STRUC\n");
	my @AtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my $nAsymmetricAtomSite = @AtomSiteList;
#print "nAsymmetricAtomSite: $nAsymmetricAtomSite$LF";
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
		my $label     = $AtomSiteList[$i]->Label();
		my $type      = $AtomSiteList[$i]->AtomNameOnly();
		my ($x,$y,$z) = $AtomSiteList[$i]->Position(1);
		my $occupancy = $AtomSiteList[$i]->Occupancy();

print "   $label ($type): ($x, $y, $z) [$occupancy]$LF";
		$out->printf("%3d %2s %7s %7.4f  %9.6f  %9.6f  %9.6f\n",
				$i+1, $type, $label, $occupancy, $x, $y, $z);
		$out->printf("                         0.000000   0.000000   0.000000\n");
	}
	$out->print("  0 0 0 0 0 0 0\n");

	$out->print("THERI 0\n");
	for(my $i = 0 ; $i < $nAsymmetricAtomSite ; $i++) {
		my $label     = $AtomSiteList[$i]->Label();
		$out->printf("%3d %7s %9.6f\n", $i+1, $label, 1.0);
	}
	$out->print("  0 0 0\n");

	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	$out->print("ATOMT\n");
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atomname = $AtomTypeList[$i]->AtomNameOnly();
		my $inf = $this->GetAtomInformation($atomname);
		$out->printf("%3d %2s  %s\n", $i+1, $atomname, $inf);
	}
	$out->print("  0 0 0 0 0 0\n");

	$out->print("MODEL 1 95\n");
	$out->print("RADII 0\n");
	$out->print("BOUND\n");
	$out->printf(" %7f  %7f  %7f  %7f  %7f  %7f\n",
			$LBoundary, $UBoundary, $LBoundary, $UBoundary, $LBoundary, $UBoundary);
	$out->print("  0   0   0   0  0\n");

	$out->print("SBOND\n");
	my $c = 0;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $name1 = $AtomTypeList[$i]->AtomNameOnly();
		next if($name1 =~ /^$Anions$/);
		for(my $j = 0 ; $j < $nAtomType ; $j++) {
			my $name2 = $AtomTypeList[$j]->AtomNameOnly();
#			next if($name2 !~ /^$Anions$/);
			$c++;
			$out->printf("%3d %2s %2s %10.5f  0  %d  1    0.00000  0\n",
				$c, $name1, $name2, $BondLength, $BeyondBoundary);
		}
	}
	$out->print("  0 0 0 0\n");

	my $h = $out->HANDLE();
	print $h <<EOT2;
HBOND 0 2
SPLAN
  0   0   0   0
ATOMS 0
BONDS 0
POLYS 1
ATOMP
 16  16   0  50  2.0   0
BONDP
  1  12  0.423  0.423  1.000 0.7690 0.7690 0.7690
POLYP
  0.500 1  1.000 0.6784 0.6784 0.6784
EOT2

#断面の透明度
	my $t = 0.5;
	$out->print("HKLPP\n");
	$out->printf("  %4.3f 1  1.000  0.0  0.0  0.0\n", $t);

	$out->print("UCOLP\n");
	$out->print("  1  1.000 0.0000 0.0000 0.0000\n");
	$out->print("ISURF\n");
	$out->printf("  0   0  96 %10f %10f %10f\n", $DensityIsoLevel, $DensityMax, $DensityMin);
	$out->print("  1          1          0\n");

	print $h <<EOT3;
BKGRC
   1.0000 1.0000 1.0000
PROJT 1
COMPS 1
DPTHQ 1 -0.5000  3.5000
SCENE
 1.000000  0.000000  0.000000  0.000000
 0.000000  1.000000  0.000000  0.000000
 0.000000  0.000000  1.000000  0.000000
 0.000000  0.000000  0.000000  1.000000
  0.000   0.000
  0.000
  1.000
LIGHT
  1.000   0.000   0.000   0.000
  0.000   1.000   0.000   0.000
  0.000   0.000   1.000   0.000
  0.000   0.000   0.000   1.000
  0.000   0.000  20.000   0.000
  0.000   0.000  -1.000
  0.100   0.100   0.100   1.000
  0.600   0.600   0.600   1.000
ATOMM
  0.250   0.250   0.250   1.000
  0.750   0.750   0.750   1.000
  1.000   1.000   1.000   1.000
128.000
BONDM
  0.250   0.250   0.250   1.000
  0.750   0.750   0.750   1.000
  1.000   1.000   1.000   1.000
128.000
POLYM
  0.250   0.250   0.250   1.000
  0.750   0.750   0.750   1.000
  1.000   1.000   1.000   1.000
128.000
HKLPM
  0.250   0.250   0.250   1.000
  0.750   0.750   0.750   1.000
  1.000   1.000   1.000   1.000
128.000
EOT3

#等値面の透明度,shiness
	$t = 0.5;
	my $s = 0.5;
	my $v = $s * 128;
	$out->print("SURFM\n");
	$out->printf("  1.000   1.000   0.000   %5.3f\n", $t);
	$out->printf("  0.000   1.000   1.000   %5.3f\n", $t);
	$out->print("  0.000   0.000   0.000   1.000\n");
	$out->printf("%7.3f\n", $v);

	print $h <<EOT4;
DELOB
  0

  0

  0

EOT4

	$hash{kVector}    = 1.0 if(!defined $hash{kVector});
	$hash{ArrowColor} = "128:128:128" if(!defined $hash{ArrowColor});
	my ($R, $G, $B)   = Utils::Split(":", $hash{ArrowColor});
	$R /= 255.0;
	$G /= 255.0;
	$B /= 255.0;
	$hash{ArrowRadius}    = 0.2 if(!defined $hash{ArrowRadius});
	$hash{PenetrateAtoms} = 1   if(!defined $hash{PenetrateAtoms});

	my @AtomSiteList = $Crystal->GetCExpandedAtomSiteList();
	my ($vx, $vy, $vz) = $AtomSiteList[0]->Velocity();
	$out->print("ARROW\n");
	if(defined $vx) {
		for(my $i = 0 ; $i < @AtomSiteList ; $i++) {
			($vx, $vy, $vz) = $AtomSiteList[$i]->Velocity();
			next if(!defined $vz);

			my $R2 = $vx*$vx + $vy*$vy + $vz*$vz;
			next if($R2 == 0);

			my $R = sqrt($R2);
			my $vx0 = $vx / $R;
			my $vy0 = $vy / $R;
			my $vz0 = $vz / $R;
			$vx *= $hash{kVector};
			$vy *= $hash{kVector};
			$vz *= $hash{kVector};
			$out->printf("   %d %5.3f %5.3f %5.3f   %5.3f  %5.3f  %5.3f  %5.3f  %5.3f 1\n",
				$i+1, $vx0, $vy0, $vz0, abs($vz), $hash{ArrowRadius}, 
				$R, $G, $B, $hash{PenetrateAtoms});
			$out->print(" $i -1\n");
		}
	}
	$out->print(" 0 0 0 0 0 0 0 0 0\n");

	print $h <<EOT5;
PLN2D
  0   0   0   0
EOT5

	return 1;
}

1;