.
#======================================================================
# L I B F U N S . P L
# doc: Wed Mar 24 11:49:13 1999
# dlm: Fri Sep 7 11:11:09 2012
# (c) 1999 A.M. Thurnherr
# uE-Info: 286 38 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Mar 24, 1999: - copied from the c-version of NR
# Mar 26, 1999: - added stuff for better [./fit]
# Sep 18, 1999: - argument typechecking
# Oct 04, 1999: - added gauss(), normal()
# Jan 25, 2001: - added f(), sgn()
# Apr 16, 2010: - added sinc()
# Sep 7, 2012: - added acosh()
require "$ANTS/libvec.pl"; # rad()
#----------------------------------------------------------------------
sub gauss(@)
{
my($x,$peak,$mean,$efs) = &antsFunUsage(4,"ffff","x, peak, mean, e-folding scale",@_);
return $peak * exp( -(($x-$mean) / $efs)**2);
}
sub normal(@)
{
my($x,$area,$mean,$sigma) = &antsFunUsage(4,"ffff","x, area, mean, stddev",@_);
my($pi) = 3.14159265358979;
return $area/(sqrt(2*$pi)*$sigma) * exp(-((($x-$mean) / $sigma)**2)/2);
}
#----------------------------------------------------------------------
# &f(lat) calculate coriolis param
#----------------------------------------------------------------------
sub f(@)
{
my($lat) = &antsFunUsage(1,"f","lat",@_);
my($Omega) = 7.292e-5; # Gill (1982)
return 2 * $Omega * sin(rad($lat));
}
#----------------------------------------------------------------------
# &sgn(v) return -1/0/+1
#----------------------------------------------------------------------
sub sgn(@)
{
my($val) = &antsFunUsage(1,"f","val",@_);
return 0 if ($val == 0);
return ($val < 0) ? -1 : 1;
}
#======================================================================
# rest of library cooked up from the diverse special function routines of NR
# Chapter 6. No attempt to clean up the code has been made.
#----------------------------------------------------------------------
# 6.1 Gamma Function et al
#----------------------------------------------------------------------
sub gammln(@)
{
my($xx) = &antsFunUsage(1,"f","xx",@_);
my($x,$y,$tmp,$ser);
my(@cof) = (76.18009172947146, -86.50532032941677,
24.01409824083091, -1.231739572450155,
0.1208650973866179e-2, -0.5395239384953e-5);
my($j);
$x = $xx;
$y = $x;
$tmp = $x + 5.5;
$tmp -= ($x+0.5) * log($tmp);
$ser = 1.000000000190015;
for ($j=0; $j<=5; $j++) {
$ser += $cof[$j] / ++$y;
}
return -$tmp + log(2.5066282746310005*$ser/$x);
}
#----------------------------------------------------------------------
# 6.2. Incomplete Gamma Function, Error Function et al
#----------------------------------------------------------------------
{ my($ITMAX)=100; my($EPS)=3.0e-7; # static vars
sub gser(@)
{
my($a,$x,$glnR) = &antsFunUsage(-2,"ff","a,x[,ref to gln]",@_);
my($gln);
my($n);
my($sum,$del,$ap);
$gln = &gammln($a);
$$glnR = $gln if (defined($glnR));
return 0 if ($x == 0);
croak("$0 (libspecfuns.pl): x<0 ($x) in &gser()\n")
if ($x < 0);
$ap = $a;
$sum = 1 / $a;
$del = $sum;
for ($n=1; $n<=$ITMAX; $n++) {
++$ap;
$del *= $x/$ap;
$sum += $del;
return $sum * exp(-$x+$a*log($x)-$gln)
if (abs($del) < abs($sum)*$EPS);
}
croak("$0 (libspecfuns.pl): a ($a) too large, " .
"ITMAX ($ITMAX) too small in &gser()\n");
}
} # end of static scope
{ my($ITMAX)=100; my($EPS)=3.0e-7; my($FPMIN)=1.0e-30; # static
sub gcf(@)
{
my($a,$x,$glnR) = &antsFunUsage(-2,"ff","a,x[,ref to gln]",@_);
my($gln);
my($i);
my($an,$b,$c,$d,$del,$h);
$gln = &gammln($a);
$$glnR = $gln if (defined($glnR));
$b = $x + 1 - $a;
croak("$0 (libspecfuns.pl): illegal params (a = x + 1) in &gcf()\n")
unless ($b);
$c = 1 / $FPMIN;
$d = 1 / $b;
$h = $d;
for ($i=1; $i<=$ITMAX; $i++) {
$an = -$i * ($i - $a);
$b += 2.0;
$d = $an * $d + $b;
$d = $FPMIN if (abs($d) < $FPMIN);
$c = $b + $an/$c;
$c = $FPMIN if (abs($c) < $FPMIN);
$d = 1 / $d;
$del= $d * $c;
$h *= $del;
last if (abs($del-1) < $EPS);
}
croak("$0 (libspecfuns.pl): a ($a) too large," .
" ITMAX ($ITMAX) too small in &gcf()\n")
if ($i > $ITMAX);
return exp(-$x + $a*log($x) - $gln) * $h;
}
} # end of static scope
sub gammq(@)
{
my($a,$x) = &antsFunUsage(2,"ff","a,x",@_);
croak("$0 (libspecfuns.pl): Invalid arguments in &gammq()\n")
if ($x < 0 || $a <= 0);
return ($x < ($a+1)) ?
1 - &gser($a,$x) :
&gcf($a,$x);
}
#----------------------------------------------------------------------
sub erfcc(@)
{
my($x) = &antsFunUsage(1,"f","x",@_);
my($t,$z,$ans);
$z = abs($x);
$t = 1/(1+0.5*$z);
$ans = $t*exp(-$z*$z-1.26551223+$t*(1.00002368+$t*(0.37409196+$t*(0.09678418+
$t*(-0.18628806+$t*(0.27886807+$t*(-1.13520398+$t*(1.48851587+
$t*(-0.82215223+$t*0.17087277)))))))));
return $x >= 0 ? $ans : 2.0-$ans;
}
{ my($warned) = 0; # static
sub erf(@)
{
my($x) = &antsFunUsage(1,"f","x",@_);
&antsInfo("(libspecfuns.pl) WARNING: using approximate erf()"),$warned=1
unless ($warned);
return 1-&erfcc($x);
}
}
#----------------------------------------------------------------------
# 6.3. Incomplete Beta Function et al
#----------------------------------------------------------------------
sub betai(@)
{
my($a,$b,$x) = &antsFunUsage(3,"fff","a,b,x",@_);
my($bt);
croak("$0 (liberrf.pl): x (=$x) out of range in betai()\n")
if ($x < 0 || $x > 1);
if ($x == 0 || $x == 1) {
$bt = 0;
} else {
$bt = exp(gammln($a+$b)-gammln($a)-gammln($b)+$a*log($x)+$b*log(1-$x));
}
if ($x < ($a+1)/($a+$b+2)) {
return $bt * betacf($a,$b,$x) / $a;
} else {
return 1 - $bt*betacf($b,$a,1-$x) / $b;
}
}
#----------------------------------------------------------------------
{ # static scope
my($MAXIT) = 100;
my($EPS) = 3.0e-7;
my($FPMIN) = 1.0e-30;
sub betacf(@)
{
my($a,$b,$x) = &antsFunUsage(3,"fff","a,b,x",@_);
my($m,$m2);
my($aa,$c,$d,$del,$h,$qab,$qam,$qap);
$qab = $a + $b;
$qap = $a + 1;
$qam = $a - 1;
$c = 1;
$d = 1 - $qab*$x/$qap;
$d = $FPMIN if (abs($d) < $FPMIN);
$d = 1 / $d;
$h = $d;
for ($m=1; $m<=$MAXIT; $m++) {
$m2 = 2 * $m;
$aa = $m*($b-$m)*$x / (($qam+$m2)*($a+$m2));
$d = 1 + $aa*$d;
$d = $FPMIN if (abs($d) < $FPMIN);
$c = 1 + $aa/$c;
$c = $FPMIN if (abs($c) < $FPMIN);
$d = 1 / $d;
$h *= $d * $c;
$aa = -($a+$m)*($qab+$m)*$x / (($a+$m2)*($qap+$m2));
$d = 1 + $aa*$d;
$d = $FPMIN if (abs($d) < $FPMIN);
$c = 1 + $aa/$c;
$c = $FPMIN if (abs($c) < $FPMIN);
$d = 1 / $d;
$del= $d * $c;
$h *= $del;
last if (abs($del-1) < $EPS);
}
croak("$0 (liberrf.pl): a or b too big, or MAXIT too small in betacf")
if ($m > $MAXIT);
return $h;
}
} # end of static scope
#----------------------------------------------------------------------
# normalized cardinal sine as used, e.g., in JAOT/polzin02
#----------------------------------------------------------------------
sub sinc($)
{
my($piX) = 3.14159265358979 * $_[0];
return $piX==0 ? 1 : sin($piX)/$piX;
}
#----------------------------------------------------------------------
# inverse hyperbolic cosine; mathworld
# - requires argument >= 1
#----------------------------------------------------------------------
sub acosh($)
{
return log($_[0] + sqrt($_[0]**2-1));
}
#----------------------------------------------------------------------
1;