#!/usr/bin/perl

BEGIN {
#use lib 'd:/Programs/Perl/lib';
#use lib '/home/tkamiya/bin/lib';
my $BaseDir = $ENV{'TkPerlDir'};
print "\n\nBaseDir: $BaseDir\n";
@INC = ("$BaseDir/lib", "$BaseDir/VNL", "d:/Programs/Perl/lib", @INC);
}

use strict;
#use warnings;
use File::Path;
use File::Basename;
use File::Find;

use Math::Matrix;
use Math::MatrixReal;

use Deps;
use Utils;
use JFile;

use MyApplication;
use Sci::Algorism;
use Sci::GeneralFileFormat;

use Crystal::CIF;
use Crystal::Crystal;
use Crystal::SpaceGroup;
use Crystal::AtomType;
use Crystal::VASP;
use Crystal::Quantum;

#===============================================
# デバッグ関係変数
#===============================================
#$PrintLevelが大きいほど、情報が詳しくなる
my $PrintLevel = 0;

#===============================================
# グローバル変数
#===============================================
my $a0 = Sci::a0() * 1.0e10; #A
my $pi = Sci::pi();

#===============================================
# 文字コード関係変数
#===============================================
# sjis, euc, jis, noconv
my $PrintCharCode      = Deps::PrintCharCode();
my $OSCharCode         = Deps::OSCharCode();
my $FileSystemCharCode = Deps::FileSystemCharCode();

my $LF        = Deps::LF();
my $DirSep    = Deps::DirSep();
my $RegDirSep = Deps::RegDirSep();

#===============================================
# Applicationオブジェクト作成
#===============================================
my $App = new MyApplication;
exit if($App->Initialize() < 0);

#$App->SetLF("<br>\n");
#$App->SetPrintCharCode("sjis");
#$App->SetDebug($Debug);
$App->SetDeleteHTMLFlag(1);

#===============================================
# スクリプト大域変数
#===============================================
my $InitialDirectory = Deps::GetWorkingDirectory();

#==========================================
# コマンドラインオプション読み込み
#==========================================
$App->AddArgument("--Action",
		"--Action=[Calc3DElectronDensity|CalcR|CalcSTO]",       '');
$App->AddArgument("--n",              "   --n=[Integer] (Def:1)",              1);
$App->AddArgument("--l",                  "--l=[Integer] (Def:1)",             1);
$App->AddArgument("--m",                  "--m=[Integer] (Def:1)",             1);
$App->AddArgument("--Z",                  "--Z=[val] (Def:1.0)",             1.0);
$App->AddArgument("--Rmax",               "--Rmax=[val] (Def:5.0)",          5.0);
$App->AddArgument("--nRMesh",             "--nRMesh=[val] (Def:201)",        201);
$App->AddArgument("--lx",                 "--lx=[value] (Def:10.0A)",        10.0);
$App->AddArgument("--ly",                 "--ly=[value] (Def:10.0A)",        10.0);
$App->AddArgument("--lz",                 "--lz=[value] (Def:10.0A)",        10.0);
$App->AddArgument("--lstep",              "--lstep=[value] (Def:0.2A)",       0.2);
$App->AddArgument("--IonName",            "--IonName=[name] (Def:Na+)",     "Na+");
$App->AddArgument("--UseExplicitFormula", "--UseExplicitFormula=[0|1] (Def:0)", 0);

exit 1 if($App->ReadArgs(0) != 1);
my $Args = $App->Args();
#my $form = new CGI;
#$Args->SetCGIForm($form);
#$Args->parseInput($WebCharCode);

#==========================================
# メイン関数スタート
#==========================================

#Utils::InitHTML("Research", $WebCharSet, "_self");

my $Debug = $Args->GetGetArg("DebugMode");
$App->SetDebug($Debug);
my $Action             = $Args->GetGetArg("Action");
my $UseExplicitFormula = $Args->GetGetArg("UseExplicitFormula");

