0
|
1 |
#======================================================================
|
|
2 |
# L I B F U N S . P L
|
|
3 |
# doc: Wed Mar 24 11:49:13 1999
|
36
|
4 |
# dlm: Fri May 11 11:40:05 2018
|
0
|
5 |
# (c) 1999 A.M. Thurnherr
|
36
|
6 |
# uE-Info: 31 77 NIL 0 0 70 2 2 4 NIL ofnI
|
0
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# HISTORY:
|
|
10 |
# Mar 24, 1999: - copied from the c-version of NR
|
|
11 |
# Mar 26, 1999: - added stuff for better [./fit]
|
|
12 |
# Sep 18, 1999: - argument typechecking
|
|
13 |
# Oct 04, 1999: - added gauss(), normal()
|
|
14 |
# Jan 25, 2001: - added f(), sgn()
|
|
15 |
# Apr 16, 2010: - added sinc()
|
3
|
16 |
# Sep 7, 2012: - added acosh()
|
20
|
17 |
# Jun 4, 2015: - added gaussRand()
|
|
18 |
# - made normal() more efficient
|
36
|
19 |
# May 11, 2018: - added Nsq()
|
0
|
20 |
|
|
21 |
require "$ANTS/libvec.pl"; # rad()
|
|
22 |
|
|
23 |
#----------------------------------------------------------------------
|
36
|
24 |
# Buoyancy-Freuquency Squared
|
|
25 |
# - based on signed buoyancy frequency => propagate sign
|
|
26 |
#----------------------------------------------------------------------
|
|
27 |
|
|
28 |
{ my(@fc);
|
|
29 |
sub Nsq(@)
|
|
30 |
{
|
|
31 |
my($N) = &antsFunUsage(1,'.','[(signed) buoyancy frequency]',\@fc,'N',@_);
|
|
32 |
return ($N < 0) ? -($N**2) : $N**2;
|
|
33 |
}
|
|
34 |
}
|
|
35 |
|
|
36 |
#----------------------------------------------------------------------
|
20
|
37 |
# gaussians/normal distribution
|
|
38 |
#----------------------------------------------------------------------
|
0
|
39 |
|
|
40 |
sub gauss(@)
|
|
41 |
{
|
|
42 |
my($x,$peak,$mean,$efs) = &antsFunUsage(4,"ffff","x, peak, mean, e-folding scale",@_);
|
|
43 |
return $peak * exp( -(($x-$mean) / $efs)**2);
|
|
44 |
}
|
|
45 |
|
|
46 |
sub normal(@)
|
|
47 |
{
|
|
48 |
my($x,$area,$mean,$sigma) = &antsFunUsage(4,"ffff","x, area, mean, stddev",@_);
|
20
|
49 |
my($sqrt2pi) = 2.506628274631;
|
|
50 |
return $area/($sqrt2pi*$sigma) * exp(-((($x-$mean) / $sigma)**2)/2);
|
0
|
51 |
}
|
|
52 |
|
|
53 |
#----------------------------------------------------------------------
|
|
54 |
# &f(lat) calculate coriolis param
|
|
55 |
#----------------------------------------------------------------------
|
|
56 |
|
|
57 |
sub f(@)
|
|
58 |
{
|
|
59 |
my($lat) = &antsFunUsage(1,"f","lat",@_);
|
|
60 |
my($Omega) = 7.292e-5; # Gill (1982)
|
|
61 |
|
|
62 |
return 2 * $Omega * sin(rad($lat));
|
|
63 |
}
|
|
64 |
|
|
65 |
#----------------------------------------------------------------------
|
|
66 |
# &sgn(v) return -1/0/+1
|
|
67 |
#----------------------------------------------------------------------
|
|
68 |
|
|
69 |
sub sgn(@)
|
|
70 |
{
|
|
71 |
my($val) = &antsFunUsage(1,"f","val",@_);
|
|
72 |
return 0 if ($val == 0);
|
|
73 |
return ($val < 0) ? -1 : 1;
|
|
74 |
}
|
|
75 |
|
|
76 |
#======================================================================
|
|
77 |
|
|
78 |
# rest of library cooked up from the diverse special function routines of NR
|
|
79 |
# Chapter 6. No attempt to clean up the code has been made.
|
|
80 |
|
|
81 |
#----------------------------------------------------------------------
|
|
82 |
# 6.1 Gamma Function et al
|
|
83 |
#----------------------------------------------------------------------
|
|
84 |
|
|
85 |
sub gammln(@)
|
|
86 |
{
|
|
87 |
my($xx) = &antsFunUsage(1,"f","xx",@_);
|
|
88 |
my($x,$y,$tmp,$ser);
|
|
89 |
my(@cof) = (76.18009172947146, -86.50532032941677,
|
|
90 |
24.01409824083091, -1.231739572450155,
|
|
91 |
0.1208650973866179e-2, -0.5395239384953e-5);
|
|
92 |
my($j);
|
|
93 |
|
|
94 |
$x = $xx;
|
|
95 |
$y = $x;
|
|
96 |
$tmp = $x + 5.5;
|
|
97 |
$tmp -= ($x+0.5) * log($tmp);
|
|
98 |
$ser = 1.000000000190015;
|
|
99 |
for ($j=0; $j<=5; $j++) {
|
|
100 |
$ser += $cof[$j] / ++$y;
|
|
101 |
}
|
|
102 |
return -$tmp + log(2.5066282746310005*$ser/$x);
|
|
103 |
}
|
|
104 |
|
|
105 |
#----------------------------------------------------------------------
|
|
106 |
# 6.2. Incomplete Gamma Function, Error Function et al
|
|
107 |
#----------------------------------------------------------------------
|
|
108 |
|
|
109 |
{ my($ITMAX)=100; my($EPS)=3.0e-7; # static vars
|
|
110 |
|
|
111 |
sub gser(@)
|
|
112 |
{
|
|
113 |
my($a,$x,$glnR) = &antsFunUsage(-2,"ff","a,x[,ref to gln]",@_);
|
|
114 |
my($gln);
|
|
115 |
my($n);
|
|
116 |
my($sum,$del,$ap);
|
|
117 |
|
|
118 |
$gln = &gammln($a);
|
|
119 |
$$glnR = $gln if (defined($glnR));
|
|
120 |
|
|
121 |
return 0 if ($x == 0);
|
|
122 |
croak("$0 (libspecfuns.pl): x<0 ($x) in &gser()\n")
|
|
123 |
if ($x < 0);
|
|
124 |
|
|
125 |
$ap = $a;
|
|
126 |
$sum = 1 / $a;
|
|
127 |
$del = $sum;
|
|
128 |
for ($n=1; $n<=$ITMAX; $n++) {
|
|
129 |
++$ap;
|
|
130 |
$del *= $x/$ap;
|
|
131 |
$sum += $del;
|
|
132 |
return $sum * exp(-$x+$a*log($x)-$gln)
|
|
133 |
if (abs($del) < abs($sum)*$EPS);
|
|
134 |
}
|
|
135 |
croak("$0 (libspecfuns.pl): a ($a) too large, " .
|
|
136 |
"ITMAX ($ITMAX) too small in &gser()\n");
|
|
137 |
}
|
|
138 |
|
|
139 |
} # end of static scope
|
|
140 |
|
|
141 |
{ my($ITMAX)=100; my($EPS)=3.0e-7; my($FPMIN)=1.0e-30; # static
|
|
142 |
|
|
143 |
sub gcf(@)
|
|
144 |
{
|
|
145 |
my($a,$x,$glnR) = &antsFunUsage(-2,"ff","a,x[,ref to gln]",@_);
|
|
146 |
my($gln);
|
|
147 |
my($i);
|
|
148 |
my($an,$b,$c,$d,$del,$h);
|
|
149 |
|
|
150 |
$gln = &gammln($a);
|
|
151 |
$$glnR = $gln if (defined($glnR));
|
|
152 |
|
|
153 |
$b = $x + 1 - $a;
|
|
154 |
croak("$0 (libspecfuns.pl): illegal params (a = x + 1) in &gcf()\n")
|
|
155 |
unless ($b);
|
|
156 |
$c = 1 / $FPMIN;
|
|
157 |
$d = 1 / $b;
|
|
158 |
$h = $d;
|
|
159 |
for ($i=1; $i<=$ITMAX; $i++) {
|
|
160 |
$an = -$i * ($i - $a);
|
|
161 |
$b += 2.0;
|
|
162 |
$d = $an * $d + $b;
|
|
163 |
$d = $FPMIN if (abs($d) < $FPMIN);
|
|
164 |
$c = $b + $an/$c;
|
|
165 |
$c = $FPMIN if (abs($c) < $FPMIN);
|
|
166 |
$d = 1 / $d;
|
|
167 |
$del= $d * $c;
|
|
168 |
$h *= $del;
|
|
169 |
last if (abs($del-1) < $EPS);
|
|
170 |
}
|
|
171 |
croak("$0 (libspecfuns.pl): a ($a) too large," .
|
|
172 |
" ITMAX ($ITMAX) too small in &gcf()\n")
|
|
173 |
if ($i > $ITMAX);
|
|
174 |
return exp(-$x + $a*log($x) - $gln) * $h;
|
|
175 |
}
|
|
176 |
|
|
177 |
} # end of static scope
|
|
178 |
|
|
179 |
sub gammq(@)
|
|
180 |
{
|
|
181 |
my($a,$x) = &antsFunUsage(2,"ff","a,x",@_);
|
|
182 |
croak("$0 (libspecfuns.pl): Invalid arguments in &gammq()\n")
|
|
183 |
if ($x < 0 || $a <= 0);
|
|
184 |
return ($x < ($a+1)) ?
|
|
185 |
1 - &gser($a,$x) :
|
|
186 |
&gcf($a,$x);
|
|
187 |
}
|
|
188 |
|
|
189 |
#----------------------------------------------------------------------
|
|
190 |
|
|
191 |
sub erfcc(@)
|
|
192 |
{
|
|
193 |
my($x) = &antsFunUsage(1,"f","x",@_);
|
|
194 |
my($t,$z,$ans);
|
|
195 |
|
|
196 |
$z = abs($x);
|
|
197 |
$t = 1/(1+0.5*$z);
|
|
198 |
$ans = $t*exp(-$z*$z-1.26551223+$t*(1.00002368+$t*(0.37409196+$t*(0.09678418+
|
|
199 |
$t*(-0.18628806+$t*(0.27886807+$t*(-1.13520398+$t*(1.48851587+
|
|
200 |
$t*(-0.82215223+$t*0.17087277)))))))));
|
|
201 |
return $x >= 0 ? $ans : 2.0-$ans;
|
|
202 |
}
|
|
203 |
|
|
204 |
{ my($warned) = 0; # static
|
|
205 |
|
|
206 |
sub erf(@)
|
|
207 |
{
|
|
208 |
my($x) = &antsFunUsage(1,"f","x",@_);
|
|
209 |
&antsInfo("(libspecfuns.pl) WARNING: using approximate erf()"),$warned=1
|
|
210 |
unless ($warned);
|
|
211 |
return 1-&erfcc($x);
|
|
212 |
}
|
|
213 |
|
|
214 |
}
|
|
215 |
|
|
216 |
#----------------------------------------------------------------------
|
|
217 |
# 6.3. Incomplete Beta Function et al
|
|
218 |
#----------------------------------------------------------------------
|
|
219 |
|
|
220 |
sub betai(@)
|
|
221 |
{
|
|
222 |
my($a,$b,$x) = &antsFunUsage(3,"fff","a,b,x",@_);
|
|
223 |
my($bt);
|
|
224 |
|
|
225 |
croak("$0 (liberrf.pl): x (=$x) out of range in betai()\n")
|
|
226 |
if ($x < 0 || $x > 1);
|
|
227 |
if ($x == 0 || $x == 1) {
|
|
228 |
$bt = 0;
|
|
229 |
} else {
|
|
230 |
$bt = exp(gammln($a+$b)-gammln($a)-gammln($b)+$a*log($x)+$b*log(1-$x));
|
|
231 |
}
|
|
232 |
if ($x < ($a+1)/($a+$b+2)) {
|
|
233 |
return $bt * betacf($a,$b,$x) / $a;
|
|
234 |
} else {
|
|
235 |
return 1 - $bt*betacf($b,$a,1-$x) / $b;
|
|
236 |
}
|
|
237 |
}
|
|
238 |
|
|
239 |
#----------------------------------------------------------------------
|
|
240 |
|
|
241 |
{ # static scope
|
|
242 |
|
|
243 |
my($MAXIT) = 100;
|
|
244 |
my($EPS) = 3.0e-7;
|
|
245 |
my($FPMIN) = 1.0e-30;
|
|
246 |
|
|
247 |
sub betacf(@)
|
|
248 |
{
|
|
249 |
my($a,$b,$x) = &antsFunUsage(3,"fff","a,b,x",@_);
|
|
250 |
my($m,$m2);
|
|
251 |
my($aa,$c,$d,$del,$h,$qab,$qam,$qap);
|
|
252 |
|
|
253 |
$qab = $a + $b;
|
|
254 |
$qap = $a + 1;
|
|
255 |
$qam = $a - 1;
|
|
256 |
$c = 1;
|
|
257 |
$d = 1 - $qab*$x/$qap;
|
|
258 |
$d = $FPMIN if (abs($d) < $FPMIN);
|
|
259 |
$d = 1 / $d;
|
|
260 |
$h = $d;
|
|
261 |
for ($m=1; $m<=$MAXIT; $m++) {
|
|
262 |
$m2 = 2 * $m;
|
|
263 |
$aa = $m*($b-$m)*$x / (($qam+$m2)*($a+$m2));
|
|
264 |
$d = 1 + $aa*$d;
|
|
265 |
$d = $FPMIN if (abs($d) < $FPMIN);
|
|
266 |
$c = 1 + $aa/$c;
|
|
267 |
$c = $FPMIN if (abs($c) < $FPMIN);
|
|
268 |
$d = 1 / $d;
|
|
269 |
$h *= $d * $c;
|
|
270 |
$aa = -($a+$m)*($qab+$m)*$x / (($a+$m2)*($qap+$m2));
|
|
271 |
$d = 1 + $aa*$d;
|
|
272 |
$d = $FPMIN if (abs($d) < $FPMIN);
|
|
273 |
$c = 1 + $aa/$c;
|
|
274 |
$c = $FPMIN if (abs($c) < $FPMIN);
|
|
275 |
$d = 1 / $d;
|
|
276 |
$del= $d * $c;
|
|
277 |
$h *= $del;
|
|
278 |
last if (abs($del-1) < $EPS);
|
|
279 |
}
|
|
280 |
croak("$0 (liberrf.pl): a or b too big, or MAXIT too small in betacf")
|
|
281 |
if ($m > $MAXIT);
|
|
282 |
return $h;
|
|
283 |
}
|
|
284 |
|
|
285 |
} # end of static scope
|
|
286 |
|
|
287 |
#----------------------------------------------------------------------
|
|
288 |
# normalized cardinal sine as used, e.g., in JAOT/polzin02
|
|
289 |
#----------------------------------------------------------------------
|
|
290 |
|
|
291 |
sub sinc($)
|
|
292 |
{
|
|
293 |
my($piX) = 3.14159265358979 * $_[0];
|
|
294 |
return $piX==0 ? 1 : sin($piX)/$piX;
|
|
295 |
}
|
|
296 |
|
|
297 |
#----------------------------------------------------------------------
|
3
|
298 |
# inverse hyperbolic cosine; mathworld
|
|
299 |
# - requires argument >= 1
|
|
300 |
#----------------------------------------------------------------------
|
|
301 |
|
|
302 |
sub acosh($)
|
|
303 |
{
|
|
304 |
return log($_[0] + sqrt($_[0]**2-1));
|
|
305 |
}
|
|
306 |
|
|
307 |
#----------------------------------------------------------------------
|
20
|
308 |
# Gaussian random numbers
|
|
309 |
# - optional argument is seed
|
|
310 |
# - http://www.design.caltech.edu/erik/Misc/Gaussian.html
|
|
311 |
# - algorithm generates 2 random numbers
|
|
312 |
# - validated with plot '<count -o samp 1-100000 | list -Lfuns -c x=gaussRand() | Hist -cs 0.05 x',100000.0*0.05/sqrt(2*3.14159265358979)*exp(-x**2/2) wi li
|
|
313 |
#----------------------------------------------------------------------
|
|
314 |
|
|
315 |
{ my($y2);
|
|
316 |
my($srand_called);
|
|
317 |
|
|
318 |
sub gaussRand(@)
|
|
319 |
{
|
|
320 |
if (@_ && !$srand_called) {
|
|
321 |
srand(@_);
|
|
322 |
$srand_called = 1;
|
|
323 |
}
|
|
324 |
|
|
325 |
if (defined($y2)) {
|
|
326 |
my($temp) = $y2;
|
|
327 |
undef($y2);
|
|
328 |
return $temp;
|
|
329 |
}
|
|
330 |
|
|
331 |
my($x1,$x2,$w);
|
|
332 |
do {
|
|
333 |
$x1 = 2 * rand() - 1;
|
|
334 |
$x2 = 2 * rand() - 1;
|
|
335 |
$w = $x1**2 + $x2**2;
|
|
336 |
} while ($w >= 1);
|
|
337 |
|
|
338 |
$w = sqrt((-2 * log($w)) / $w);
|
|
339 |
$y2 = $x2 * $w;
|
|
340 |
return $x1 * $w;
|
|
341 |
}
|
|
342 |
|
|
343 |
}
|
|
344 |
|
|
345 |
#----------------------------------------------------------------------
|
0
|
346 |
|
|
347 |
1;
|