.
--- a/HISTORY
+++ b/HISTORY
@@ -1,11 +1,24 @@
#======================================================================
# H I S T O R Y
# doc: Thu May 7 13:12:05 2015
-# dlm: Thu May 7 13:12:52 2015
+# dlm: Fri Jun 12 07:24:38 2015
# (c) 2015 A.M. Thurnherr
-# uE-Info: 10 10 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 24 13 NIL 0 0 70 0 2 4 NIL ofnI
#======================================================================
May 7, 2015:
- - V6.0 [ants.pl] [.hg/hgrc] published for release of LADCPproc V1.2 (Slocum/Explorer processing
+ - V6.0 [ants.pl] [.hg/hgrc]
+ - published for release of LADCPproc V1.2 (Slocum/Explorer processing)
+
+May 17, 2015:
+ - V6.1 [ants.pl] [.hg/hgrc]
+ - added $skip to cFFT to for LADCP_w_spec and binpgrams
+ - NOT YET PUBLISHED
+May 18, 2015:
+ - added &antsFindParam() to [antsutils.pl] for LADCP_w_regrid
+ - added pulse length to [libLADCP.pl] for LADCP_VKE
+ - added rec indices to [lfit.pl] for binpgrams
+
+Jun 12, 2015:
+ - added &
--- a/ants.pl
+++ b/ants.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S . P L
# doc: Fri Jun 19 14:01:06 1998
-# dlm: Thu Oct 30 09:33:41 2014
+# dlm: Sun May 17 20:18:01 2015
# (c) 1998 A.M. Thurnherr
-# uE-Info: 22 55 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 17 34 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -13,12 +13,13 @@
# Jul 5, 2006: - removed `basename`
# Jul 19, 2006: - added error if exec($ANTS_PERL) fails
# Sep 24, 2012: - added support for $ANTSLIB
-# Oct 29, 2014: - added $antsLibVersion with compile-time version check
+# Oct 29, 2014: - added $antsLibVersion with compile-time version check (V6.0)
+# May 17, 2015: - updated to V6.1
exec($ENV{ANTS_PERL},$0,@ARGV),die("$ENV{ANTS_PERL}: $!")
if (defined($ENV{ANTS_PERL}) && $^X ne $ENV{ANTS_PERL});
-$antsLibVersion = 6.0;
+$antsLibVersion = 6.1;
die(sprintf("$0: obsolete library V%.1f; V%.1f required\n",
$antsLibVersion,$antsMinLibVersion))
if (!defined($antsMinLibVersion) || $antsMinLibVersion>$antsLibVersion);
--- a/antsexprs.pl
+++ b/antsexprs.pl
@@ -2,8 +2,8 @@
# A N T S E X P R S . P L
# (c) 2005 Andreas Thurnherr
# doc: Sat Dec 31 18:35:33 2005
-# dlm: Sat Mar 10 16:28:46 2012
-# uE-Info: 207 38 NIL 0 0 70 2 2 4 NIL ofnI
+# dlm: Fri May 15 20:12:34 2015
+# uE-Info: 183 56 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -42,6 +42,7 @@
# May 22, 2011: - made it work
# Feb 20, 2012: - BUG: quoting had not been implemented
# Mar 10, 2012: - added ${field..field} syntax to edit exprs
+# May 15, 2015: - BUG: -S did not work with :: %PARAMs
$DEBUG = 0;
@@ -80,8 +81,10 @@
$expr =~ /^\$?([\w\.]+)\s*>$/ ||
$expr =~ /^>\$?([\w\.]+)$/);
- if ($expr =~ /^(%?[\w\.]+):/ || $expr =~ /^(\$\d+):/) { # old -G syntax
+ $expr =~ s/::/QquOte/g; # new-style :: %PARAMs
+ if ($expr =~ /^(%?[\w\.]+):/ || $expr =~ /^(\$\d+):/) { # old -G syntax
my($fname) = $1; my($range) = $';
+ $fname =~ s/QquOte/::/g;
if ($range =~ /(.*)\.\.(.*)/) {
my($min) = ($1 eq '*') ? -1e99 : $1;
my($max) = ($2 eq '*') ? 1e99 : $2;
@@ -112,7 +115,8 @@
}
print(STDERR "-G AddrExpr = $expr\n") if ($DEBUG);
}
-
+ $expr =~ s/QquOte/::/g;
+
my($relop) = '<|<=|>|>=|!=|~=|<>|=='; # relational ops
my($comparee) = '-?%?\$?[\w\.\+\-]+'; # nums, fields, PARAMs
my($numvar) = '^[\w\.]+$'; # fields
@@ -176,7 +180,7 @@
$expr =~ s{\$\+$1}{(AnTsDtArEf\[$fnr\]-AnTsDtArEf0\[$fnr\])};
}
$expr =~ s{%%}{ AnTsPeRcEnT }g; # escape
- while ($expr =~ /%([\w\.]+)/) { # %PARAMs
+ while ($expr =~ /%([\w\.:]+)/) { # %PARAMs
my($p) = $1;
croak("$0: Undefined PARAM %$p\n")
unless defined($P{$p});
--- a/antsusage.pl
+++ b/antsusage.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S U S A G E . P L
# doc: Fri Jun 19 13:43:05 1998
-# dlm: Thu Mar 5 12:58:27 2015
+# dlm: Fri May 15 18:46:22 2015
# (c) 1998 A.M. Thurnherr
-# uE-Info: 161 0 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 161 75 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -158,6 +158,7 @@
# Jul 30, 2014: - added special args to -U)sage output
# Jan 30, 2015: - added &antsFunOpt()
# Jan 31, 2015: - made it work
+# May 15, 2015: - changed (()) semantics to expand only to existing files
# NOTES:
# - ksh expands {}-arguments with commas in them!!! Use + instead
@@ -371,11 +372,13 @@
sprintf("$pref%%0%dd$suff",length($1)) : "$pref%d$suff";
if ($2 > $1) {
for (my($i)=$1; $i<=$2; $i++) {
- push(@exp,sprintf($fmt,$i));
+ my($f) = sprintf($fmt,$i);
+ push(@exp,$f) if (-f $f);
}
} else {
for (my($i)=$1; $i>=$2; $i--) {
- push(@exp,sprintf($fmt,$i));
+ my($f) = sprintf($fmt,$i);
+ push(@exp,$f) if (-f $f);
}
}
} else {
--- a/antsutils.pl
+++ b/antsutils.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S U T I L S . P L
# doc: Fri Jun 19 23:25:50 1998
-# dlm: Tue Jul 22 20:35:50 2014
+# dlm: Fri Jun 12 07:31:08 2015
# (c) 1998 A.M. Thurnherr
-# uE-Info: 91 56 NIL 0 0 70 10 2 4 NIL ofnI
+# uE-Info: 103 66 NIL 0 0 70 10 2 4 NIL ofnI
#======================================================================
# Miscellaneous auxillary functions
@@ -99,6 +99,8 @@
# Sep 5, 2013: - FINALLY: added $pi
# May 23, 2014: - made ismember understand "123,1-10"
# Jul 22, 2014: - removed support for antsFnrNegativeOk
+# May 18, 2015: - added antsFindParam()
+# Jun 21, 2015: - added antsParam(), modified antsRequireParam()
# fnr notes:
# - matches field names starting with the string given, i.e. "sig" is
@@ -589,14 +591,42 @@
return @params;
} # sub antsfunusage()
+#----------------------------------------------------------------------
+
sub antsRequireParam($)
{
my($pn) = @_;
+ my($pv) = antsParam($pn);
croak("$0: required PARAM $pn not set\n")
- unless (defined($P{$pn}));
- return $P{$pn};
+ unless defined($pv);
+ return $pv;
+}
+
+
+sub antsFindParam($) # find parameter using RE (e.g. antsFindParam('dn\d\d'))
+{
+ my($re) = @_;
+ foreach my $k (keys(%P)) {
+ return ($k,$P{$k}) if ($k =~ /^$re$/);
+ }
+ return (undef,undef);
}
+sub antsParam($) # get parameter value for any ::-prefix
+{
+ my($pn) = @_;
+ my($nfound,$val);
+ foreach my $k (keys(%P)) {
+ next unless ($k eq $pn) || ($k =~ /::$pn$/);
+ $val = $P{$k};
+ $nfound++;
+ }
+ croak("$0: %PARAM $pn ambiguous\n")
+ if ($nfound > 1);
+ return $val;
+}
+
+#----------------------------------------------------------------------
{ my($term); # STATIC
--- a/fft.pl
+++ b/fft.pl
@@ -1,9 +1,9 @@
#======================================================================
# F F T . P L
# doc: Fri Mar 12 09:20:33 1999
-# dlm: Mon Jul 24 14:58:04 2006
+# dlm: Sun May 17 19:50:39 2015
# (c) 1999 A.M. Thurnherr
-# uE-Info: 241 36 NIL 0 0 72 66 2 4 NIL ofnI
+# uE-Info: 43 50 NIL 0 0 72 66 2 4 NIL ofnI
#======================================================================
# Notes:
@@ -40,6 +40,7 @@
# Feb 9, 2005: - added &phase_pos(), &phase_neg()
# Jul 1, 2006: - Version 3.3 [HISTORY]
# Jul 24, 2006: - modified to use $PRACTICALLY_ZERO
+# May 15, 2015: - added $skip argument to cFFT()
#----------------------------------------------------------------------
# FOUR1 routine
@@ -172,34 +173,35 @@
sub cFFT($$$) { return cFFT_bufR(\@ants_,@_); }
-sub cFFT_bufR($$$) # ($bufR, $rfnr, $ifnr, [$N]) => @coeff
+sub cFFT_bufR($$$) # ($bufR, $rfnr, $ifnr, [$N[, $skip]]) => @coeff
{
- my($bufR,$fnr,$ifnr,$N) = @_;
- my(@data,$i,$lastSet);
+ my($bufR,$fnr,$ifnr,$N,$skip) = @_;
+ my(@data,$i,$r,$lastSet);
unless ($N) { # $N not set
for ($N=1; $N <= $#ants_; $N <<= 1) {} # next greater pwroftwo
&antsInfo("(fft.pl) N set to $N")
unless ($N == $#ants_+1);
}
- for ($i=0; $i<$N && $i<=$#ants_; $i++) { # PAD
- last if (numberp($bufR->[$i][$fnr]) &&
- (isnan($ifnr) || numberp($bufR->[$i][$ifnr])));
+ for ($i=0,$r=$skip; $i<$N && $r<=$#ants_; $i++,$r++) { # PAD
+ last if (numberp($bufR->[$r][$fnr]) &&
+ (isnan($ifnr) || numberp($bufR->[$r][$ifnr])));
$data[2*$i] = 0;
$data[2*$i+1] = 0;
}
$lastSet = $i - 1;
&antsInfo("(fft.pl) WARNING: $i initial non-numbers padded with 0es!!!"),
$padded=1 if ($i);
- while ($i<$N && $i<=$#ants_) { # fill
- $i++,next unless (numberp($bufR->[$i][$fnr]) && # skip non-numbers
- (isnan($ifnr) || numberp($bufR->[$i][$ifnr])));
- croak("$0: (fft.pl) $lastSet, $i can't handle missing values ($bufR->[$lastSet+1][$fnr])!\n")
+ while ($i<$N && $r<=$#ants_) { # fill
+ $i++,$r++,next
+ unless (numberp($bufR->[$r][$fnr]) && # skip non-numbers
+ (isnan($ifnr) || numberp($bufR->[$r][$ifnr])));
+ croak("$0: (fft.pl) $lastSet, $i can't handle missing values ($bufR->[$lastSet+$skip+1][$fnr])!\n")
if ($lastSet != $i-1); # missing values
- $data[2*$i] = $bufR->[$i][$fnr]; # real
- $data[2*$i+1] = isnan($ifnr) ? 0 : $bufR->[$i][$ifnr]; # imag
+ $data[2*$i] = $bufR->[$r][$fnr]; # real
+ $data[2*$i+1] = isnan($ifnr) ? 0 : $bufR->[$r][$ifnr]; # imag
$lastSet = $i;
- $i++;
+ $i++; $r++;
}
&antsInfo("(fft.pl) WARNING: %d final non-numbers padded with 0es!!!",$i-$lastSet-1),
$padded=1 if ($i > $lastSet+1);
--- a/lfit.pl
+++ b/lfit.pl
@@ -1,9 +1,9 @@
#======================================================================
# L F I T . P L
# doc: Sat Jul 31 11:24:47 1999
-# dlm: Thu Jan 5 12:53:11 2012
+# dlm: Mon May 18 10:07:47 2015
# (c) 1999 A.M. Thurnherr
-# uE-Info: 19 60 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 20 61 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# LFIT routine from Numerical Recipes adapted to ANTS
@@ -17,6 +17,7 @@
# Jan 5, 2012: - BUG: non-numeric x/y were not handled correctly;
# this was only easily apparent when the last
# record contained non-numeric values
+# May 18, 2015: - added firstRec, lastRec parameters (V6.1)
# Notes:
# - x,y,sig are field numbers for data in $ants_
@@ -29,9 +30,11 @@
require "$ANTS/covsrt.pl";
require "$ANTS/gaussj.pl";
-sub lfit($$$$$$$)
+sub lfit(@)
{
- my($xfnr,$yfnr,$sig,$aR,$iaR,$covarR,$funcsR) = @_;
+ my($xfnr,$yfnr,$sig,$aR,$iaR,$covarR,$funcsR,$firstRec,$lastRec) = @_;
+ $firstRec = 0 unless defined($firstRec);
+ $lastRec = $#ants_ unless defined($lastRec);
my($i,$j,$k,$l,$m,$mfit); # int
my($ym,$wt,$sum,$sig2i,$chisq); # float
@@ -49,7 +52,7 @@
}
$beta[$j][1] = 0;
}
- for ($i=0; $i<=$#ants_; $i++) {
+ for ($i=$firstRec; $i<=$lastRec; $i++) {
next if ($antsFlagged[$i]);
next unless numberp($ants_[$i][$xfnr]) && numberp($ants_[$i][$yfnr]);
&$funcsR($i,$xfnr,\@afunc);
@@ -84,7 +87,7 @@
for ($j=0,$l=1;$l<=$#{$aR};$l++) {
$aR->[$l]=$beta[++$j][1] if ($iaR->[$l]);
}
- for ($i=0; $i<=$#ants_; $i++) {
+ for ($i=$firstRec; $i<=$lastRec; $i++) {
next if ($antsFlagged[$i]);
next unless numberp($ants_[$i][$xfnr]) && numberp($ants_[$i][$yfnr]);
&$funcsR($i,$xfnr,\@afunc);
--- a/libLADCP.pl
+++ b/libLADCP.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B L A D C P . P L
# doc: Wed Jun 1 20:38:19 2011
-# dlm: Wed Jun 26 11:47:45 2013
+# dlm: Mon May 18 09:49:19 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 211 0 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 49 12 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -22,6 +22,7 @@
# Oct 25, 2012: - renamed T_SM() to T_ASM()
# Jun 26. 2013: - added T_w_z()
# - added parameter checks to processing-specific corrections
+# May 18, 2015: - added pulse length to T_w() and T_w_z()
require "$ANTS/libvec.pl";
require "$ANTS/libfuns.pl";
@@ -40,16 +41,17 @@
#----------------------------------------------------------------------
#------------------------------------------------------------------------------
-# T_ravg(k,blen)
+# T_ravg(k,blen[,plen])
# - correct for range averaging due to finite pulse and finite receive window
-# - this version assumes bin-length == pulse-length
+# - when called with 2 arguments, bin-length == pulse-length assumed
#------------------------------------------------------------------------------
-sub T_ravg($$)
+sub T_ravg(@)
{
- my($k,$blen) =
- &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <pulse/bin-length[m]>',@_);
- return 1 / sinc($k*$blen/2/$PI)**4;
+ my($k,$blen,$plen) =
+ &antsFunUsage(-2,'ff','<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]]',@_);
+ $plen = $blen unless defined($plen);
+ return 1 / sinc($k*$blen/2/$PI)**2 / sinc($k*$plen/2/$PI)**2;
}
#-------------------------------------------------------------
@@ -175,7 +177,7 @@
}
#----------------------------------------------------------------------
-# T_w(k,blen,dz,range_max)
+# T_w(k,blen,plen,dz,range_max)
# - vertical-velocity method of Thurnherr (IEEE 2011)
# - range_max == 0 disables tilt correction
#----------------------------------------------------------------------
@@ -183,18 +185,18 @@
{ my(@fc);
sub T_w(@)
{
- my($k,$blen,$dz,$range_max) =
- &antsFunUsage(4,'ffff',
- '[vertical wavenumber[rad/s]] <ADCP bin size[m]> <output grid resolution[m]> <range max[m]>',
- \@fc,'k',undef,undef,undef,@_);
- croak("T_w($k,$blen,$dz,$range_max): bad parameters\n")
- unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);
- return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
+ my($k,$blen,$plen,$dz,$range_max) =
+ &antsFunUsage(5,'fffff',
+ '[vertical wavenumber[rad/s]] <ADCP bin length[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>',
+ \@fc,'k',undef,undef,undef,undef,@_);
+ croak("T_w($k,$blen,$plen,$dz,$range_max): bad parameters\n")
+ unless ($k>=0 && $blen>0 && $plen>0 && $dz>0 && $range_max>=0);
+ return T_ravg($k,$blen,$plen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
}
}
#----------------------------------------------------------------------
-# T_w_z(k,blen,dz,range_max)
+# T_w_z(k,blen,plen,dz,range_max)
# - vertical-velocity method of Thurnherr (IEEE 2011)
# - first differencing of gridded shear to calculate dw/dz
# - NB: grid-scale differentiation assumed
@@ -204,13 +206,13 @@
{ my(@fc);
sub T_w_z(@)
{
- my($k,$blen,$dz,$range_max) =
- &antsFunUsage(4,'ffff',
- '[vertical wavenumber[rad/s]] <ADCP bin size[m]> <output grid resolution[m]> <range max[m]>',
- \@fc,'k',undef,undef,undef,@_);
- croak("T_w_z($k,$blen,$dz,$range_max): bad parameters\n")
- unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);
- return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz);
+ my($k,$blen,$plen,$dz,$range_max) =
+ &antsFunUsage(5,'fffff',
+ '[vertical wavenumber[rad/s]] <ADCP bin size[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>',
+ \@fc,'k',undef,undef,undef,undef,@_);
+ croak("T_w_z($k,$blen,$plen,$dz,$range_max): bad parameters\n")
+ unless ($k>=0 && $blen>0 && $plen>0 && $dz>0 && $range_max>=0);
+ return T_ravg($k,$blen,$plen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz);
}
}
--- a/libfuns.pl
+++ b/libfuns.pl
@@ -1,9 +1,9 @@
#======================================================================
# 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
+# dlm: Thu Jun 4 17:56:37 2015
# (c) 1999 A.M. Thurnherr
-# uE-Info: 286 38 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 306 13 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -14,10 +14,14 @@
# Jan 25, 2001: - added f(), sgn()
# Apr 16, 2010: - added sinc()
# Sep 7, 2012: - added acosh()
+# Jun 4, 2015: - added gaussRand()
+# - made normal() more efficient
require "$ANTS/libvec.pl"; # rad()
#----------------------------------------------------------------------
+# gaussians/normal distribution
+#----------------------------------------------------------------------
sub gauss(@)
{
@@ -28,8 +32,8 @@
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);
+ my($sqrt2pi) = 2.506628274631;
+ return $area/($sqrt2pi*$sigma) * exp(-((($x-$mean) / $sigma)**2)/2);
}
#----------------------------------------------------------------------
@@ -287,5 +291,43 @@
}
#----------------------------------------------------------------------
+# Gaussian random numbers
+# - optional argument is seed
+# - http://www.design.caltech.edu/erik/Misc/Gaussian.html
+# - algorithm generates 2 random numbers
+# - 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
+#----------------------------------------------------------------------
+
+{ my($y2);
+ my($srand_called);
+
+sub gaussRand(@)
+{
+ if (@_ && !$srand_called) {
+ srand(@_);
+ $srand_called = 1;
+ }
+
+ if (defined($y2)) {
+ my($temp) = $y2;
+ undef($y2);
+ return $temp;
+ }
+
+ my($x1,$x2,$w);
+ do {
+ $x1 = 2 * rand() - 1;
+ $x2 = 2 * rand() - 1;
+ $w = $x1**2 + $x2**2;
+ } while ($w >= 1);
+
+ $w = sqrt((-2 * log($w)) / $w);
+ $y2 = $x2 * $w;
+ return $x1 * $w;
+}
+
+}
+
+#----------------------------------------------------------------------
1;