package Sci;
use Common;
@ISA = qw(Common);

#公開したいサブルーチン
#@EXPORT   = qw(erfc tan);
@EXPORT_OK = qw($pi $pi2 $torad $todeg $h $h_bar $hbar $c $e $me $mp $mn
				$u0 $e0 $e2_4pie0 $a0 $kB $NA $basee $R $F $g 
				$HartreeToeV $RyToeV $KToeV $eVToK $JToeV $eVToJ
				$G $DayToSecond $SecondToDay $AstronomicalUnit $AU
				mod log10
				tan dtan dcos dsin acos asin atan dacos dasin datan
				tanh cosh sinh torad todeg 
				XRDDtoQ CircleArea 
				Gaussian Lorentzian erf erfc Gamma LegendrePolynomial
				eVTonm nmToeV
				FuncLinear FuncMinus FuncLog10 FuncLog10Abs FuncLn FuncLnAbs
				FD FDT FDE FermiDiracDistribution PoweredFermiDirac 
				FermiDiracIntegralByNumericalIntegration
				FermiDiracIntegral
				);

use strict;
use Sci::Algorism;

#BEGIN { $Exporter::Verbose = 1 }

our $pi          = 3.14159265358979323846; # 3.1415926535;
our $pi2         = $pi + $pi;
our $torad       = 0.01745329251944; # rad/deg";
our $todeg       = 57.29577951472;   # deg/rad";
our $h           = 6.6260755e-34;    # Js";
our $h_bar       = 1.05459e-34;      # "Js";
our $hbar        = $h_bar;
our $c           = 2.99792458e8;     # m/s";
our $e           = 1.60218e-19;      # C";
our $me          = 9.1093897e-31;    # kg";
our $mp          = 1.6726231e-27;    # kg";
our $mn          = 1.67495e-27;      # kg";
our $u0          = 4*3.14*1e-7;      # . "Ns<sup>2</sup>C<sup>-2</sup>";
our $e0          = 8.854418782e-12;  # C<sup>2</sup>N<sup>-1</sup>m<sup>-2</sup>";
our $e2_4pie0    = 2.30711e-28;      # Nm<sup>2</sup>";
our $a0          = 5.29177e-11;      # m";
our $kB          = 1.380658e-23;     # JK<sup>-1</sup>";
our $NA          = 6.0221367e23;     # mol<sup>-1</sup>";
our $basee       = 2.71828183;
our $R           = 8.31451;          # J/K/mol";
our $F           = 96485.3;          # C/mol";
our $g           = 9.81;             # m/s2";
our $HartreeToeV = 27.2116; #27.1481;          # eV";
#     $me * $e*$e*$e / (4.0*$pi*$e0)^2 / $h_bar/$h_bar;
our $HTV         = $HartreeToeV;
$HTV =~ s/\s+.*$//;
our $RyToeV      = $HTV / 2.0;
our $KToeV       = $kB / $e;
our $eVToK       = $e / $kB;
our $JToeV       = 1.0 / $e;
our $eVToJ       = $e;
our $Debye            = 3.33564e-30;   # Cm, for dipole (D)

our $G                = 6.67259e-11;   #Nm2/kg2 
our $DayToSecond      = 60 * 60 * 24;  #s
our $SecondToDay      = 1.0 / $DayToSecond;
our $AstronomicalUnit = 1.49597870e11; #m
our $AU               = $AstronomicalUnit;

sub a0    { return $a0; }
sub u0    { return $u0; }
sub e0    { return $e0; }
sub pi    { return $pi; }
sub ub    { return "9.27402e-24 J/T";  }
sub un    { return "5.050787e-27 J/T"; }
sub NA    { return $NA; }
sub R     { return $R;  }
sub kB    { return $kB; }
sub F     { return $F;  }
sub g     { return $g;  }
sub c     { return $c;  }
sub me    { return $me; }
sub mn    { return $mn; }
sub mp    { return $mp; }
sub e     { return $e;  }
sub h     { return $h;  }
sub h_bar { return $h_bar;  }
sub hbar  { return $h_bar;  }
sub G     { return $G; }
sub Debye { return $Debye; }

