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

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

use strict;
use File::Basename;

use Sci::Science;

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

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 PrintVRMLHeader
{
	my ($this, $out, $color, $Ver) = @_;
	return $this->PrintVRMLHeaderV1($out, $color) if($Ver == 1.0);
	return $this->PrintVRMLHeaderV2($out, $color);
}

sub PrintVRMLHeaderV1
{
	my ($this, $out, $color) = @_;
	$color = "1 1 1" if($color eq '');

	print $out "#VRML V1.0 ascii\n";
	print $out "\n";
	print $out "DEF BackgroundColor Info {\n";
	print $out "	string \"$color\"\n";
	print $out "}\n";

	return 1;
}

sub PrintVRMLHeaderV2
{
	my ($this, $out, $color) = @_;
	$color = "1 1 1" if($color eq '');

	print $out <<EOT;
#VRML V2.0 utf8

Background {
  groundAngle [ 0.0 ]
  groundColor [ 1.0 1.0 1.0 ]
  skyAngle    [ 0.0 ]
  skyColor    [ 1.0 1.0 1.0 ]
}

EOT
	return 1;
}

sub PrintCylinder
{
	my ($this, $out, $x0, $y0, $z0, $x1, $y1, $z1, 
		$cellbarradius, $scale, $Ver) = @_;

#print "scale: $scale<br>\n";
	$x0 *= $scale;
	$y0 *= $scale;
	$z0 *= $scale;
	$x1 *= $scale;
	$y1 *= $scale;
	$z1 *= $scale;
	$cellbarradius *= $scale;
	my $centerx = ($x0 + $x1) / 2.0;
	my $centery = ($y0 + $y1) / 2.0;
	my $centerz = ($z0 + $z1) / 2.0;
	my $centerxm = -$centerx;
	my $centerym = -$centery;
	my $centerzm = -$centerz;
	my $dx = $x1 - $x0;
	my $dy = $y1 - $y0;
	my $dz = $z1 - $z0;
	my $r = $dx*$dx + $dy*$dy + $dz*$dz;
	$r = sqrt($r);
#	my $cosQx = $dz / sqrt($dy*$dy + $dz*$dz);
#	my $Qx = Sci::acos($cosQx);
#	my $cosQz = $dy / sqrt($dx*$dx + $dy*$dy);
#	my $Qz = Sci::acos($cosQz);
	my $cosQ = $dy / $r;
	my $Q = Sci::acos($cosQ);

	my $dxm = -$dx;
	my $dzm = -$dz;
	my $Qm = -$Q;

	my $centerstr  = sprintf("%8.4f %8.4f %8.4f", $centerx,  $centery,  $centerz);
	my $centermstr = sprintf("%8.4f %8.4f %8.4f", $centerxm, $centerym, $centerzm);
	my $rotationstr  = sprintf("%8.4f %8.4f %8.4f   %8.4f", 
					$dz, 0, $dxm, $Q);
	my $rotationmstr = sprintf("%8.4f %8.4f %8.4f   %8.4f", 
					$dz, 0, $dxm, $Qm);
	$cellbarradius = sprintf("%8.4f", $cellbarradius);
	$r = sprintf("%8.4f", $r);

	if($Ver == 1.0) {
		print $out <<EOT1;
Separator {
	Transform {
		translation $centerstr
		rotation $rotationstr
	}
	Cylinder {
		parts   ALL
		radius  $cellbarradius
		height  $r
	}
	Transform {
		translation $centermstr
		rotation $rotationmstr
	}
}
EOT1
	}
	else {
		print $out <<EOT2;
Transform {
	translation $centerstr
	rotation $rotationstr
	children [
		Shape {
			appearance Appearance {
				material Material {}
			}
			geometry Cylinder {
				radius  $cellbarradius
				height  $r
			}
		}
	]
}
EOT2
	}
	return 1;
}

sub PrintSphere
{
	my ($this, $out, $rx, $ry, $rz, $r, $color, $scale, $rscale, $Ver) = @_;

	$rscale = 1.0 if($rscale <= 0.0);
	$rx *= $scale;
	$ry *= $scale;
	$rz *= $scale;
	$r  *= $scale*$rscale;
	my $rxm = -$rx;
	my $rym = -$ry;
	my $rzm = -$rz;
	
	$rx  = sprintf("%8.4f", $rx);
	$ry  = sprintf("%8.4f", $ry);
	$rz  = sprintf("%8.4f", $rz);
	$rxm = sprintf("%8.4f", $rxm);
	$rym = sprintf("%8.4f", $rym);
	$rzm = sprintf("%8.4f", $rzm);
	$r   = sprintf("%8.4f", $r);

	if($Ver == 1.0) {
		print $out <<EOT1;
Separator {
	Material {
		diffuseColor $color
	}
	Transform {
		translation $rx $ry $rz
	}
	Sphere {
		radius $r
	}
	Transform {
		translation $rxm $rym $rzm
	}
}
EOT1
	}
	else {
		print $out <<EOT2;
Transform {
	translation $rx $ry $rz
	children [
		Shape {
			appearance Appearance {
				material Material {
					diffuseColor $color
				}
			}
			geometry Sphere {
				radius $r
			}
		}
	]
}
EOT2
	}

	return 1;
}

