1 #====================================================================== |
1 #====================================================================== |
2 # L I B F U N S . P L |
2 # L I B F U N S . P L |
3 # doc: Wed Mar 24 11:49:13 1999 |
3 # doc: Wed Mar 24 11:49:13 1999 |
4 # dlm: Fri Sep 7 11:11:09 2012 |
4 # dlm: Thu Jun 4 17:56:37 2015 |
5 # (c) 1999 A.M. Thurnherr |
5 # (c) 1999 A.M. Thurnherr |
6 # uE-Info: 286 38 NIL 0 0 72 2 2 4 NIL ofnI |
6 # uE-Info: 306 13 NIL 0 0 72 2 2 4 NIL ofnI |
7 #====================================================================== |
7 #====================================================================== |
8 |
8 |
9 # HISTORY: |
9 # HISTORY: |
10 # Mar 24, 1999: - copied from the c-version of NR |
10 # Mar 24, 1999: - copied from the c-version of NR |
11 # Mar 26, 1999: - added stuff for better [./fit] |
11 # Mar 26, 1999: - added stuff for better [./fit] |
12 # Sep 18, 1999: - argument typechecking |
12 # Sep 18, 1999: - argument typechecking |
13 # Oct 04, 1999: - added gauss(), normal() |
13 # Oct 04, 1999: - added gauss(), normal() |
14 # Jan 25, 2001: - added f(), sgn() |
14 # Jan 25, 2001: - added f(), sgn() |
15 # Apr 16, 2010: - added sinc() |
15 # Apr 16, 2010: - added sinc() |
16 # Sep 7, 2012: - added acosh() |
16 # Sep 7, 2012: - added acosh() |
|
17 # Jun 4, 2015: - added gaussRand() |
|
18 # - made normal() more efficient |
17 |
19 |
18 require "$ANTS/libvec.pl"; # rad() |
20 require "$ANTS/libvec.pl"; # rad() |
19 |
21 |
|
22 #---------------------------------------------------------------------- |
|
23 # gaussians/normal distribution |
20 #---------------------------------------------------------------------- |
24 #---------------------------------------------------------------------- |
21 |
25 |
22 sub gauss(@) |
26 sub gauss(@) |
23 { |
27 { |
24 my($x,$peak,$mean,$efs) = &antsFunUsage(4,"ffff","x, peak, mean, e-folding scale",@_); |
28 my($x,$peak,$mean,$efs) = &antsFunUsage(4,"ffff","x, peak, mean, e-folding scale",@_); |
26 } |
30 } |
27 |
31 |
28 sub normal(@) |
32 sub normal(@) |
29 { |
33 { |
30 my($x,$area,$mean,$sigma) = &antsFunUsage(4,"ffff","x, area, mean, stddev",@_); |
34 my($x,$area,$mean,$sigma) = &antsFunUsage(4,"ffff","x, area, mean, stddev",@_); |
31 my($pi) = 3.14159265358979; |
35 my($sqrt2pi) = 2.506628274631; |
32 return $area/(sqrt(2*$pi)*$sigma) * exp(-((($x-$mean) / $sigma)**2)/2); |
36 return $area/($sqrt2pi*$sigma) * exp(-((($x-$mean) / $sigma)**2)/2); |
33 } |
37 } |
34 |
38 |
35 #---------------------------------------------------------------------- |
39 #---------------------------------------------------------------------- |
36 # &f(lat) calculate coriolis param |
40 # &f(lat) calculate coriolis param |
37 #---------------------------------------------------------------------- |
41 #---------------------------------------------------------------------- |
285 { |
289 { |
286 return log($_[0] + sqrt($_[0]**2-1)); |
290 return log($_[0] + sqrt($_[0]**2-1)); |
287 } |
291 } |
288 |
292 |
289 #---------------------------------------------------------------------- |
293 #---------------------------------------------------------------------- |
|
294 # Gaussian random numbers |
|
295 # - optional argument is seed |
|
296 # - http://www.design.caltech.edu/erik/Misc/Gaussian.html |
|
297 # - algorithm generates 2 random numbers |
|
298 # - 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 |
|
299 #---------------------------------------------------------------------- |
|
300 |
|
301 { my($y2); |
|
302 my($srand_called); |
|
303 |
|
304 sub gaussRand(@) |
|
305 { |
|
306 if (@_ && !$srand_called) { |
|
307 srand(@_); |
|
308 $srand_called = 1; |
|
309 } |
|
310 |
|
311 if (defined($y2)) { |
|
312 my($temp) = $y2; |
|
313 undef($y2); |
|
314 return $temp; |
|
315 } |
|
316 |
|
317 my($x1,$x2,$w); |
|
318 do { |
|
319 $x1 = 2 * rand() - 1; |
|
320 $x2 = 2 * rand() - 1; |
|
321 $w = $x1**2 + $x2**2; |
|
322 } while ($w >= 1); |
|
323 |
|
324 $w = sqrt((-2 * log($w)) / $w); |
|
325 $y2 = $x2 * $w; |
|
326 return $x1 * $w; |
|
327 } |
|
328 |
|
329 } |
|
330 |
|
331 #---------------------------------------------------------------------- |
290 |
332 |
291 1; |
333 1; |