package POVRay;
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 { shift->{'SampleName'}; }

sub SetTexture
{
	my ($this, $texture) = @_;
	return $this->{'Texture'} = $texture;
}

sub Texture { return shift->{'Texture'}; }

sub PrintHeader
{
	my ($this, $out, $color) = @_;
	$color = "1.0,1.0,1.0" if($color eq '');
	my $RelativeLightIntensity = 0.7;
	if($this->Texture() eq 'Glass') {
		$RelativeLightIntensity = 1.7;
	}

	print $out <<EOT;
//POV-Ray file by Perl Crystal::POVRay.pm

#include "colors.inc"
#include "shapes.inc"
#include "colors.inc"
#include "metals.inc"
#include "glass.inc"
#include "skies.inc"
#include "stars.inc"

//sky_sphere { S_Cloud2 }

background { color rgb <$color> }

camera { 
  orthographic
  location <0.0, 0.0, -50.000000> 
  direction <0.0, 0.0, 1> 
  up  <0.0, 11.989000, 0.0> 
  right <15.985294, 0.0, 0.0> 
  angle 41.147919 
  look_at <0, 0, 0> 
} 

light_source {<-20.0, 20.0, -60.0> color <1.0, 1.0, 1.0>*$RelativeLightIntensity}
light_source {< 20.0, 20.0, -60.0> color <1.0, 1.0, 1.0>*$RelativeLightIntensity}

EOT
	return 1;
}

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

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

	print $out <<EOT;
translation $centerstr
rotation $rotationstr
radius  $cellbarradius
height  $r
EOT
	return 1;
}

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

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

	print $out <<EOT;
sphere { <$rx,$ry,$rz> $r
  material { AtomMaterial }
  pigment  { AtomFil }
}

EOT

	return 1;
}

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

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

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

	print $out <<EOT;
translation $x $y $z
string ["$text"]
family "$FontFamily"
style  "BOLD"
size $size
EOT

	return 1;
}

sub PrintCellRodTermination
{
	my ($this, $out, $x, $y, $z, $scale) = @_;
	$x = sprintf("%8.4f", $x*$scale);
	$y = sprintf("%8.4f", $y*$scale);
	$z = sprintf("%8.4f", $z*$scale);

	print $out <<EOT;
sphere { <$x,$y,$z> CellRodRadius
  material {LineMaterial}
  pigment { CellRodCollor }
}

EOT
}

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

#print "scale: $scale<br>\n";
	$x0 = sprintf("%8.4f", $x0*$scale);
	$y0 = sprintf("%8.4f", $y0*$scale);
	$z0 = sprintf("%8.4f", $z0*$scale);
	$x1 = sprintf("%8.4f", $x1*$scale);
	$y1 = sprintf("%8.4f", $y1*$scale);
	$z1 = sprintf("%8.4f", $z1*$scale);

	print $out <<EOT;
cylinder { <$x0,$y0,$z0> <$x1,$y1,$z1>, CellRodRadius open
  material { LineMaterial }
  pigment  { CellRodCollor }
}
EOT
	return 1;
}

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

	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->PrintCellRod($out, 0, 0, 0, $a11, $a12, $a13, $scale);
	$this->PrintCellRod($out, 0+$a21, 0+$a22, 0+$a23, 
		     $a11+$a21, $a12+$a22, $a13+$a23, $scale);
	$this->PrintCellRod($out, 0+$a31, 0+$a32, 0+$a33, 
		     $a11+$a31, $a12+$a32, $a13+$a33, $scale);
	$this->PrintCellRod($out, 0+$a21+$a31, 0+$a22+$a32, 0+$a23+$a33, 
		     $a11+$a21+$a31, $a12+$a22+$a32, $a13+$a23+$a33, $scale);
	$this->PrintCellRod($out, 0, 0, 0, $a21, $a22, $a23, $scale);
	$this->PrintCellRod($out, 0+$a11, 0+$a12, 0+$a13, 
		     $a21+$a11, $a22+$a12, $a23+$a13, $scale);
	$this->PrintCellRod($out, 0+$a31, 0+$a32, 0+$a33, 
		     $a21+$a31, $a22+$a32, $a23+$a33, $scale);
	$this->PrintCellRod($out, 0+$a11+$a31, 0+$a12+$a32, 0+$a13+$a33, 
		     $a21+$a11+$a31, $a22+$a12+$a32, $a23+$a13+$a33, $scale);
	$this->PrintCellRod($out, 0, 0, 0, $a31, $a32, $a33, $scale);
	$this->PrintCellRod($out, 0+$a11, 0+$a12, 0+$a13, 
		     $a31+$a11, $a32+$a12, $a33+$a13, $scale);
	$this->PrintCellRod($out, 0+$a21, 0+$a22, 0+$a23, 
		     $a31+$a21, $a32+$a22, $a33+$a23, $scale);
	$this->PrintCellRod($out, 0+$a11+$a21, 0+$a12+$a22, 0+$a13+$a23, 
		     $a31+$a11+$a21, $a32+$a12+$a22, $a33+$a13+$a23, $scale);

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

	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;
sphere { <$rx,$ry,$rz> AtomRad_${label}
  material { AtomMaterial }
  pigment  { AtomFil_${label} }
}

EOT

	return 1;
}

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

