use lib 'd:/Programs/Perl/lib';
use Sci qw(Gamma);
use Sci::Algorism;

for(my $r = 0.2 ; $r < 5.0 ; $r += 0.2) {
	my $G  = Gamma($r);
#	my $G2 = c_gamma($r);
	my $xi = -2;
	my $F  = &FermiDiracIntegral($xi, $r);
	my $F2 = &FermiDiracIntegralByNumericalIntegration($xi, $r);
#	printf "r=%6.2f: G=%8.4f/%8.4f F=%8.4f\n", $r, $G, $G2, $F;
	printf "r=%6.2f: G=%8.6f F=%8.6f/%8.6f\n", $r, $G, $F, $F2;
}

for(my $x = -5.0 ; $x <= 10.0 ; $x += 0.5) {
	my $r1 = 2.5;
	my $r2 = $r1 - 1.0;
	my $F1 = &FermiDiracIntegral($x, $r1);
	my $F2 = &FermiDiracIntegral($x, $r2);
	printf "x=%6.2f: F($r1)=%8.4f F($r2)=%8.4f R=%8.4f\n", $x, $F1, $F2, $F1 / $F2;
}

# 熱電変換工学−基礎と応用−、編集委員長：阪田亮、REALIZE INC.、2001
sub PoweredFermiDirac
{
	my ($x, $x0, $r) = @_;
	if(!defined $r or $r == 0.0) {
		return 1.0 / (1.0 + exp($x-$x0));
	}
	return $x**$r / (1.0 + exp($x-$x0));
}

# Calculation of the Fermi-Dirac integral FermiIntegral(r, x)

sub FermiDiracIntegralByNumericalIntegration
{
	my ($xi, $r, $h, $xmax) = @_;
	$xmax = 5 * abs($xi)      if(!defined $xmax);
	$h    = $xmax / 5001 if(!defined $h);
	my @y;
	my $i = 0;
	for(my $x = 0 ; $x <= $xmax ; $x += $h) {
		$y[$i] = &PoweredFermiDirac($x, $xi, $r);
#print "$i: $x, $y[$i]\n";
		$i++;
	}
	$i = 2 * int($i / 2) - 1;
#print "n=$i\n";
	return Algorism::SinglePointIntegrateByConstantStepSimpson($i, $h, \@y);
}

sub FermiDiracIntegral
{
	my ($xi, $r) = @_;
	return &DegenerateFermiDirac($xi, $r)    if($xi >= 20.0); # degenerate case
	return &NonDegenerateFermiDirac($xi, $r) if($xi < -0.1);  # non-degenerate case
	return &IntermediateFermiDirac($xi, $r);                  # intermediate case
}

# degenerate case
my @DegenerateFD_z = (0.822467, 0.947033, 0.985551, 0.996233, 0.999040,
					  0.999758, 0.999939, 0.999985, 0.999996, 0.999999);
sub DegenerateFermiDirac
{
	my ($xi, $r) = @_;
	
	my $gi = $xi**($r+1.0) / ($r+1.0);
	my $t = $r;
	my $gg = $t / ($xi*$xi);
	my $gs = $gg * $DegenerateFD_z[0];
	$t -= 1.0;
	for(my $k = 1 ; $k <= 9 ; $k++) {
		$gg *= ($t * ($t-1.0) / ($xi*$xi));
		$gs += ($gg * $DegenerateFD_z[$k]);
		$t -= 2.0;
	}
	
	return $gi * (1.0 + 2.0 * ($r+1.0) * $gs);
}

# non-degenerate case
sub NonDegenerateFermiDirac
{
	my ($xi, $r) = @_;

	my $sum = 0.0;
	my $s = 1.0;
	my $gg;
	do {
		$gg = exp($s * $xi) / $s**($r+1.0);
		$sum += $gg;
		$s += 1.0;
		$gg = exp($s*$xi) / $s**($r+1.0);
		$sum -= $gg;
		$s += 1.0;
	} while(abs($gg) > 1.0e-8);
	
	return $sum * Gamma($r+1.0);
}

sub IntermediateFermiDirac
{
	my ($xi, $r) = @_;

	return &head($xi, $r) + &tail($xi, $r);
}

# Series expansion of the head Hr(xi)
sub head
{
	my ($xi, $r) = @_;
	
	my @a;
	my $g = exp($xi);
	my $gp = $g + 1.0;
	my $gm = $g - 1.0;
	my $gs = $g * $g;
	my $gx = $gp;
	$a[0] = 1.0 / $gx / ($r + 2.0);
	$gx *= $gp;
	$a[1] = $gm / $gx / 2.0 / ($r + 3.0);
	$gx *= $gp;
	$a[2] = ($gs - 4.0*$g + 1.0) / $gx / 6.0 / ($r + 4.0);
	$gx *= $gp;
	$a[3] = $gm * ($gs - 10.0*$g + 1.0) / $gx / 24.0 / ($r + 5.0);
	$gg = 1.0;
	my $sum = 1.0 / ($r + 1.0);
	
	for(my $i = 0 ; $i <= 3 ; $i++) {
		$gg *= 0.1;
		$sum -= $a[i] * $gg;
	}
	
	return $sum * $g / $gp * 0.1**($r+1.0);
}

# Numerical integration of the tail Tr(xi)
sub tail
{
	my ($xi, $r) = @_;

	my $sum = 0.0;
	my $n;
	if($r < 0.0) {
		$n = 128;
	}
	else {
		$n = 64;
	}
	$sum += &FermiDiracSimpson($xi, $r, 0.1, 1.0, $n);

	if($r < 0.0) {
		$n = 512;
	}
	else {
		$n = 256;
	}
	my $lim = 17.0 + 4.5 * $r + $xi;
	$sum += &FermiDiracSimpson($xi, $r, 1.0, $lim, $n);
	return $sum;
}

sub FermiDiracSimpson
{
	my ($xi, $r, $a, $b, $n) = @_;

	my $s2 = 0.0;
	my $s4 = 0.0;
	my $dx = ($b - $a) / $n;
	my $x = $a;
	my $sa = &PoweredFermiDirac($x, $xi, $r);

	for(my $i = 1 ; $i <= $n / 2.0 ; $i++) {
		$x += $dx;
		$s4 += &PoweredFermiDirac($x, $xi, $r);
		$x += $dx;
		$s2 += &PoweredFermiDirac($x, $xi, $r);
	}
	$x = $b;
	$sb = &PoweredFermiDirac($x, $xi, $r);
	return ($sa - $sb + 2.0*$s2 + 4.0*$s4) * $dx / 3.0;
}

# Gamma function
sub c_gamma
{
	my ($x) = @_;

	my $q = 1.0;
	while($x < 5.0) {
		$q *= $x;
		$x += 1.0;
	}
	my $xs = $x * $x;
	my $xx = $x;
	my $p = 0.0833333 / $xx;
	$xx *= $xs;
	$p -= (2.77778e-3 / $xx);
	$xx *= $xs;
	$p += (7.93651e-4 / $xx);
	my $y = ($x - 0.5) * log($x) - $x + $p;
	return 2.50663 * exp($y) / $q;
}