#!/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 Utils;
use JFile;

use Sci qw($a0 $pi);
use Sci::Algorism;
use Crystal::Quantum;

my $a0InA = $a0 * 1.0e10;

#===============================================
# グローバル変数
#===============================================
my $OutFile = "out.csv";

my $Z = 5;
my $n = 4;
my $l = 3;
my $m = 0;
my ($rmin, $Rmax, $nRMesh) = (0.0, 7.0, 7001);
my $RStep = ($Rmax - $rmin) / ($nRMesh-1);

my $Quantum = new Quantum;

#==========================================
# メイン関数スタート
#==========================================
&CalcR();
exit;

#==========================================
# Subroutines
#==========================================
sub CalcR
{
	my $OrbName = $Quantum->GetOrbName($n, $l, $m);

	print("OutFile: [$OutFile]\n");
	print("  Rmax  =$Rmax\n");
	print("  nRMesh=$nRMesh\n");
	print("  Z=$Z\n");
	print("  n,l,m=$n,$l,$m\n");
	print("  Orbital: $OrbName\n");

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

	$out->print("r(A),RA(r),r2RA2(r),SA(r),R(r),r2R2(r),S(r)\n");
	my $SA = 0.0;
	my $S  = 0.0;
	for(my $i = 0 ; $i < $nRMesh ; $i++) {
		my $r     = $rmin + $i * $RStep;
		my $rho   = 2.0 * $r / $a0InA * $Z / $n;

		my $RrA   = $Quantum->CalRadialWaveFunction($Z, $n, $l, $m, $r);
		my $DensA = $r**2 * $RrA*$RrA;
		$SA += $DensA * $RStep;

		my $Rr     = $Quantum->CalRadialWaveFunction2($Z, $n, $l, $m, $r);
#		my $L      = Sci::LaguerrePolynomial2($n+$l, 2.0*$l+1, $rho);
#		my $nl1    = Sci::Factorial($n - $l - 1);
#		my $nl     = Sci::Factorial($n + $l);
#		my $c      = -(2.0 * $Z / $n / $a0InA)**1.5 * sqrt($nl1 / 2.0 / $n / $nl);
#		my $Rr     = $c * $rho**$l * $L * exp(-$rho / 2.0);
#		my $RrNorm = $Rr / $nl;
##		my $RrNorm = ($n * $a0InA / 2.0 / $Z)**(1.5+$l) * $Rr;
##		my $Rr     = $c * $r**$l * $L * exp(-$rho / 2.0);
##		my $RrNorm = $Rr;
#		my $Dens   = $r**2 * $RrNorm*$RrNorm;
		my $Dens   = $r**2 * $Rr*$Rr;
		$S += $Dens * $RStep;

		$out->printf("%12.6g,%12.6g,%12.6g,%12.6g,%12.6g,%12.6g,%12.6g\n", $r, $RrA, $DensA, $SA, $Rr, $Dens, $S);
	}
print "SA = $SA\n";
print "S  = $S\n";

	$out->Close();
}