sub todeg { return $todeg; }
sub torad { return $torad; }
sub HartreeToeV { return $HartreeToeV; }
sub RyToeV      { return $RyToeV; }
sub eVToHartree { my $a = $HartreeToeV; $a =~ s/\s+.*$//; return 1.0/$a; }
sub eVToRy      { return 1.0/$RyToeV; }
sub JToeV       { return $JToeV; }

sub DayToSecond      { return $DayToSecond; }
sub SecondToDay      { return $SecondToDay; }
sub AstronomicalUnit { return $AstronomicalUnit; }
sub eVTonm           { return ($h / ($_[0]*$e)) * $c * 1e9; }
sub nmToeV           { 
	my ($wl) = @_;
	return 0.0 if($wl == 0.0);
	return ($h / ($wl*$e)) * $c * 1e9; 
}

sub mod
{
	my ($val, $max) = @_;
	if($val > $max) {
		$val -= int($val / $max) * $max;
	}
	elsif($val < 0.0) {
		$val -= (-1 + int($val / $max)) * $max;
	}
	return $val;
}
sub tan    { return sin($_[0]) / cos($_[0]); }
sub dtan   { return tan($_[0]*$torad);  }
sub dcos   { return cos($_[0]*$torad);  }
sub dsin   { return sin($_[0]*$torad);  }
sub asin{
	my $sinQ = $_[0];
	my $cosQ = sqrt(1.0 - $sinQ*$sinQ);
	return atan2($sinQ, $cosQ);
}
sub acos {
	my ($cosa) = @_;
	return $pi/2.0 if(abs($cosa) < 0.000001);
	return 0.0 if(abs($cosa - 1.0) < 0.000001);
	return $pi if(abs($cosa + 1.0) < 0.000001);
	my $sina = 1.0 - $cosa*$cosa;
	return 0.0 if($sina <= 0.0);
	$sina = sqrt($sina);
	return atan2($sina, $cosa);
}
sub atan       { return atan2($_[0], 1.0);    }
sub dacos      { return acos($_[0])*$todeg;   }
sub dasin  { return asin($_[0])*$todeg; }
sub datan      { return atan($_[0])*$todeg;   }

sub cosh   { return (exp($_[0])+exp(-$_[0]))/2.0 }
sub sinh   { return (exp($_[0])-exp(-$_[0]))/2.0 }
sub tanh   { return sinh($_[0]) / cosh($_[0]);   }

sub log10      { return log($_[0]) * 0.4342944819; }

sub CircleArea { return $_[0]*$_[0]*$pi; }
sub XRDDtoQ    { return asin(1.5405/2/$_[0])*$todeg/2.0; }

sub FuncLinear     { return $_[0]; }
sub FuncMinus      { return -$_[0]; }
sub FuncLog10      { return log10($_[0]) if($_[0] > 0);         return -5; }
sub FuncLog10Abs   { return log10(abs($_[0])) if($_[0] != 0.0); return -5; }
sub FuncLn         { return log($_[0]) if($_[0] > 0);           return -5; }
sub FuncLnAbs      { return log(abs($_[0])) if($_[0] != 0);     return -5; }
sub FuncInverse    { return 1 / $_[0] if($_[0] != 0);           return 1000; }
sub FuncPower10    { return 10**($_[0]); }
sub FuncPower10Abs { return 10**(abs($_[0])); }
sub FuncExp        { return exp($_[0]); }
sub FuncExpAbs     { return exp(abs($_[0])); }

sub Factorial {
	my ($n) = @_;
	return 1 if($n <= 1);
	my $y = 1;
	for(my $i = 2 ; $i <= $n ; $i++) {
		$y *= $i;
	}
	return $y;
}

sub HPG
{
	my ($n, $x, $x0, $b) = @_;

	my $s = ($x - $x0) / $b;
	my $G = &Gaussian($s, 0.0, 1.0);
	my $H = &HermitePolynomial($n, $s);
	my $k = sqrt(2.0**$n * sqrt($pi) * Sci::Factorial($n));
	return $k * $H * $G;
}

sub HermitePolynomial
{
	my ($n, $x) = @_;

	return 1.0    if($n == 0);
	return 2.0*$x if($n == 1);
	return 2.0 * ($x * &HermitePolynomial($n-1, $x) - ($n-1) * &HermitePolynomial($n-2, $x));
}

sub HermitePolynomial2
{
	my ($n, $x) = @_;

	my $y = 0.0;
	for(my $m = 0 ; $m <= $n / 2 ; $m++) {
		my $n1  = Sci::Factorial($n);
		my $m1  = Sci::Factorial($m);
		my $mn1 = Sci::Factorial($n - 2*$m);
		$y += $n1 * (-1)**$m * (2.0*$x)**($n-2*$m) / $m1 / $mn1;
	}
	return $y;
}

sub LegendrePolynomial
{
	my ($n, $x) = @_;
	
	if($n == 0) {
		return 1.0;
	}
	elsif($n == 1) {
		return $x;
	}
	my $n1 = $n - 1;
	return 1.0 / $n * ((2.0*$n1 + 1.0) * &LegendrePolynomial($n-1, $x) - $n1 * &LegendrePolynomial($n-2, $x));
}

#n: Legendre多項式の次数
sub LegendrePolynomial_prev
{
	my ($x, $n) = @_;
	
	return undef if($n < 0);
	return 1.0 if($n == 0);
	return $x if($n == 1);
	my $wx = $x;
	my $wn = $n;
	my $p0 = 1.0;
	my $p1 = $wx;
	my $p2;
	do {
		$wn -= 1.0;
		my $w = $wn / ($wn + 1.0);
		$p2 = $wx * $p1 + $w * ($wx * $p1 - $p0);
		$p0 = $p1;
		$p1 = $p2;
	} while ($wn >= 2.0);
	return $p2;
}

sub LegendrePolynomial2
{
	my ($l, $m, $x) = @_;
#print "m=$m, l=$l\n";
	if($m < 0) {
		my $lmm1 = Sci::Factorial($l+$m);
		my $lmp1 = Sci::Factorial($l-$m);
		my $c = (-1)**(-$m) * $lmm1 / $lmp1;
		return $c * &LegendrePolynomial2($l, -$m, $x);
	}

	if($l < $m or $l < 0) {
		return 0.0;
	}
	elsif($l == $m) {
		my $m21 = Sci::Factorial(2*$m - 1);
		return $m21 * (1.0 - $x*$x)**($m/2.0);
	}

	my $c1 = (2.0*$l - 1.0) / ($l - $m);
	my $c2 = ($l + $m - 1.0) / ($l - $m);
	return $c1 * $x * &LegendrePolynomial2($l-1, $m, $x) - $c2 * &LegendrePolynomial2($l-2, $m, $x);
}

sub LaguerrePolynomial
{
	my ($n, $x) = @_;
	
	if($n == 0) {
		return 1.0;
	}
	elsif($n == 1) {
		return -$x + 1.0;
	}
	my $n1 = $n - 1;
	return (2.0*$n1 + 1.0 - $x) * &LaguerrePolynomial($n-1, $x) - $n1*$n1 * &LaguerrePolynomial($n-2, $x);
}

sub LaguerrePolynomial2
{
	my ($n, $k, $x) = @_;

	my $f = 0.0;
	for(my $m = 0 ; $m <= $n-$k ; $m++) {
		my $n1   = Sci::Factorial($n);
		my $m1   = Sci::Factorial($m);
		my $mk1  = Sci::Factorial($m+$k);
		my $nmk1 = Sci::Factorial($n-$m-$k);
		$f += (-1)**($m+$k) * $n1*$n1 / ($m1 * $mk1 * $nmk1) * $x**$m;
	}
	return $f;
}

sub MethfesselPaxton
{
	my ($N, $x, $w) = @_;
	$w = 1.0 if(!defined $w);
	my $s = $x / $w;

	my $y = 0.0;
	for(my $n = 0 ; $n <= $N ; $n++) {
		my $An = (-1)**$n / Sci::Factorial($n) / 4.0**$n / sqrt($pi);
		my $y0 = $An / $w * Sci::HermitePolynomial(2*$n, $s) * exp(-($s*$s));
		$y += $y0;
#print "N=$N, x=$x: n=$n: An=$An, y0=$y0, y=$y\n";
	}
	return $y;
}

sub Round
{
	my ($val, $digit) = @_;
$val = 0.0 if(!defined $val);
$digit = 0.0 if(!defined $digit);

	$digit = int($digit+0.1);
	$digit = 30 if($digit > 30);
	$digit = 4 if($digit <= 0);

	my $k = 1;
	for(my $i = 0 ; $i < $digit ; $i++) {
		$k *= 10;
	}
	return int($val*$k + 0.1) / $k;
}

sub Polynomial
{
	my ($x, $m, @c) = @_;
	my $v = $c[0];
	for(my $i = 1 ; $i < $m+1 ; $i++) {
		$v += $c[$i] * ($x ** $i)
	}
	return $v;
}

sub Gaussian
{
	my ($x, $x0, $whalf) = @_;
#A = 1/whalf * sqrt(ln2 / pi)
	my $A = 0.469718639 / $whalf;
#a = whalf / sqrt(ln2)
	my $a = $whalf / 0.832554611;
	my $X = ($x - $x0) / $a;
	return $A * exp(-$X*$X);
}

sub Lorentzian
{
	my ($x, $x0, $whalf) = @_;
#A = 1/whalf/pi
	my $A = 1.0 / $whalf / $pi;
	my $X = ($x - $x0) / $whalf;
	return $A / (1.0 + $X*$X);
}

sub GaussLorentz
{
	my ($x, $x0, $C0, $WL, $GFraction, $GWRatio) = @_;
	my $y = $C0 * $GFraction * Gaussian($x, $x0, $WL*$GWRatio);
	$y += $C0 * (1.0 - $GFraction) * Lorentzian($x, $x0, $WL);
	return $y;
}

sub GaussLorentzArea 
{
	my ($C0, $WL, $GFraction, $GWRatio) = @_;
	my $y = $C0 * $GFraction * sqrt($pi) * $WL*$GWRatio / 0.832554611; #Gaussian($x, $x0, $WL*$GWRatio);
	$y += $C0 * (1.0 - $GFraction) * $pi * $WL; #Lorentzian($x, $x0, $WL);
	return $y;
}

my @erfc_a = (1.0, 0.0705230784, 0.0422820123, 0.92705272e-2, 0.1520143e-3, 0.2765672e-3);
sub erfc
{
	my ($x) = @_;
	return 1.0 - erf($x);

	my $wx = $x;
	return undef if($wx < 0.0);
	return 0.0 if($wx ==  0.0);
	return 1.0 if($wx >= 10.0);
	my $w = 0.430638e-4;
	for(my $i = 5 ; $i >= 0 ; $i--) {
		$w = $w * $wx + $erfc_a[$i];
	}
	for(my $i = 0 ; $i < 4 ; $i++) {
		$w *= $w;
	}
	return 1.0 / $w;
}

sub erf
{
	my ($x) = @_;

	my $wx = $x;
	return undef if($wx < 0.0);
	return 0.0 if($wx ==  0.0);
	return 1.0 if($wx >= 10.0);
	my $w = 0.430638e-4;
	for(my $i = 5 ; $i >= 0 ; $i--) {
		$w = $w * $wx + $erfc_a[$i];
	}
	for(my $i = 0 ; $i < 4 ; $i++) {
		$w *= $w;
	}
	return 1.0 - 1.0 / $w;
}

my @Gamma_a = (1.0,         -0.577191652, 0.988205891, -0.897056937, 
			   0.918206857, -0.756704078, 0.482199394, -0.193527818);
sub Gamma
{
	my ($x) = @_;

	my $wx = $x;
	return undef if($wx <= 0.0 or $wx >= 34.0);
	my $w1 = 1.0;
	if($wx > 1.0) {
		do {
			$wx -= 1.0;
			$w1 *= $wx;
		} while($wx > 1.0);
	}
	my $w = 0.035868343;
	for(my $i = 7 ; $i >= 0 ; $i--) {
		$w = $w * $wx + $Gamma_a[$i];
	}
	return $w1 * $w / $wx;
}

sub ElectroneVTonm
{
	my ($ElectronEnergy) = @_;
	return sqrt($h * $h / 2.0 / $me / $ElectronEnergy / $e) * 1.0e9; # in nm
}

#============================================================
# Vector
#============================================================

sub NormalizeVector
{
	my ($pv1, $l) = @_;
	$l = 1.0 if(!defined $l);

	my $L2 = &InnerProduct($pv1, $pv1);
	return undef if($L2 < 1.0e-15);

	my $L  = sqrt($L2);
	my $n = @$pv1;
	my @NormalizeV;
	for(my $i = 0 ; $i < $n ; $i++) {
		$NormalizeV[$i] = $pv1->[$i] * $l / $L;
	}
	return \@NormalizeV;
}

sub VectorProduct
{
	my ($pv1, $pv2, $k) = @_;
	$k = 1.0 if(!defined $k);
	my $n = @$pv1;
	return (($pv1->[1]*$pv2->[2] - $pv1->[2]*$pv2->[1])*$k,
			($pv1->[2]*$pv2->[0] - $pv1->[0]*$pv2->[2])*$k,
			($pv1->[0]*$pv2->[1] - $pv1->[1]*$pv2->[0])*$k);
}

sub ScalarProduct
{
	my ($pv1, $pv2) = @_;
	return &InnerProduct($pv1, $pv2);
}

sub InnerProduct
{
	my ($pv1, $pv2) = @_;
	my $n = @$pv1;
	my $val = 0.0;
	for(my $i = 0 ; $i < $n ; $i++) {
		$val += $pv1->[$i] * $pv2->[$i];
	}
	return $val;
}

sub CalVectorAngle
{
	my ($pv1, $pv2) = @_;
	my $l1 = sqrt(InnerProduct($pv1, $pv1));
	return undef if(abs($l1) < 1.0e-15);
	my $l2 = sqrt(InnerProduct($pv2, $pv2));
	return undef if(abs($l2) < 1.0e-15);
	my $cosQ = InnerProduct($pv1, $pv2) / $l1 / $l2;
	return $todeg * acos($cosQ);
}


#============================================================
# 熱電変換工学−基礎と応用−、編集委員長：阪田亮、REALIZE INC.、2001
#============================================================
sub DiffFermiDiracDistribution
{
	my ($w, $x0, $x) = @_;
	my $s = ($x - $x0) / $w;
# f = 1.0 / (1.0 + exp(s))
	return exp($s) / $w / (1.0 + exp($s)) / (1.0 + exp($s));
}


sub FermiDiracDistribution
{
	my ($x, $x0) = @_;
	return &FD($x, $x0);
}

# 電子のFermi-Dirac分布
sub FEe
{
	my ($E, $EF, $T) = @_;

	if($T == 0.0) {
		return 0.0 if($E > $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($E-$EF) / $T / $KToeV) );
}

# 正孔のFermi-Dirac分布
sub FEh
{
	my ($E, $EF, $T) = @_;

	if($T == 0.0) {
		return 0.0 if($E < $EF);
		return 0.5 if($E == $EF);
		return 1.0;
	}
	return 1.0 / ( 1.0 + exp(($EF-$E) / $T / $KToeV) );
}

sub FDTeV
{
	my ($T, $EF, $E) = @_;
	if($T <= 0.0) {
		return 1.0 if($E < $EF);
		return 0.0;
	}

	my $kBT = $kB * $T * $JToeV;
	my $ee =  ($E - $EF) / $kBT;
	return 1.0 / (1.0 + exp($ee));
}

sub FDT
{
	my ($T, $EF, $E) = @_;
	if($T <= 0.0) {
		return 1.0 if($E < $EF);
		return 0.0;
	}

	my $kBT = $kB * $T;
	my $ee =  ($E - $EF) / $kBT;
	return 1.0 / (1.0 + exp($ee));
}

sub FDE
{
	my ($EW, $EF, $E) = @_;
	if($EW <= 0.0) {
		return 1.0 if($E < $EF);
		return 0.0;
	}

	my $ee =  ($E - $EF) / $EW;
	return 1.0 / (1.0 + exp($ee));
}

sub FD
{
	my ($x, $x0) = @_;
	return 1.0 / (1.0 + exp($x-$x0));
}

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);
	my $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;
	my $sb = &PoweredFermiDirac($x, $xi, $r);
	return ($sa + $sb + 2.0*$s2 + 4.0*$s4) * $dx / 3.0;
}

my $KBBRadiation = 2.0 * $h * $c * $c;
my $kexpBBRadiation = $h * $c / $kB;
sub BlackbodyRadiationBynm
{
	my ($l, $T) = @_; # nm

	my $lambda  = $l * 1.0e-9;
	my $lambda2 = $lambda * $lambda;
	my $ex = exp($kexpBBRadiation / $lambda / $T);
	return $KBBRadiation / $lambda2 / $lambda2 / $lambda / ($ex - 1.0);
}

#============================================================
#　Sci::Gammaと同じ関数。使わない
#============================================================
# 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;
}

1;