#print "scale: $scale<br>\n";
	$x0 = sprintf("%8.4f", $x0*$scale);
	$y0 = sprintf("%8.4f", $y0*$scale);
	$z0 = sprintf("%8.4f", $z0*$scale);
	$x1 = sprintf("%8.4f", $x1*$scale);
	$y1 = sprintf("%8.4f", $y1*$scale);
	$z1 = sprintf("%8.4f", $z1*$scale);

	print $out <<EOT;
cylinder { <$x0,$y0,$z0> <$x1,$y1,$z1>, BondRad open
  material { BondMaterial }
  pigment  { BondFil }
}

EOT
}

sub SaveFile
{
	my ($this, $Crystal, $filename, $nx, $ny, $nz, 
			$MaxBondLength, $nMaxBonds, $Scale, $rScale, $texture) = @_;
	my $LF = "<br>\n";
	$this->SetTexture($texture);

	my $Texture = 'T_Glass1';
	my $Ior     = 1.2;

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

# Header
	$this->PrintHeader($out, "1.0,1.0,1.0");

	if($this->Texture() eq 'Glass') {

		print $out <<EOTPRIMITIVES1;
//==============================================
//  Declarations for primitives
//==============================================

#declare AtomMaterial = material {
  texture  { $Texture }
  interior { ior $Ior }
}

#declare BondMaterial = material {
  texture  { $Texture }
  interior { ior $Ior }
}

#declare PolyMaterial = material {
  texture  { $Texture }
  interior { ior $Ior }
}

#declare LineMaterial = material {
  texture  { $Texture }
  interior { ior $Ior }
}

EOTPRIMITIVES1

	}
	else {
		print $out <<EOTPRIMITIVES2;
//==============================================
//  Declarations for primitives
//==============================================

#declare AtomMaterial = material {
  texture{
    finish {
      ambient  0.100000
      diffuse  1.000000
      specular 1.000000
    }
  }
}

#declare BondMaterial = material {
  texture{
    finish {
      ambient  0.100000
      diffuse  1.000000
      specular 1.000000
    }
  }
}

#declare PolyMaterial = material {
  texture{
    finish {
      ambient  0.100000
      diffuse  1.000000
      specular 0.000000
    }
  }
}

#declare LineMaterial = material {
  texture{
    finish {
      ambient  0.100000
      diffuse  1.000000
      specular 0.000000
    }
  }
}

EOTPRIMITIVES2
	
	}

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

	print $out <<EOTPROTOATOM;
//===================================================
//   Prototype declaretion for atoms
//===================================================
EOTPROTOATOM

	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 = sprintf("%7.4f", $AtomRadius[$iAtom-1] * $Scale * $rScale);

		if($this->SetTexture($texture) eq 'Glass') {
			print $out <<EOTATOM;
#declare AtomFil_${label} = color rgbf <$color, 0.9>;
#declare AtomRim_${label} = color rgbf <$color, 0.9>;
#declare AtomRad_${label} = $radius;

EOTATOM
		}
		else {
			print $out <<EOTATOM;
#declare AtomFil_${label} = color rgb <$color>;
#declare AtomRim_${label} = color rgb <$color>;
#declare AtomRad_${label} = $radius;

EOTATOM
		}
	}

	my $BondFilColor  = "0.9,0.9,0.9";
	my $BondRimColor  = "0.0,0.0,0.0";
	my $ScaledBondRadius = sprintf("%8.4f", $BondRadius * $Scale);
	print $out <<EOTPROTOBOND;
//===================================================
//   Prototype declaretion for chemical bonds
//===================================================
#declare BondFil = color rgb <$BondFilColor>;
#declare BondRim = color rgb <$BondRimColor>;
#declare BondRad = $ScaledBondRadius;

EOTPROTOBOND

	print $out <<EOTPROTO;
//===================================================
//   Prototype declaretion for crystal components
//===================================================

//polyhedron colors, sizes
#declare polyrad    = 0.050000;
#declare polysphere = 0.200000;

//crystal edges, faces
#declare edgecolfrt = color rgb <1.000000, 1.000000, 0.000000>;
#declare edgecolbck = color rgb <0.000000, 0.000000, 1.000000>;
#declare edgerad = 0.050000;
#declare xtlfacecol = color rgbf <0.000000, 1.000000, 1.000000, 0.500000>;

//unit-cell edges, faces
#declare CellRodCollor = color rgb <0.000000, 0.000000, 0.000000>;
#declare CellRodRadius = 0.050000;
#declare UnitCellFaceColor = color rgbf <1.000000, 1.000000, 0.000000, 0.500000>;

//crystal axes
#declare axescol = color rgb <0.000000, 0.000000, 0.000000>;
#declare axesrad = 0.050000;

//cavity colors
#declare cavfacecol = color rgb <1.000000, 1.000000, 0.000000>;
#declare cavendcol = color rgbf <1.000000, 0.000000, 0.000000, 0.0>;

EOTPROTO

# Start drawing
	print $out "union {\n";

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

EOTLATT
	$this->PrintLattice($out, $Crystal, $CellBarRadius, $Scale);
	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 $Color = $AtomColor[$iAtom-1];
		$nAtomTypes{$atomname}++;
		my $label = $SiteLabel;

		my $occscale = (0.3 + $occupancy) / 1.3;
		$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 <<EOTBONDS;
//===================================================
//   Draw Chemical bonds
//===================================================
EOTBONDS

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

        no_shadow
        scale 1.5
}

EOTEND

	close($out);

	return ($filename);
}

1;