my $ret = 0;
if($Action =~ /Calc3DElectronDensity/i) {
#	&CalR();
	&Make3DDensityMap();
}
elsif($Action =~ /CalcR/i) {
	&CalcR();
}
elsif($Action =~ /CalcSTO/i) {
	&CalcSTO();
}
else {
	$App->print("Error: Invald Action: $Action\n");
}

#Utils::EndHTML();

exit $ret;

#===============================================
# スクリプト終了
#===============================================

#==========================================
# &Subroutines
#==========================================

sub CalcSTO
{
	my $IonName = $Args->GetGetArg("IonName");
	$IonName = "Na+" unless(defined $IonName);

	$App->print("Ion Name: $IonName\n");
	my ($AtomicNumber, $AtomName, $Charge) = AtomType::GetAtomInformation($IonName);
	$App->print("Atom Name: $AtomName\n");
	$App->print("Charge   : $Charge\n");
	$App->print("Atomic Number: $AtomicNumber\n");

#	my ($zStar, $nStar, $n1, $Rmax) = Quantum->CalSTO($IonName);
#return;

	my $nElectron = $AtomicNumber - $Charge;
	my $ResCharge = $nElectron;
	my $n1s = 0.0;
	if($ResCharge > 2.0) {
		$n1s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n1s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n2s = 0.0;
	if($ResCharge > 2.0) {
		$n2s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n2s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n2p = 0.0;
	if($ResCharge > 6.0) {
		$n2p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n2p = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n3s = 0.0;
	if($ResCharge > 2.0) {
		$n3s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n3s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n3p = 0.0;
	if($ResCharge > 6.0) {
		$n3p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n3p = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n3d = 0.0;
	my $n4s = 0.0;
	if($Charge == 0.0) {
		if($ResCharge > 2.0) {
			$n4s = 2.0;
			$ResCharge -= 2.0;
		}
		else {
			$n4s = $ResCharge;
			$ResCharge = 0.0;
		}
		if($ResCharge > 10.0) {
			$n3d = 10.0;
			$ResCharge -= 10.0;
		}
		else {
			$n3d = $ResCharge;
			$ResCharge = 0.0;
		}
	}
	else {
		if($ResCharge > 10.0) {
			$n3d = 10.0;
			$ResCharge -= 10.0;
		}
		else {
			$n3d = $ResCharge;
			$ResCharge = 0.0;
		}
		if($ResCharge > 2.0) {
			$n4s = 2.0;
			$ResCharge -= 2.0;
		}
		else {
			$n4s = $ResCharge;
			$ResCharge = 0.0;
		}
	}
	my $n4p = 0.0;
	if($ResCharge > 6.0) {
		$n4p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n4p = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n4d = 0.0;
	if($ResCharge > 10.0) {
		$n4d = 10.0;
		$ResCharge -= 10.0;
	}
	else {
		$n4d = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n5s = 0.0;
	if($ResCharge > 2.0) {
		$n5s = 2.0;
		$ResCharge -= 2.0;
	}
	else {
		$n5s = $ResCharge;
		$ResCharge = 0.0;
	}
	my $n5p = 0.0;
	if($ResCharge > 6.0) {
		$n5p = 6.0;
		$ResCharge -= 6.0;
	}
	else {
		$n5p = $ResCharge;
		$ResCharge = 0.0;
	}

	my $nTotal = $n1s + $n2s + $n2p + $n3s + $n3p + $n3d + $n4s + $n4p + $n4d + $n5s + $n5p;
$App->print("Total: $nTotal\n");
$App->print(" 1s: $n1s\n");
$App->print(" 2s: $n2s\n");
$App->print(" 2p: $n2p\n");
$App->print(" 3s: $n3s\n");
$App->print(" 3p: $n3p\n");
$App->print(" 3d: $n3d\n");
$App->print(" 4s: $n4s\n");
$App->print(" 4p: $n4p\n");
$App->print(" 4d: $n4d\n");
$App->print(" 5s: $n5s\n");
$App->print(" 5p: $n5p\n");
	if($ResCharge > 0.0) {
		$App->print("Error: Ne=$AtomicNumber - $Charge is not supported.\n");
		exit;
	}
	if($nTotal ne $nElectron) {
		$App->print("Error: Inconsistent Ne (Ntotal=$nTotal != Ne=$nElectron).\n");
		exit;
	}

	my @ne = ($n1s, $n2s+$n2p, $n3s+$n3p+$n3d, $n4s+$n4p+$n4d, $n5s+$n5p);
	my $MaxN = 0.0;
	for(my $i = 0 ; $i < @ne ; $i++) {
		if($ne[$i] > 0.0) {
			$MaxN = $i+1;
		}
	}
$App->print("Max n: $MaxN\n");
	my @nStar = (undef, 1, 2, 3, 3.7, 4.0, 4.2);
	my $nStar = $nStar[$MaxN];
$App->print("n*   : $nStar\n");

	my $Is_d = 0;
	if(($n3d > 0.0 and $n4s == 0.0) or ($n4d > 0.0 and $n5s == 0.0)) {
		$Is_d = 1;
	}

#(i)  z*=z-sとする。ここでzは核電荷、sは遮蔽定数
#(ii) 電子を(1s)(2s,2p)(3s,3p)(3d)･････という群にわける。
#	my @nGroup = ($n1s, $n2s+$n2p, $n3s+$n3p, $n3d, $n4s+$n4p, $n4d, $n5s+$n5p);
	my @nGroup = ($n1s, $n2s+$n2p, $n3s+$n3p+$n3d, $n4s+$n4p+$n4d, $n5s+$n5p);
	my $nMax = 0;
	for(my $i = 0 ; $i < @nGroup ; $i++) {
		if($nGroup[$i] > 0.0) {
			$nMax = $i;
		}
	}
$App->print("n(max) for groups: $nMax\n");

	my $stotal = 0.0;
#(iii) 考えている電子より主量子数nの大きい電子からのsへの寄与はゼロとする。
#(iv)  考えている電子の属する群の他の電子からはそれぞれ0.35とする。ただし、1s群の場合は0.30とする。
	my $s;
	my $n = $nGroup[$nMax]-1.0;
	my $k;
	if($nGroup[$nMax] > 1.0) {
		if($nMax == 0) { # 1s群
			$k = 0.30;
		}
		else {
			$k = 0.35;
		}
		$s = $n * $k;
	}
$App->print("s(outmost)  = $n * $k = $s\n");
	$stotal += $s;

#(v) 考えている電子がsまたはp電子なら、それよりnの一つ小さい電子からは0.85ずつ
#ただし、d電子ならば、それより内側の電子からの寄与はすべて1.00とする。
	$n = $nGroup[$nMax-1];
	if($nMax >= 1 and $nGroup[$nMax-1] > 1.0) {
#		if($nMax == 3) { # d群
		if($Is_d) { # d群
			$k = 1.0;
		}
		else {
			$k = 0.85;
		}
		$s = $n * $k;
	}
$App->print("s(n(max)-1) = $n * $k = $s\n");
	$stotal += $s;

#もっと内側の電子からは1.00ずつ寄与する。ただし、d電子ならば、それより内側の電子からの寄与はすべて1.00とする。
	$k = 1.0;
	$n = 0.0;
	for(my $i = $nMax - 2 ; $i >= 0 ; $i--) {
			$n += $nGroup[$i] * 1.0;
	}
	$s = $n * $k;
	$stotal += $s;
$App->print("s(inner)    = $n * $k = $s\n");
$App->print("s(total)    = $stotal\n");

	my $zStar = $AtomicNumber - $stotal;
	my $Rmax  = $nStar*$nStar / $zStar * Sci::a0() * 1.0e10;
	my $n1    = $nStar - 1;

$App->print("Z    : $AtomicNumber\n");
$App->print("Z*   : $zStar\n");
$App->print("R(r) = N * r^$n1 * exp[-($zStar/$nStar a0) * r]\n");
$App->print("     a0=$a0 angstrom\n");
$App->print("Rmax = (a0/Z*)(n*)^2 = $Rmax (A)\n");
}

sub GetOrbName
{
	my ($n, $l, $m) = @_;
	my @names1 = ('s'); 						# for l=0
	my @names2 = ('py', 'pz', 'px');				# for l=1, m=-1, 0, 1
	my @names3 = ('dxy', 'dyz', 'd3z2-r2', 'dzx', 'dx2-y2');	# for l=2, m=-2,-1,0,1,2
	my @names4 = ('fy(3x2-y2)', 'fxyz', 'fy(5z2-r2)', 'fz(5z2-3r2)',
		      'fx(5z2-r2)', 'fz(x2-y2)', 'fx(x2-3y2)');		# for l=3, m=-3,-2,-1,0,1,2,3
	my @names = (\@names1, \@names2, \@names3, \@names4);

#print "names=@names\n";
#print "l=$l  m+l=", $m+$l, "\n";
	my $pNames = $names[$l];
#print "pNames=$pNames\n";
	my $name = $pNames->[$m+$l];
#print "name=$name\n";
	return "$n$name";
}

sub CalRadialWaveFunction
{
	my ($Z, $n, $l, $m, $x, $y, $z, $cx, $cy, $cz) = @_;

	my $dx = $x - $cx;
	my $dy = $y - $cy;
	my $dz = $z - $cz;
	my $r2 = $dx*$dx + $dy*$dy + $dz*$dz;
	my $r  = sqrt($r2);
	$r = 0.001 if($r < 0.001);

	if($UseExplicitFormula) {
		return Quantum->CalRadialWaveFunction($Z, $n, $l, $m, $r);
	}
	return Quantum->CalRadialWaveFunction2($Z, $n, $l, $m, $r);




	my $k = 1.0 / sqrt($pi) * sqrt(($Z/$a0)^3);
	my $Z32 = $Z**1.5;
	my $zeta = $Z / $n;
	my $rho  = 2.0 * $zeta * $r;
	my $Rr;
	if($n == 1 and $l == 0) {	# 1s
		$Rr = $Z32 * 2.0 * exp(-$rho/2.0);
	}
	elsif($n == 2 and $l == 0) {	# 2s
		$Rr = 0.5 * sqrt(2.0) * $Z32 * (2.0-$rho) * exp(-$rho/2.0);
	}
	elsif($n == 2 and $l == 1) {	# 2p
		$Rr = 0.5 * sqrt(6.0) * $Z32 * $rho * exp(-$rho/2.0);
	}
	elsif($n == 3 and $l == 0) {	# 3s
		$Rr = 1.0/9.0 * sqrt(3.0) * $Z32 * (6.0-6.0*$rho+$rho*$rho) * exp(-$rho/2.0);
	}
	elsif($n == 3 and $l == 1) {	# 3p
		$Rr = 1.0/9.0 * sqrt(6.0) * $Z32 * (4.0-$rho)*$rho * exp(-$rho/2.0);
#print "rho=$rho R=$Rr\n";
	}
	elsif($n == 3 and $l == 2) {	# 3d
		$Rr = 1.0/9.0 * sqrt(30.0) * $Z32 * $rho*$rho * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 0) {	# 4s
		$Rr = 1.0/96.0 * $Z32 * (24.0-36.0*$rho+12.0*$rho*$rho-$rho*$rho*$rho) * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 1) {	# 4p
		$Rr = 1.0/32.0 * sqrt(15.0) * $Z32 * (20.0-10.0*$rho+$rho*$rho)*$rho * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 2) {	# 4d
		$Rr = 1.0/96.0 * sqrt(5.0) * $Z32 * (6.0-$rho)*$rho*$rho * exp(-$rho/2.0);
	}
	elsif($n == 4 and $l == 3) {	# 4f
		$Rr = 1.0/96.0 * sqrt(35.0) * $Z32 * $rho*$rho*$rho * exp(-$rho/2.0);
	}
	else {
		return undef;
	}
#print "x=$x  r=$r rho=$rho R=$Rr\n";
	return $Rr;
}

sub CalWaveFunction
{
	my ($Z, $n, $l, $m, $x, $y, $z, $cx, $cy, $cz) = @_;

	my $dx = $x - $cx;
	my $dy = $y - $cy;
	my $dz = $z - $cz;
	my $r2 = $dx*$dx + $dy*$dy + $dz*$dz;
	my $r  = sqrt($r2);
	if($r < 0.001) {
		$r = 0.001;
	}

	my $cosQ   = $dz / $r;
	my $Q      = Sci::acos($cosQ);
	my $sinQ   = sin($Q);
	my $cosPhi = 0.0;
	if(abs($sinQ) < 0.001) {
		$cosPhi = 1.0;
	}
	else {
		$cosPhi = $dx / $r / $sinQ;
	}
	my $Phi    = Sci::acos($cosPhi);
	   $Phi    = -$Phi if($dy < 0.0);

	my $Rr = &CalRadialWaveFunction($Z, $n, $l, $m, $x, $y, $z, $cx, $cy, $cz);
	unless(defined $Rr) {
		return undef;
	}

	if(!$UseExplicitFormula) {
		my $Quantum = new Quantum;

#		my $Theta = $Quantum->CalTheta($l, $m, $Q);
		my $Theta = $Quantum->CalTheta2($l, $m, $Q);
		my $Phi   = $Quantum->CalPhi($m, $Phi);
		return $Rr * $Theta * $Phi;
	}

	my $k = 1.0 / sqrt($pi) * sqrt(($Z/$a0)^3);
#	my $sinPhi = 1.0 - $cosPhi*$cosPhi;
#	$sinPhi = 0.0 if($sinPhi < 0.0);
#	$sinPhi = sqrt($sinPhi);
	my $sinPhi = sin($Phi);
	my $rho    = $Z/$a0 * $r;

	my $val = 0.0;
	if($l == 0) {	# s
		$val = sqrt(1.0/4.0/$pi) * $Rr;
	}
	elsif($l == 1 and $m == 0) {	# pz
		$val = sqrt(3.0/4.0/$pi) * $cosQ * $Rr;
#print "cosQ=$cosQ  val=$val\n";
	}
	elsif($l == 1 and $m == -1) {	# px
		$val = sqrt(3.0/4.0/$pi) * $sinQ * $sinPhi * $Rr;
	}
	elsif($l == 1 and $m == +1) {	# py
		$val = -sqrt(3.0/4.0/$pi) * $sinQ * $cosPhi * $Rr;
	}
	elsif($l == 2 and $m == 2) {	# dx2-y2
		$val = sqrt(15.0/16.0/$pi) * $sinQ*$sinQ*($cosPhi*$cosPhi-$sinPhi*$sinPhi) * $Rr;
	}
	elsif($l == 2 and $m == 1) {	# dzx
		$val = -sqrt(15.0/4.0/$pi) * $sinQ*$cosQ*$cosPhi * $Rr;
	}
	elsif($l == 2 and $m == 0) {	# dz2
		$val = sqrt(5.0/16.0/$pi) * (3.0*$cosQ*$cosQ-1.0) * $Rr;
	}
	elsif($l == 2 and $m == -1) {	# dyz
		$val = -sqrt(15.0/4.0/$pi) * $sinQ*$cosQ*$sinPhi * $Rr;
	}
	elsif($l == 2 and $m == -2) {	# dxy
		$val = sqrt(15.0/4.0/$pi) * $sinQ*$sinQ*$sinPhi*$cosPhi * $Rr;
	}
	elsif($l == 3 and $m == 3) {	# fx(x2-3y2)
		$val = -sqrt(35.0/32.0/$pi) * $sinQ*$sinQ*$sinQ * 
				$cosPhi * ($cosPhi*$cosPhi-3.0*$sinPhi*$sinPhi) * $Rr;
	}
	elsif($l == 3 and $m == 2) {	# fz(x2-y2)
		$val = sqrt(105.0/16.0/$pi) * $sinQ*$sinQ*$cosQ * 
				($cosPhi*$cosPhi-$sinPhi*$sinPhi) * $Rr;
	}
	elsif($l == 3 and $m == 1) {	# fx(5z2-r2)
		$val = -sqrt(21.0/32.0/$pi) * $sinQ * 
				$cosPhi * (5.0*$cosQ*$cosQ-1.0) * $Rr;
	}
	elsif($l == 3 and $m == 0) {	# fz(5z2-3r2)
		$val = sqrt(7.0/16.0/$pi) * $cosQ*(5.0*$cosQ*$cosQ-3.0) * $Rr;
	}
	elsif($l == 3 and $m == -1) {	# fy(5z2-r2)
		$val = sqrt(21.0/31.0/$pi) * $sinQ*(5.0*$cosQ*$cosQ-1.0) * 
				$cosPhi * $Rr;
	}
	elsif($l == 3 and $m == -2) {	# fxyz
		$val = sqrt(105.0/4.0/$pi) * $sinQ*$sinQ*$cosQ * 
				$sinPhi*$cosPhi * $Rr;
	}
	elsif($l == 3 and $m == -3) {	# fy(3x2-y2)
		$val = -sqrt(35.0/32.0/$pi) * $sinQ*$sinQ*$sinQ * 
				$sinPhi*(3.0*$cosPhi*$cosPhi-$sinPhi*$sinPhi) * $Rr;
	}

	return $val;
}

sub Make3DDensityMap
{
$App->printf("\nCalc 3D Electron Density map\n");

	my $OutFile = $Args->GetGetArg(0);
	   $OutFile = "EDens.vasp" unless(defined $OutFile);
	my $Z       = $Args->GetGetArg("Z");
	   $Z = 1.0 unless(defined $Z);
	my $n       = $Args->GetGetArg("n");
	   $n = 1 unless(defined $n);
	my $l       = $Args->GetGetArg("l");
	   $l = 0 unless(defined $l);
	my $m       = $Args->GetGetArg("m");
	   $m = 0 unless(defined $m);

	my $OrbName = Quantum->GetOrbName($n, $l, $m);

	$App->print("  Z=$Z\n");
	$App->print("  n,l,m=$n,$l,$m\n");
	$App->print("  Orbital: $OrbName\n");
	$OutFile =~ s/{Orb}/$OrbName/i;
	$App->print("OutFile: [$OutFile]\n");

	my $out = new JFile($OutFile, "w");
	if(!defined $out) {
		$App->print("Error: Can not write to [$OutFile].\n");
		return -1;
	}

	my $lx      = $Args->GetGetArg("lx"); # A
	   $lx = 10.0 unless(defined $lx);
	my $ly      = $Args->GetGetArg("ly");
	   $ly = 10.0 unless(defined $ly);
	my $lz      = $Args->GetGetArg("lz");
	   $lz = 10.0 unless(defined $lz);
	my $step    = $Args->GetGetArg("lstep"); # A
	   $step = 0.2 unless(defined $step);
#	my $nx = int($lx / $step + 1.0);
#	my $ny = int($ly / $step + 1.0);
#	my $nz = int($lz / $step + 1.0);
	my $nx = int($lx / $step + 0.1);
	my $ny = int($ly / $step + 0.1);
	my $nz = int($lz / $step + 0.1);
$App->print("Range: $lx $ly $lz\n");
$App->print("nMesh: $nx $ny $nz\n");

	$out->print("Electron density for Z=$Z, n=$n, l=$l, m=$m\n");
	$out->printf("   %12.6f\n", $lx);
	$out->printf("     %12.6f   %12.6f   %12.6f\n", $lx / $lx, 0.0,       0.0);
	$out->printf("     %12.6f   %12.6f   %12.6f\n", 0.0,       $ly / $lx, 0.0);
	$out->printf("     %12.6f   %12.6f   %12.6f\n", 0.0,       0.0,       $lz / $lx);
	$out->printf("  1\n");
	$out->print("Direct\n");
	$out->print("  0.500000  0.500000  0.500000\n");
	$out->print("\n");
	$out->printf("  %3d  %3d  %3d\n", $nx, $ny, $nz);

	my $Count = 0;
	for(my $ix = 0 ; $ix < $nx ; $ix++) {
	for(my $iy = 0 ; $iy < $ny ; $iy++) {
	for(my $iz = 0 ; $iz < $nz ; $iz++) {
		my $x = $ix * $step;
		my $y = $iy * $step;
		my $z = $iz * $step;
#print "(x,y,z) = ($x, $y, $z)\n";

#中心は(0.5$lx, 0.5$ly, 0.5$lz)
		my $f = &CalWaveFunction($Z, $n, $l, $m, $x, $y, $z, 0.5*$lx, 0.5*$ly, 0.5*$lz);
		unless(defined $f) {
			$App->print("Invalid n,l,m($n,$l,$m).\n");
			$out->Close();
			return -2;
		}

		$out->printf(" %11.5g", $f);
		$Count++;
		$out->printf("\n") if($Count % 10 == 0);
	}
	}
	}

	$out->Close();
}

sub CalcR
{
$App->printf("\nCalc Radial Wave Function\n");

	my $OutFile = $Args->GetGetArg(0);
	   $OutFile = "EDens.vasp" unless(defined $OutFile);
	my $Z       = $Args->GetGetArg("Z");
	   $Z = 1.0 unless(defined $Z);
	my $n       = $Args->GetGetArg("n");
	   $n = 1 unless(defined $n);
	my $l       = $Args->GetGetArg("l");
	   $l = 0 unless(defined $l);
	my $m       = $Args->GetGetArg("m");
	   $m = 0 unless(defined $m);
	my $Rmax    = $Args->GetGetArg("Rmax");
	   $Rmax = 0 unless(defined $Rmax);
	my $nRMesh  = $Args->GetGetArg("nRMesh");
	   $nRMesh = 0 unless(defined $nRMesh);

	my $OrbName = Quantum->GetOrbName($n, $l, $m);

	$App->print("  Rmax  =$Rmax\n");
	$App->print("  nRMesh=$nRMesh\n");
	$App->print("  Z=$Z\n");
	$App->print("  n,l,m=$n,$l,$m\n");
	$App->print("  Orbital: $OrbName\n");
	$OutFile =~ s/{Orb}/${OrbName}-R/i;
	$OutFile =~ s/\.vasp/\.csv/i;
	$App->print("OutFile: [$OutFile]\n");

	my $out = new JFile($OutFile, "w");
	if(!defined $out) {
		$App->print("Error: Can not write to [$OutFile].\n");
		return -1;
	}
	$out->print("r(A),R(r),r^2R(r)^2\n");

	my $rmin = 0.0;
	
	my $RStep = ($Rmax - $rmin) / ($nRMesh-1);
	for(my $i = 0 ; $i < $nRMesh ; $i++) {
		my $r = $rmin + $i * $RStep;
#print "#i=$i r=$r\n";
		my $Rr = &CalRadialWaveFunction($Z, $n, $l, $m, $r, 0.0, 0.0, 0.0, 0.0, 0.0);
		unless(defined $Rr) {
			return undef;
		}
		my $Dens = $r*$r * $Rr*$Rr;
#$App->printf("%3d: %12.6g  %12.6g  %12.6g\n", $i, $r, $Rr, $Dens);
		$out->printf("%12.6g,%12.6g,%12.6g\n", $r, $Rr, $Dens);
	}

	$out->Close();
}