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

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

use strict;
use File::Basename;

use Sci::Science;
use Sci::EnergyBand;
use GraphData;
use MyTk::GraphFrameArray;

use Crystal::Crystal;
use Crystal::AtomType;
use Crystal::AtomSite;
use Crystal::CIF;

#===============================================
# スクリプト大域変数
#===============================================


#===============================================
# 変数取得関数
#===============================================
BEGIN {}

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

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

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

sub DESTROY
{
	my $this = shift;
}

#============================================================
#　継承クラスで定義しなおす関数
#============================================================
sub CheckFileType
{
	my ($path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
#	my $in = new JFile($path, "r");
#	return undef unless($in);
#	my $line1 = $in->ReadLine();
#	$in->Close();
	if($filename =~ /BandStr\.castep$/i) {
		return "CASTEP Band Output file";
	}
	return undef;
}

sub ReadCASTEPBandOutput
{
	my ($this, $path) = @_;
	my ($drive, $dir, $filename, $ext, $lastdir, $filebody)
		= Deps::SplitFilePath($path);
	my $SampleName = $filebody;

print "  Read [$path].\n";
	my $in = new JFile($path, "r");
	unless($in) {
		print "  Error: Can not ead [$path].\n";
		return undef;
	}

	my $line = $in->SkipTo("number of  electrons");
	my ($nTotalElectrons) = ($line =~ /:\s*(\d+)/);
	$line = $in->SkipTo("net charge of system");
	my ($NetCharge) = ($line =~ /:\s*(\d+)/);
	$line = $in->SkipTo("number of\\s+up\\s+spins");
	my ($nUpSpins) = ($line =~ /:\s*(\d+)/);
	$line = $in->SkipTo("number of\\s+down\\s+spins");
	my ($nDownSpins) = ($line =~ /:\s*(\d+)/);
	$line = $in->SkipTo("\\+  Fermi energy for");
	my ($EF) = ($line =~ /\s([\+\-\d\.eE]+)\s*eV/);
print "nTotalElectrons=$nTotalElectrons\n";
print "NetCharge =$NetCharge\n";
print "nUpSpins  =$nUpSpins\n";
print "nDownSpins=$nDownSpins\n";
print "EF=$EF\n";

#番度数、
	my $nLevels = 0;
	my $nBand = 0;
	while(1) {
		$line = $in->SkipTo("Spin=");
		last if($in->eof());
		my ($Spin, $kpt, $kx, $ky, $kz, $kptgrp)
			= ($line =~ /Spin=(\d+)\s+kpt=\s+(\d+)\s\(\s*([\d\+\-\.]+)\s*([\d\+\-\.]+)\s*([\d\+\-\.]+)\s*\)\s+kpt\-group=\s*(\d)+[\s\+]/);
		last if(not defined $Spin);
#print "Spin=$Spin kpt=$kpt ($kx $ky $kz) kpt-group=$kptgrp\n";

#  + ----------------------------------------------------------------- +
		$in->ReadLine();
#  +                                                                   +
		$in->ReadLine();
		$nLevels = 0;
		while(1) {
#  +           1      -18.605352                                       +
			$line = $in->ReadLine();
			last if($in->eof() or not defined $line);
			last if($line =~ /\+ -----/);
			my ($iLevel, $Energy) = ($line =~ /\+\s+(\S+)\s+(\S+)/);
			last if(not defined $iLevel);
#print "$iLevel: $Energy\n";
			$nLevels++;
		}
		$nBand++;
	}
print "nBands : $nBand\n";
print "nLevels: $nLevels\n";

	my @BandArray;
	for(my $i = 0 ; $i < $nLevels ; $i++) {
		$BandArray[$i] = new EnergyBand;
	}

	$in->rewind();

	my @BandBoundaryDistance;
	my @BoundaryPositions;
	my $iBand;
	my $iLevel;
	my ($prevkx, $prevky, $prevkz);
	my $prevkptgrp;
	my $kDistance = 0.0;
	while(1) {
#  +  Spin=1 kpt=    1 (  0.000000  0.000000  0.000000) kpt-group=  1  +
		$line = $in->SkipTo("Spin=");
		last if($in->eof());
		my ($Spin, $kpt, $kx, $ky, $kz, $kptgrp)
			= ($line =~ /Spin=(\d+)\s+kpt=\s+(\d+)\s\(\s*([\d\+\-\.]+)\s*([\d\+\-\.]+)\s*([\d\+\-\.]+)\s*\)\s+kpt\-group=\s*(\d)+[\s\+]/);
		last if(not defined $Spin);
print "Spin=$Spin kpt=$kpt ($kx $ky $kz) kpt-group=$kptgrp\n";

		my $dk = 0.0;
		if(defined $prevkx) {
			$dk = ($kx-$prevkx)*($kx-$prevkx) + ($ky-$prevky)*($ky-$prevky) + ($kz-$prevkz)*($kz-$prevkz);
			$dk = sqrt($dk);
			$kDistance += $dk;
			if($kptgrp > $prevkptgrp) {
				push(@BandBoundaryDistance, $kDistance);
				push(@BoundaryPositions, "$kx $ky $kz");
				$prevkptgrp = $kptgrp;
			}
		}
		else {
			push(@BandBoundaryDistance, 0.0);
			push(@BoundaryPositions, "$kx $ky $kz");
			$prevkptgrp = $kptgrp;
		}
		$prevkx = $kx;
		$prevky = $ky;
		$prevkz = $kz;
#  + ----------------------------------------------------------------- +
		$in->ReadLine();
#  +                                                                   +
		$in->ReadLine();
		$iLevel = 0;
		while(1) {
#  +           1      -18.605352                                       +
			$line = $in->ReadLine();
			last if($in->eof() or not defined $line);
			last if($line =~ /\+ -----/);
			my ($iLevel, $Energy) = ($line =~ /\+\s+(\S+)\s+(\S+)/);
			last if(not defined $iLevel);
#print "$iLevel: $Energy\n";

			$Energy -= $EF;
			my $ne = 2.0;
			$ne = 0.0 if($Energy > 0.0);

			$BandArray[$iLevel-1]->AddByKDEN($kx, $ky, $kz, $kDistance, $Energy, $ne);
		}
		$iBand++;
	}
	$in->Close();

	$this->SetDataArray(new GraphDataArray);
	my $Data = new GraphData;
	$Data->SetTitle($SampleName);

	my $band = $BandArray[0];
	for(my $i = 0 ; ; $i++) {
		$band = $BandArray[$i];
		last unless($band);
		$band->CalMinMax();
		$Data->SetXDataArray($i, $band->pDistance());
		$Data->SetYDataArray($i, $band->pEnergy());
		$Data->{"Ne$i"} = $band->pNe();
		my ($ymin, $ymax) = $band->GetYMinMax();
#		my $Ne = $band->Ne(0); #Ne($i);
#		my $s = "Band $i ($ymin - $ymax) Ne=$Ne";
		my $s = "Band $i ($ymin - $ymax)";
		$Data->SetXName(0, "k");
		$Data->SetYName($i, $s);
	}

	$Data->CalMinMax();
	$this->DataArray()->AddGraphData($Data);

	return $path;
}

sub Read
{
	my ($this, $filename) = @_;
	$this->ClearAll();
	my $FileType = $this->{'FileType'} = CheckFileType($filename);
	return undef unless($FileType);

	$this->SetFileName($filename);
	if($FileType eq "CASTEP Band Output file") {
		return $this->ReadCASTEPBandOutput($filename);
	}
	return undef;
}

#===============================================
# 一般関数
#===============================================
sub SaveCARFile
{
	my ($this, $Crystal, $filename, $IsChooseRandomly) = @_;
	my $LF = "<br>\n";

#ファイル作製開始
	unless(open(OUT,">$filename")) {
		print "Can not write to $filename.$LF$LF";
		return;
	}
	print OUT "!BIOSYM archive 3\n";
	print OUT "PBC=ON\n";

	my ($sec, $min, $hour, $mday, $mon, $year, $wday, $yday, $isdst) 
		= localtime(time());
	my @wDayStr = ("Sun", "Mon", "Tue", "Wed", "Thu", "Fri", "Sat");
	my @MonStr  = ("Jan", "Feb", "Mar", "Apr", "May", "Jun", 
		       "Jul", "Aug", "Sep", "Oct", "Nov", "Dec");
	printf OUT "!DATE %3s %3s %2d %2d:%2d:%2d %4d\n",
		$wDayStr[$wday], $MonStr[$mon], $mday, 
		$hour, $min, $sec, $year+1900;

	my $SPGName = $Crystal->SPGNameByOutputMode();
	$SPGName = 'P1' if($SPGName eq 'P');
	my ($a,$b,$c,$alpha,$beta,$gamma) 
		= $Crystal->LatticeParametersByOutputMode(0);
	printf OUT "PBC  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f  %8.4f (%s)\n",
		$a, $b, $c, $alpha, $beta, $gamma, $SPGName;

	my @ExpandedAtomSiteList = $Crystal->GetCExpandedAtomSiteListByOutputMode();

#Occupancyが1のサイト
	my %AtomCount;
	my $count = 0;
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		next if($occupancy < 0.9999);

		$AtomCount{$atomname}++;
		$count++;

		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
		printf OUT "%-5s %14.9f %14.9f %14.9f %-4s %-7s%-7s %-2s %6.3f\n",
			"$atomname$AtomCount{$atomname}", $xc, $yc, $zc, 
			'XXXX', '1', 'xx', $atomname, $charge;
	}
#Occupancyが1未満のサイト
	if(@ExpandedAtomSiteList > $count and !$IsChooseRandomly) {
		print "<H3>***Warning***: One line should be removed from $filename.</H3>\n";
		print OUT "*** Shared atoms: Remove this line for calculation: OCC=";
		for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
			my $atom      = $ExpandedAtomSiteList[$i];
			my $occupancy = $atom->Occupancy();
			next if($occupancy >= 0.9999);
			printf OUT "%lg ", $occupancy;
		}
		print OUT "\n";
	}
	for(my $i = 0 ; $i < @ExpandedAtomSiteList ; $i++) {
		my $atom      = $ExpandedAtomSiteList[$i];
		my $atomname  = $atom->AtomNameOnly();
		my $charge    = $atom->Charge();
		my ($x,$y,$z) = $atom->Position(1);
		my $occupancy = $atom->Occupancy();
		next if($occupancy >= 0.9999);

		$AtomCount{$atomname}++;

		my ($xc,$yc,$zc) = $Crystal->FractionalToCartesian($x,$y,$z);
		printf OUT "%-5s %14.9f %14.9f %14.9f %-4s %-7s%-7s %-2s %6.3f\n",
			"$atomname$AtomCount{$atomname}", $xc, $yc, $zc, 
			'XXXX', '1', 'xx', $atomname, $charge;
	}

	print OUT "end\n";
	print OUT "end\n";
	print OUT "\n";

#	&ListPartialOccupancyAtoms();

	close(OUT);
}

1;