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