package VRML;
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 PrintX3DVRMLEncodingHeader
{
	my ($this, $out, $color) = @_;
	$color = "1 1 1" if($color eq '');

	print $out <<EOT;
#X3D V3.0 utf8
PROFILE Immersive

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 PrintCellRod
{
	my ($this, $out, $x0, $y0, $z0, $x1, $y1, $z1, 
		$cellbarradius, $scale, $Ver) = @_;
	return $this->PrintCylinder($out, $out, $x0, $y0, $z0, $x1, $y1, $z1, 
			$cellbarradius, $scale, $Ver) if($Ver == 1.0);

#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 $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);


	print $out <<EOT;
CellRod {
	CellRodTranslation $centerstr
	CellRodRotation    $rotationstr
	CellRodLength      $r
}

EOT
	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";

	if($Ver == 1.0) {
		$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);
	}
	else {
		$this->PrintCellRod($out, 0, 0, 0, $a11, $a12, $a13, $cellbarradius, $scale, $Ver);
		$this->PrintCellRod($out, 0+$a21, 0+$a22, 0+$a23, 
			     $a11+$a21, $a12+$a22, $a13+$a23, $cellbarradius, $scale, $Ver);
		$this->PrintCellRod($out, 0+$a31, 0+$a32, 0+$a33, 
			     $a11+$a31, $a12+$a32, $a13+$a33, $cellbarradius, $scale, $Ver);
		$this->PrintCellRod($out, 0+$a21+$a31, 0+$a22+$a32, 0+$a23+$a33, 
			     $a11+$a21+$a31, $a12+$a22+$a32, $a13+$a23+$a33, 
			     				      $cellbarradius, $scale, $Ver);

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

		$this->PrintCellRod($out, 0, 0, 0, $a31, $a32, $a33, $cellbarradius, $scale, $Ver);
		$this->PrintCellRod($out, 0+$a11, 0+$a12, 0+$a13, 
			     $a31+$a11, $a32+$a12, $a33+$a13, $cellbarradius, $scale, $Ver);
		$this->PrintCellRod($out, 0+$a21, 0+$a22, 0+$a23, 
			     $a31+$a21, $a32+$a22, $a33+$a23, $cellbarradius, $scale, $Ver);
		$this->PrintCellRod($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 PrintAtom
{
	my ($this, $out, $label, $rx, $ry, $rz, $scale) = @_;

	$rx  = sprintf("%8.4f", $rx*$scale);
	$ry  = sprintf("%8.4f", $ry*$scale);
	$rz  = sprintf("%8.4f", $rz*$scale);

	print $out <<EOT;
${label}Atom {
	${label}AtomTranslation $rx $ry $rz

}

EOT
	return 1;
}

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

#print "scale: $scale<br>\n";
	$x0 *= $scale;
	$y0 *= $scale;
	$z0 *= $scale;
	$x1 *= $scale;
	$y1 *= $scale;
	$z1 *= $scale;
	my $centerx = ($x0 + $x1) / 2.0;
	my $centery = ($y0 + $y1) / 2.0;
	my $centerz = ($z0 + $z1) / 2.0;
	my $dx = $x1 - $x0;
	my $dy = $y1 - $y0;
	my $dz = $z1 - $z0;
	my $dxm = -$dx;

	my $r = $dx*$dx + $dy*$dy + $dz*$dz;
	$r = sprintf("%7.4f", sqrt($r));
	my $cosQ = $dy / $r;
	my $Q = Sci::acos($cosQ);

	my $centerstr  = sprintf("%8.4f %8.4f %8.4f", $centerx,  $centery,  $centerz);
	my $rotationstr  = sprintf("%8.4f %8.4f %8.4f   %8.4f", $dz, 0, $dxm, $Q);

	print $out <<EOT;
${plabel}Bond {
	${plabel}BondTranslation $centerstr
	${plabel}BondRotation    $rotationstr
	${plabel}BondLength      $r
}

EOT
}

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

#ファイル作製開始
	my $out;
	unless(open($out,">$filename")) {
		print "Can not write to $filename.$LF$LF";
		return;
	}

# Header
	if($Ver == 1.0) {
		$this->PrintVRMLHeaderV1($out, "0.3 0.3 0.3");
	}
	elsif($Ver == 2.0) {
		$this->PrintVRMLHeaderV2($out, "0.3 0.3 0.3");
	}
	else {
		$this->PrintX3DVRMLEncodingHeader($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 @AsymmetricAtomSiteList = $Crystal->GetCAsymmetricAtomSiteList();
	my @ExpandedAtomSiteList    = $Crystal->GetCExpandedAtomSiteListByOutputMode();
	my $nExpandedAtomSite = @ExpandedAtomSiteList;
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(1);

	my $ColorOffset = 0.3;
	my @AtomColor;
	my @AtomRadius;

	for(my $i = 0 ; $i < @AtomTypeList ; $i++) {
		my $AtomType = $AtomTypeList[$i];
		my $atomname = $AtomType->AtomNameOnly();
		$AtomRadius[$i]    = $AtomType->IonRadius(); #AtomicRadius();
#print "r($atomname)=$AtomRadius[$i]<br>\n";
#		$AtomRadius[$i] = 1.4 if($AtomRadius[$i] <= 0 and $atomname eq 'O');
		$AtomRadius[$i] = 0.6 if($AtomRadius[$i] <= 0);
		if($AtomColor[$i] 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($i % 3 == 0);
			$green /= 3 if($i % 3 == 1);
			$blue  /= 3 if($i % 3 == 2);
			$red   = sprintf("%7.4f", $red);
			$green = sprintf("%7.4f", $green);
			$blue  = sprintf("%7.4f", $blue);
			$AtomColor[$i] = "$red $green $blue";
		}
	}

	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);

	my $ScaledCellBarRadius = sprintf("%8.4f", $CellBarRadius * $Scale);

	if($Ver == 1.0) {
	}
	else {
		print $out <<EOTPROTO;
#===================================================
#   Prototype declarations for primitives
#===================================================
PROTO MyCylinder [
	field SFColor    MyColor        0.5  0.5  0.5
	field SFFloat    MyRadius       0.0424
	field SFFloat    MyLength       1.0
]
{
      Shape {
        appearance Appearance {
          material Material {
            diffuseColor IS MyColor
          }
        }
        geometry Cylinder {
          radius IS MyRadius
          height IS MyLength
        }
      }
}

PROTO MySphere [
	field SFColor    MyColor        0.5  0.5  0.5
	field SFFloat    MyRadius       0.0424
]
{
      Shape {
        appearance Appearance {
          material Material {
            diffuseColor IS MyColor
          }
        }
	geometry Sphere {
          radius IS MyRadius
        }
      }
}

#===================================================
#   Prototype declarations for crystal components
#===================================================
PROTO CellRod [
	field SFVec3f    CellRodTranslation 0.0  0.0  0.0
	field SFRotation CellRodRotation    1.0  0.0  0.0    0.0
	field SFColor    CellRodColor       0.9  0.9  0.9
	field SFFloat    CellRodRadius      $ScaledCellBarRadius
	field SFFloat    CellRodLength      1.0
]
{
  Transform {
    translation IS CellRodTranslation
    rotation    IS CellRodRotation
    children [
      MyCylinder {
          MyColor  IS CellRodColor
          MyRadius IS CellRodRadius
          MyLength IS CellRodLength
      }
    ]
  }
}

EOTPROTO

		for(my $i = 0 ; $i < @AsymmetricAtomSiteList ; $i++) {
			my $site = $AsymmetricAtomSiteList[$i];
			my $atomname = $site->AtomNameOnly();
			my $label  = $site->Label();
			my $iAtom  = $Crystal->FindiAtomType($atomname);
			my $color  = $AtomColor[$iAtom-1];
			my $radius = $AtomRadius[$iAtom-1] * $Scale;

			print $out <<EOTATOM
PROTO ${label}Atom [
	field SFVec3f ${label}AtomTranslation 0.0  0.0  0.0
	field SFColor ${label}AtomColor       $color
	field SFFloat ${label}AtomRadius      $radius
]
{
  Transform {
    translation IS ${label}AtomTranslation
    children [
	MySphere {
	    MyColor  IS ${label}AtomColor
	    MyRadius IS ${label}AtomRadius
	}
    ]
  }
}

EOTATOM
		}

		my $BondColor  = "0.9 0.9 0.9";
		my $ScaledBondRadius = sprintf("%8.4f", $BondRadius * $Scale);
		for(my $i = 0 ; $i < @AsymmetricAtomSiteList ; $i++) {
			my $site1  = $AsymmetricAtomSiteList[$i];
			my $label1 = $site1->Label();
			for(my $j = $i ; $j < @AsymmetricAtomSiteList ; $j++) {
				my $site2  = $AsymmetricAtomSiteList[$j];
				my $label2 = $site2->Label();
				my $plabel;
				if($label1 gt $label2) {
					$plabel = "$label2$label1";
				}
				else {
					$plabel = "$label1$label2";
				}
				print $out <<PROTOBOND;
PROTO ${plabel}Bond [
	field SFVec3f    ${plabel}BondTranslation 0.0  0.0  0.0
	field SFRotation ${plabel}BondRotation    1.0  0.0  0.0    0.0
	field SFColor    ${plabel}BondColor       $BondColor
	field SFFloat    ${plabel}BondRadius      $ScaledBondRadius
	field SFFloat    ${plabel}BondLength      1.0
]
{
  Transform {
    translation IS ${plabel}BondTranslation
    rotation    IS ${plabel}BondRotation
    children [
      MyCylinder {
          MyColor  IS ${plabel}BondColor
          MyRadius IS ${plabel}BondRadius
          MyLength IS ${plabel}BondLength
      }
    ]
  }
}

PROTOBOND
			}
		}
	}

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

EOT11
	}
	else {
	}

	print $out <<EOTLATT;
#===================================================
#   Unit cell
#===================================================

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

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


	print $out <<EOTATOMS;
#===================================================
#   Draw Atoms
#===================================================
EOTATOMS

	my $inx = int($nx-0.01) + 1;
	my $iny = int($ny-0.01) + 1;
	my $inz = int($nz-0.01) + 1;
	my %nAtomTypes;
	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 $SiteLabel = $AtomSite->Label();
#print "<b>Label: $SiteLabel</b><br>\n";
		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   = $AtomRadius[$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 $label = $SiteLabel;

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

	print $out <<EOTATOMS;
#===================================================
#   Draw Chemical bonds
#===================================================
EOTATOMS

	my @alist = $ShowCrystal->GetCAsymmetricAtomSiteList();
	my $count = 0;
	for(my $i = 0 ; $i < @alist ; $i++) {
		my $atom0 = $alist[$i];
		my $SiteLabel0 = $atom0->Label();
#print "<b>SiteLabel0: $SiteLabel0</b> - \n";
		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 $SiteLabel1 = $atom1->Label();
#print "<b>SiteLabel1: $SiteLabel1</b><br>\n";
			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);
			if($Ver == 1.0) {
				$this->PrintCylinder($out, $xc0, $zc0, $yc0,
					$xc1, $zc1, $yc1, $BondRadius, $Scale, $Ver);
			}
			else {
				my $plabel;
				if($SiteLabel0 gt $SiteLabel1) {
					$plabel = "$SiteLabel1$SiteLabel0";
				}
				else {
					$plabel = "$SiteLabel0$SiteLabel1";
				}
				$this->PrintChemicalBond($out, $plabel, $xc0, $zc0, $yc0,
					$xc1, $zc1, $yc1, $Scale);
			}
			$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/\..*?$//;
	if($Ver > 2.0) {
		$NameFile = "$NameFile-name.x3dv";
	}
	else {
		$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");
	}
	elsif($Ver == 2.0) {
		$this->PrintVRMLHeaderV2($out, "0.3 0.3 0.3");
	}
	else {
		$this->PrintX3DVRMLEncodingHeader($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;