sub PrintText
{
	my ($this, $out, $x, $y, $z, $text, $size, $scale, $Ver) = @_;
	my $FontFamily = "SERIF"; #"TYPEWRITER"

	$x *= $scale;
	$y *= $scale;
	$z *= $scale;
	my $xm = -$x;
	my $ym = -$y;
	my $zm = -$z;
	$size *= $scale;

	if($Ver == 1.0) {
	}
	else {
		$y -= $size * 1.0;
	}
	
	$x  = sprintf("%8.4f", $x);
	$y  = sprintf("%8.4f", $y);
	$z  = sprintf("%8.4f", $z);
	$xm = sprintf("%8.4f", $xm);
	$ym = sprintf("%8.4f", $ym);
	$zm = sprintf("%8.4f", $zm);
	$size = sprintf("%9.4f", $size);

	if($Ver == 1.0) {
		print $out <<EOT1;
Separator {
	Transform {
		translation $x $y $z
	}
	FontStyle {
		size $size
	}
	AsciiText {
		string "$text"
	}
	Transform {
		translation $xm $ym $zm
	}
}
EOT1
	}
	else {
		print $out <<EOT2;
Transform {
	translation $x $y $z
	children [
		Shape {
			appearance Appearance {
				material Material {}
			}
			geometry Text {
				string ["$text"]
				fontStyle FontStyle {
					family "$FontFamily"
					style  "BOLD"
					size $size
				}
			}
		}
	]
}
EOT2
	}

	return 1;
}

sub PrintLattice
{
	my ($this, $out, $Crystal, $cellbarradius, $scale, $Ver) = @_;

	my ($a11,$a13,$a12,
	    $a31,$a33,$a32,
	    $a21,$a23,$a22) = $Crystal->LatticeVectors();
#print "Latt: $a11,$a12,$a13<br>\n";
#print "Latt: $a21,$a22,$a23<br>\n";
#print "Latt: $a31,$a32,$a33<br>\n";

	$this->PrintCylinder($out, 0, 0, 0, $a11, $a12, $a13, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a21, 0+$a22, 0+$a23, 
			     $a11+$a21, $a12+$a22, $a13+$a23, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a31, 0+$a32, 0+$a33, 
			     $a11+$a31, $a12+$a32, $a13+$a33, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a21+$a31, 0+$a22+$a32, 0+$a23+$a33, 
			     $a11+$a21+$a31, $a12+$a22+$a32, $a13+$a23+$a33, 
			     				      $cellbarradius, $scale, $Ver);

	$this->PrintCylinder($out, 0, 0, 0, $a21, $a22, $a23, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a11, 0+$a12, 0+$a13, 
			     $a21+$a11, $a22+$a12, $a23+$a13, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a31, 0+$a32, 0+$a33, 
			     $a21+$a31, $a22+$a32, $a23+$a33, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a11+$a31, 0+$a12+$a32, 0+$a13+$a33, 
			     $a21+$a11+$a31, $a22+$a12+$a32, $a23+$a13+$a33, 
			     				      $cellbarradius, $scale, $Ver);

	$this->PrintCylinder($out, 0, 0, 0, $a31, $a32, $a33, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a11, 0+$a12, 0+$a13, 
			     $a31+$a11, $a32+$a12, $a33+$a13, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a21, 0+$a22, 0+$a23, 
			     $a31+$a21, $a32+$a22, $a33+$a23, $cellbarradius, $scale, $Ver);
	$this->PrintCylinder($out, 0+$a11+$a21, 0+$a12+$a22, 0+$a13+$a23, 
			     $a31+$a11+$a21, $a32+$a12+$a22, $a33+$a13+$a23, 
			     				      $cellbarradius, $scale, $Ver);

	return 1;
}

sub SaveFile
{
	my ($this, $Crystal, $filename, $nx, $ny, $nz, 
			$MaxBondLength, $nMaxBonds, $Scale, $rScale, $VRMLVersion) = @_;
	my $Ver = 2.0;
	$Ver = 1.0 if($VRMLVersion =~ /1/);
	my $LF = "<br>\n";

#ファイル作製開始
	my $out;
	unless(open($out,">$filename")) {
		print "Can not write to $filename.$LF$LF";
		return;
	}
	if($Ver == 1.0) {
		$this->PrintVRMLHeaderV1($out, "0.3 0.3 0.3");
	}
	else {
		$this->PrintVRMLHeaderV2($out, "0.3 0.3 0.3");
	}
	
	my $SampleName = $this->SampleName();
	my $SPGName = $Crystal->SPGNameByOutputMode();
	my $iSPG    = $Crystal->iSPGByOutputMode();
	my @AtomTypeList = $Crystal->GetCAtomTypeList();
	my $nAtomType = @AtomTypeList;
	my @ExpandedAtomSiteList    = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nExpandedAtomSite = @ExpandedAtomSiteList;
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(1);

	my $Max = $a*$nx;
	$Max = $b*$ny if($Max < $b*$ny);
	$Max = $c*$nz if($Max < $c*$nz);
	$Max = 1.0 if($Max <= 0.0);
	$Scale = $Scale / $Max * 10.0;

	my $CellBarRadius = 0.05;
	my $BondRadius = 0.1;

	my $centerx = -(-$nx+0.98 + 1.02) * $a / 2.0 * $Scale;
	my $centery = -(-$ny+0.98 + 1.02) * $b / 2.0 * $Scale;
	my $centerz = -(-$nz+0.98 + 1.02) * $c / 2.0 * $Scale;
#print "aa ($nx): $centerx, ", $a*$Scale, " (", -$nx+0.98+1.02, ")$LF";
	$centerx = sprintf("%8.4f", $centerx);
	$centery = sprintf("%8.4f", $centery);
	$centerz = sprintf("%8.4f", $centerz);

	if($Ver == 1.0) {
		print $out <<EOT11;
Separator {
	Transform {
		translation $centerx $centerz $centery
	}

EOT11
	}
	else {
	}

	$this->PrintLattice($out, $Crystal, $CellBarRadius, $Scale, $Ver);
	print $out "\n";

	my $ShowCrystal = new Crystal;
	$ShowCrystal = $Crystal;
	$ShowCrystal->SetnAsymmetricAtomSite(0);
	$ShowCrystal->ExpandCoordinates();

	my $ColorOffset = 0.3;
	my $inx = int($nx-0.01) + 1;
	my $iny = int($ny-0.01) + 1;
	my $inz = int($nz-0.01) + 1;
	my %nAtomTypes;
	my @AtomColor;
	for(my $ix = -$inx ; $ix <= $inx ; $ix++) {
	for(my $iy = -$iny ; $iy <= $iny ; $iy++) {
	for(my $iz = -$inz ; $iz <= $inz ; $iz++) {

	for(my $i = 0 ; $i < $nExpandedAtomSite ; $i++) {
		my $AtomSite  = $ExpandedAtomSiteList[$i];
		my $atomname  = $AtomSite->AtomNameOnly();
		my $charge    = $AtomSite->Charge();
		my ($x,$y,$z) = $AtomSite->Position(1);
#print "$atomname: ($x,$y,$z)+($ix,$iy,$iz)$LF";
		$x += $ix;
		$y += $iy;
		$z += $iz;
		next if($x < -$nx+0.98 or 1.02 < $x);
		next if($y < -$ny+0.98 or 1.02 < $y);
		next if($z < -$nz+0.98 or 1.02 < $z);
#print "$atomname: ($x,$y,$z)$LF";

		my $occupancy = $AtomSite->Occupancy();
		my ($xc, $yc, $zc) = $Crystal->FractionalToCartesian($x, $y, $z);

		my $iAtom = $Crystal->FindiAtomType($atomname);
		my $AtomType  = $AtomTypeList[$iAtom-1];
		my $Radius    = $AtomType->IonRadius(); #AtomicRadius();
#print "r($atomname)=$Radius<br>\n";
#		$Radius = 1.4 if($Radius <= 0 and $atomname eq 'O');
		$Radius = 0.6 if($Radius <= 0);
		my $Color = $AtomColor[$iAtom-1];
		if($Color eq '') {
			my $red   = $ColorOffset + rand(1.0-$ColorOffset);
			my $green = $ColorOffset + rand(1.0-$ColorOffset);
			my $blue  = $ColorOffset + rand(1.0-$ColorOffset);
			$red   /= 3 if($iAtom - 1 % 3 == 0);
			$green /= 3 if($iAtom - 1 % 3 == 1);
			$blue  /= 3 if($iAtom - 1 % 3 == 2);
			$red   = sprintf("%7.4f", $red);
			$green = sprintf("%7.4f", $green);
			$blue  = sprintf("%7.4f", $blue);
			$Color = $AtomColor[$iAtom-1] = "$red $green $blue";
		}
		$nAtomTypes{$atomname}++;
		my $label = $atomname . $nAtomTypes{$atomname};

		my $occscale = (0.3 + $occupancy) / 1.3;
		$this->PrintSphere($out, $xc, $zc, $yc, $Radius*$occscale, 
					$Color, $Scale, $rScale, $Ver);
		$ShowCrystal->AddAtomSite($label, $atomname, $x, $y, $z, $occupancy);
#print "$label($atomname): ($x,$y,$z)$LF";
	}
	
	}
	}
	}
	print $out "\n";

	my @alist = $ShowCrystal->GetCAsymmetricAtomSiteList();
	my $count = 0;
	for(my $i = 0 ; $i < @alist ; $i++) {
		my $atom0 = $alist[$i];
		my $name0 = $atom0->AtomNameOnly();
		my ($x0,$y0,$z0) = $atom0->Position(0);
		my ($xc0, $yc0, $zc0) = $ShowCrystal->FractionalToCartesian($x0, $y0, $z0);

		for(my $j = $i+1 ; $j < @alist ; $j++) {
			my $atom1 = $alist[$j];
			my $name1 = $atom1->AtomNameOnly();
			my ($x1,$y1,$z1) = $atom1->Position(0);
			
			my $dis = $ShowCrystal->GetInterAtomicDistance(
						$x0, $y0, $z0, $x1, $y1, $z1);
#print "($i,$j) ($name0,$name1) ($x0,$y0,$z0)-($x1,$y1,$z1) $dis$LF";
			next if($dis < 0.1 or $dis > $MaxBondLength);
#print "...Show ($MaxBondLength)$LF";
			my ($xc1, $yc1, $zc1) = $ShowCrystal->FractionalToCartesian($x1, $y1, $z1);
			$this->PrintCylinder($out, $xc0, $zc0, $yc0,
				$xc1, $zc1, $yc1, $BondRadius, $Scale, $Ver);

			$count++;
			if($count > $nMaxBonds) {
				print "Number of bonds reaches maximum ($nMaxBonds).$LF";
				last;
			}
		}
		last if($count > $nMaxBonds);
	}

	print $out "\n";
	print $out "}\n";
	print $out "\n";
	close($out);


#原子名のVRMLファイルを作成
	my $NameFile = $filename;
	$NameFile =~ s/\..*?$//;
	$NameFile = "$NameFile-name.wrl";
	my $out;
	unless(open($out,">$NameFile")) {
		print "Can not write to $NameFile.$LF$LF";
		return;
	}

	if($Ver == 1.0) {
		$this->PrintVRMLHeaderV1($out, "0.3 0.3 0.3");
	}
	else {
		$this->PrintVRMLHeaderV2($out, "0.3 0.3 0.3");
	}

	my $TotalRad = 0.0;
	my $MaxRad = 0.0;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom  = $AtomTypeList[$i];
		my $atomname = $atom->AtomNameOnly();
		my $Radius    = $atom->IonRadius(); #AtomicRadius();
#		$Radius = 1.4 if($Radius <= 0 and $atomname eq 'O');
		$Radius = 0.6 if($Radius <= 0);
		$TotalRad += $Radius;
		$MaxRad = $Radius if($Radius > $MaxRad);
	}
	my $ShiftY = $MaxRad * $nAtomType / 2.0 * 0.7 ;
	$ShiftY = sprintf("%8.4f", $ShiftY);

	if($Ver == 1.0) {
		print $out <<EOT2;
Separator {
	Transform {
		translation 0 $ShiftY 0
	}

EOT2
	}

	my $xs = 0.0;
	my $ys = 0.0;
	my $zs = 0.0;
	for(my $i = 0 ; $i < $nAtomType ; $i++) {
		my $atom  = $AtomTypeList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $Radius    = $atom->IonRadius(); #AtomicRadius();
print "r($atomname)=$Radius<br>\n";
#		$Radius = 1.4 if($Radius <= 0 and $atomname eq 'O');
		$Radius = 0.6 if($Radius <= 0);
		my $Color = $AtomColor[$i];

		$this->PrintSphere($out, $xs, $ys, $zs, $Radius, $Color, $Scale, 
						$rScale, $Ver);
		$this->PrintText(  $out, $xs+$MaxRad*1.5, $ys+2.0*$Scale, $zs, 
						$atomname, 2.0, $Scale, $Ver);
		$ys -= $MaxRad * 2 * 0.7;
	}

	print $out "}\n";
	print $out "\n";
	close($out);

	return ($filename, $NameFile);
}

1;