.
new file mode 100755
--- /dev/null
+++ b/INDEX
@@ -0,0 +1,213 @@
+======================================================================
+ I N D E X
+ doc: Wed Jun 18 09:46:58 1997
+ dlm: Tue Oct 16 09:12:43 2012
+ (c) 1997 Andreas Thurnherr
+ uE-Info: 122 0 NIL 0 0 72 2 2 4 NIL ofnI
+======================================================================
+
+NOTES:
+ - for documentation see [doc/INDEX]
+ - utilities not in ubtest are marked with !!!
+
+=Utilities=
+
+-Oceanographic-
+
+[ARGOSfloatdrift] split ARGOS float track into subsurface/surface portions
+[CMrecavg] record-average velocity & stderr
+[CTD_w] calculate w from isopycnal displacement
+[dist] distance on Earth's surface
+[distn] station distance / mid point btw 2 profiles
+[distcircle] calc geographic coords a given distance from point
+[distfrom] calculate distance to particular station
+[dreckon] dead reckoning (calc geo coords from length displacement)
+[dynmodes]!!! calculate dynamical modes
+[expandtrack] build interpolate track between stations
+[gamma_n] calculate neutral densities
+[gvel] calculate geostrophic velocities
+[interpstn] interpolate between 2 stations
+[isopycnal_TS] calc TS props on isopycnal
+[LADCPfs]!!! calculate LADCP finestructure from [binpgrams] output
+[N]!!! calculate buoyancy frequency
+[NMEA2latlon] extract lat/lon info from NMEA strings
+[NODCsplit]!!! split NODC_SD2 file into components
+[repeatstations]!!! find repeat stations
+[Ri]!!! calculate Richardson number
+[Rrho]!!! stability ratio profile
+[thorpe77]!!! calc Thorpe fluctuations
+[tidalphasediff] calc M2 & K1 phase differences between two profiles
+[trackdist] calculate distance for track-files
+[tracksplit] split data set into linear tracks
+[veldiff] compare velocity profiles
+ [.interp.linear] linear interpolation
+ [.interp.poly] polygonal interpolation
+ [.interp.spline] spline interpolation
+ [.interp.nnbr] nearest-neighbor
+ [.interp.ADCP] RDI ADCP velocity sampling (triangular instrument response)
+[waterdepth]!!! get water depth from Smith and Sandwell topography
+[xover] get crossover stations
+[xpgrams]!!! extract pgrams from [binpgrams] [LADCPfs] output for plotting
+[yoyo] splits yoyo file into individual casts
+
+-Model Fitting-
+
+[lfit] linear least squares
+[match]!!! match two data sets --- UNFINISHED
+ [.interp.linear] linear interpolation
+ [.interp.poly] polygonal interpolation
+ [.interp.spline] spline interpolation
+ [.interp.nnbr] nearest-neighbor
+[mfit]!!! robust linear regression
+[lmfit] non-linear least-squares (Levenberg-Marquardt)
+ [.lmfit.poly] fit polynomial
+ [.lmfit.gauss] fit Gaussian
+ [.lmfit.exp] fit exponential
+[lsfit] generalized linear least-squares
+ [.lsfit.poly] fit polynomial
+ [.lsfit.bilin] bi-linear fit
+
+-Stats-
+
+[avg] avg/stddev/stderr/ndata/outliers/mean-sq amp/absolute mag/rms/gaussian avg
+[avgr] avg correlation coefficient
+[bootstrap] calc bootstrap confidence interval
+[Hist] histogram
+[histeq] histogram equalization
+[max] maximum
+[median] median
+[min] minimum
+[Minmax] min/max combo
+
+-Calculations-
+
+[abc] PERL floating point replacement for expr(1)/bc(1)
+[count] as in 1, 2, 3, ...
+[fdiff] first differencing
+[fillgaps] general interpolation of missing data
+ [.nminterp.linear]
+ [.nminterp.spline]
+[integrate] integrate data
+[sum] calc sum of field
+!!![wdiff] window differencing
+
+-Filters-
+
+[fftfilt] simplistic frequency domain digital filter
+[fgauss]!!! running Gaussian filter
+[fmean]!!! running mean filter
+[fmedian]!!! median filter (OBSOLETE)
+[rfilt]!!! general running filter (mean, median, min, max, ...)
+
+-Spectral Methods-
+
+[binpgrams]!!! calc peridograms in (possibly overlapping) bins
+[lagcorr] (lag) correlation / autocorrelation
+[pgram] calc periodograms (power spectra)
+
+-Data Selection/Handling-
+
+[ANTS2mat] Matlab export
+[bindata] bin data & calc stats
+[binmerge]!!! merge data into binned file (e.g. from [binpgrams])
+[Cat] ants version of cat(1)
+[crossings]!!! get records when specific field val is crossed
+[data]! get info about data
+[ded] data editor
+[distfromtrack]!!! calculate minimum perpendicular distance to piece-wise linear track
+[Echo] ants versions of echo(1)
+[fileavg]!!! avg, stddev files record-by-record
+[filediff]! calc difference between files record-by-record
+[filestats]!!! calc arbitrary stats of files record-by-record
+[fnr] extracts field# from header or from Layout file
+[gextrema] select extrema
+[glist] select records from list
+[glmax]!!! find local maxima
+[gmonot]!!! select monotonic records (remove spurious inversions)
+[gpoly] select points inside a polygon
+[how]!!! extract command sequence from ANTS header
+[importCNV] read ASCII and binary SeaBird CNV files
+[importfixed] import fixed-size records
+[list] list data and metadata
+[listNC] list netcdf(1) data and metadata
+[merge] merge files by matching numeric values
+[NCode] encode ANTS to netcdf
+[resample] resample data
+ [.interp.linear] linear interpolation
+ [.interp.poly] polygonal interpolation
+ [.interp.spline] spline interpolation
+ [.interp.nnbr] nearest-neighbor
+ [.interp.fillnan] fill with nans
+[reverse] reverse order of records in file
+[scantext] extract numbers from free-format text files
+[Sort] sort data
+[Split] split data file according to field value
+[splitCNV] split ASCII SeaBird CNV file
+[splitNC] extract multiple ANTS files from single netCDF
+[stackplots] plot multiple files with offsets
+[subsample] nearest-neighbors; does not need monotonic x field
+[varbindata] bin data in variable-sized bins
+
+-Misc-
+
+[fixfmt]!!! fix path in ANTS header & nan format
+[tile]!!! tile eps-files by 8
+[tabulate]!!! make visually pleasing table from data
+
+-Debugging-
+
+[argtest] test argument typechecking stuff
+
+
+=ANTS Libraries=
+
+[ants.pl] common stuff
+ [antsio.pl] input/output handling
+ [antsusage.pl] usage handling
+ [antsutils.pl] miscellaneous
+[antsfilters.pl] library for data filtering
+[antsintegrate.pl] integrator
+
+=Libraries (for -L)=
+
+[libCPT.pl] GMT cpt file parsing
+[libEOS83.pl] equation of state
+[libfuns.pl] chapter 6 from NR: special functions; own stuff
+[libgamma.pl] Jacket + McDougall gamma_n
+[libGHP.pl] Gregg-Henyey-Polzin parameterization
+[libGM.pl] Garrett & Munk spectrum
+[libLADCP.pl] LADCP-related funs
+[libNODC.pl]!!! NODC SD2 conversions
+[libPOSIX.pl] POSIX functions
+[libstats.pl] significance tests
+[libubtest.pl] testing stuff
+[libvec.pl] vector stuff, including distance on globe
+[libWOCE.pl] WOCE conversions
+[libconv.pl] conversions
+
+=Numerical Recipes Routines=
+
+[covsrt.pl] for [lmfit]
+[fft.pl] FFT
+[gaussj.pl] Gaussian eliminiation
+[lfit.pl] linear least squares
+[mrqcof.pl] for [lmfit]
+[mrqmin.pl] for [lmfit]
+[nrutil.pl] aux utilities (vector/matrix/macros)
+[pearsn.pl] for [corr]
+[polint.pl] for [.interp.poly]
+[pythag.pl] pythagoras
+
+=GMT Utilities=
+
+[adjustROI]!!! make ROI compatible with grd file
+[CPTcolor]!!! extract colour from cpt-file
+[CPTcontours]!!! extract contour-levels from cpt-file
+[grdjoin]!!! join multiple compatible grd files
+[lmcont]!!! linear contouring of (quasi-)monotonic data
+[mkCPT] makecpt(l) replacement
+[rectangulate]!!! construct GMT multiseg rectangles from x/y/z/width info
+[psbath]!!! build bathymetric map; similar to pscoast
+[psROI]!!! psxy frontent to plot regions of interest
+[psSamp]!!! generate GMT-compatible PS output for pie-wedge (repeat-)station data
+
--- 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: Wed Jul 5 15:37:12 2006
+# dlm: Mon Sep 24 12:41:50 2012
# (c) 1998 A.M. Thurnherr
-# uE-Info: 18 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 25 76 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -12,15 +12,24 @@
# Jul 3, 2006: - added support for ANTS_PERL
# Jul 5, 2006: - removed `basename`
# Jul 19, 2006: - added error if exec($ANTS_PERL) fails
+# Sep 24, 2012: - added support for $ANTSLIB
exec($ENV{ANTS_PERL},$0,@ARGV),die("$ENV{ANTS_PERL}: $!")
if (defined($ENV{ANTS_PERL}) && $^X ne $ENV{ANTS_PERL});
-($ANTS) = ($0 =~ m{^(.*)/[^/]*$}) unless defined($ANTS);
-
-require "$ANTS/antsusage.pl";
-require "$ANTS/antsio.pl";
-require "$ANTS/antsutils.pl";
-require "$ANTS/antsexprs.pl";
+if (defined($ANTSLIB)) { # new style (V5)
+ require "$ANTSLIB/antsusage.pl";
+ require "$ANTSLIB/antsio.pl";
+ require "$ANTSLIB/antsutils.pl";
+ require "$ANTSLIB/antsexprs.pl";
+ $ANTS = $ANTSLIB; # backward compatibility
+} elsif (defined($ANTS)) { # old style
+ require "$ANTS/antsusage.pl";
+ require "$ANTS/antsio.pl";
+ require "$ANTS/antsutils.pl";
+ require "$ANTS/antsexprs.pl";
+} else {
+ die("neither \$ANTS nor \$ANTSLIB defined\n");
+}
1;
--- a/antsio.pl
+++ b/antsio.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S I O . P L
# doc: Fri Jun 19 19:22:51 1998
-# dlm: Thu Apr 26 09:01:50 2012
+# dlm: Wed Oct 10 12:36:29 2012
# (c) 1998 A.M. Thurnherr
-# uE-Info: 200 107 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 598 0 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -577,23 +577,25 @@
$OEparamsOnly = 1;
for (my($if)=my($of)=0; $if<@antsOutExprs; $if++,$of++) {
- if ($antsOutExprs[$if] =~ m{^%([\w\.]+)$}) { # %PARAM
+# if ($antsOutExprs[$if] =~ m{^%([\w\.]+)$}) { # %PARAM
+ if ($antsOutExprs[$if] =~ m{^%([^=]+)$}) { # %PARAM
$ofn[$of] = $1;
$OEparam[$of] = 1;
- } elsif ($antsOutExprs[$if] =~ m{^[\w\.]+$}) { # single field
- undef($OEparamsOnly);
- $ofn[$of] = $antsOutExprs[$if];
- $OEfield[$of] = &outFnr($antsOutExprs[$if]);
} elsif ($antsOutExprs[$if] eq \'$@\') { # all fields
undef($OEparamsOnly);
for (my($i)=0; $i<@ofn_buf; $i++,$of++) {
$ofn[$of] = $ofn_buf[$i];
$OEfield[$of] = $i;
}
+# } elsif ($antsOutExprs[$if] =~ m{^[\w\.]+$}) { # single field
+ } elsif ($antsOutExprs[$if] =~ m{^[^=]+$}) { # single field
+ undef($OEparamsOnly);
+ $ofn[$of] = $antsOutExprs[$if];
+ $OEfield[$of] = &outFnr($antsOutExprs[$if]);
} else { # expression
undef($OEparamsOnly);
my($expr);
- ($ofn[$of],$expr) = ($antsOutExprs[$if] =~ m{^([\w\.]*)=(.*)$});
+ ($ofn[$of],$expr) = ($antsOutExprs[$if] =~ m{^([^=]*)=(.*)$});
croak("$0: cannot parse -F $antsOutExprs[$if]\n")
unless defined($expr);
my(@tmp) = @antsLayout;
--- 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: Mon Feb 13 19:57:03 2012
+# dlm: Mon Oct 29 12:08:54 2012
# (c) 1998 A.M. Thurnherr
-# uE-Info: 441 0 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 153 62 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -150,6 +150,7 @@
# resulted in always extending the layout, even when the field already
# existed)
# Feb 13, 2012: - antsNewFieldOpt simplified by using 2nd arg to fnrNoErr
+# Oct 29, 2012: - diabled "no file" messages on special args
# NOTES:
# - ksh expands {}-arguments with commas in them!!! Use + instead
@@ -326,19 +327,19 @@
for (my($i)=$1; $i<=$2; $i++) {
my($f) = sprintf($fmt,$i);
if (-f $f) { push(@exp,$f); }
- else { &antsInfo("$ARGV[$ai]: no file <$f>"); }
+# else { &antsInfo("$ARGV[$ai]: no file <$f>"); }
}
} else {
for (my($i)=$1; $i>=$2; $i--) {
my($f) = sprintf($fmt,$i);
if (-f $f) { push(@exp,$f); }
- else { &antsInfo("$ARGV[$ai]: no file <$f>"); }
+# else { &antsInfo("$ARGV[$ai]: no file <$f>"); }
}
}
} else {
my($f) = "$pref$range$suff";
if (-f $f) { push(@exp,$f); }
- else { &antsInfo("$ARGV[$ai]: no file <$f>"); }
+# else { &antsInfo("$ARGV[$ai]: no file <$f>"); }
}
@exp = ($ARGV[$ai]) # make sure it *was* special arg
unless (@exp);
--- 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: Thu May 31 09:13:03 2012
+# dlm: Wed Oct 24 09:56:52 2012
# (c) 1998 A.M. Thurnherr
-# uE-Info: 157 9 NIL 0 0 70 10 2 4 NIL ofnI
+# uE-Info: 259 43 NIL 0 0 70 10 2 4 NIL ofnI
#======================================================================
# Miscellaneous auxillary functions
@@ -94,6 +94,7 @@
# - BUG: fnrNoErr disregarded exact flag for external layouts
# May 16, 2012: - adapted to V5.0
# May 31, 2012: - changed ismember() semantics for use in psSamp
+# Jun 12, 2012: - added &compactList()
# fnr notes:
# - matches field names starting with the string given, i.e. "sig" is
@@ -255,7 +256,7 @@
return $num;
}
-sub log10 { my $n = shift; return log($n)/log(10); } # c.v. perlfunc(1)
+sub log10 { my $n = shift; return ($n>0) ? log($n)/log(10) : nan; } # c.v. perlfunc(1)
#----------------------------------------------------------------------
@@ -463,6 +464,43 @@
}
#----------------------------------------------------------------------
+# deal with lists of numbers
+#----------------------------------------------------------------------
+
+sub compactList(@)
+{
+ my(@out);
+ my($seqStart);
+ my($lv) = -9e99;
+
+ foreach my $v (@_) {
+ if (numberp($v)) {
+ if ($v == $lv+1) { # we're in a sequence
+ $seqStart = $lv # record beginning value
+ unless defined($seqStart);
+ } elsif (defined($seqStart)) { # we've just completed a sequence
+ pop(@out);
+ push(@out,"$seqStart-$lv");
+ push(@out,$v);
+ undef($seqStart);
+ } else { # not in a sequence
+ push(@out,$v);
+ }
+ $lv = $v;
+ } else {
+ push(@out,$v);
+ $lv = -9e99;
+ }
+ }
+ if (defined($seqStart)) { # list ends with a sequence
+ pop(@out);
+ push(@out,"$seqStart-$lv");
+ }
+
+ return @out;
+}
+
+#----------------------------------------------------------------------
# Misc funs
#----------------------------------------------------------------------
new file mode 100644
--- /dev/null
+++ b/libGHP.pl
@@ -0,0 +1,46 @@
+#======================================================================
+# L I B G H P . P L
+# doc: Fri Sep 7 09:56:08 2012
+# dlm: Mon Oct 22 13:10:47 2012
+# (c) 2012 A.M. Thurnherr
+# uE-Info: 11 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Sep 7, 2012: - created
+# Oct 22, 2012: - cosmetics
+
+require "$ANTS/libfuns.pl"; # arccosh
+require "$ANTS/libGM.pl"; # GM_N0
+
+#----------------------------------------------------------------------------
+# h(R_omega) correction factor for different shear/strain (R_omega) ratios
+# - this version from Kunze et al. (2006)
+# - R_omega:
+# = (N^2-omega^2)/N^2 (omega^2+f^2)/(omega^2-f^2) Polzin et al. 1995
+# = (omega^2+f^2)/(omega^2-f^2) Kunze et al. 2006
+# - the Kunze et al formulation is presumably an approximation valid for
+# omega<<N
+#----------------------------------------------------------------------------
+
+sub h1($)
+{
+ my($R_omega) = @_;
+ return 3*($R_omega+1) / (2*sqrt(2)*$R_omega*sqrt($R_omega-1));
+}
+
+#----------------------------------------------------------------------------
+# j(f,N) correction factor for latitude
+# - this version from Kunze et al. (2006)
+#----------------------------------------------------------------------------
+
+sub j(@)
+{
+ my($f,$N) = @_;
+ return 0 if ($f == 0);
+ $N = $GM_N0 unless defined($N);
+ my($f30) = &f(30);
+ return $f*acosh($N/$f) / ($f30*acosh($GM_N0/$f30));
+}
+
+1;
--- a/libGM.pl
+++ b/libGM.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B G M . P L
# doc: Sun Feb 20 14:43:47 2011
-# dlm: Sun Apr 1 11:29:53 2012
+# dlm: Fri Sep 14 12:35:40 2012
# (c) 2011 A.M. Thurnherr
-# uE-Info: 53 1 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 61 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -15,12 +15,23 @@
# affects only long wavelengths
# - return nan for omega outside internal-wave band
# Mar 29, 2012: - re-wrote using definition of B(omega) from Munk (1981)
+# Aug 23, 2012: - cosmetics?
+# Sep 7, 2012: - made N0, E0, b, jstar global
require "$ANTS/libEOS83.pl";
my($pi) = 3.14159265358979;
#======================================================================
+# Global Constants
+#======================================================================
+
+$GM_N0 = 5.24e-3; # rad/s # reference stratification (from Gregg + Kunze, 1991)
+$GM_E0 = 6.3e-5; # dimensionless # spectral level (Munk 1981)
+$GM_b = 1300; # m # pycnocline e-folding scale
+$GM_jstar = 3; # dimless # peak mode number
+
+#======================================================================
# Vertical velocity spectral density
#
# Units: K.E. per frequency per wavenumber [m^2/s^2*1/s*1/m = m^3/s]
@@ -47,13 +58,9 @@
{
my($j,$N,$omega) = @_;
- my($b) = 1300; #m # stratification e-folding scale (Munk 81)
- my($N0) = 5.2e-3; #rad/s # extrapolated to surface value (Munk 81)
-
-# print(STDERR "omega = $omega, N = $N\n");
return defined($omega)
- ? $pi / $b * sqrt(($N**2 - $omega**2) / ($N0**2 - $omega**2)) * $j
- : $pi * $j * $N / ($b * $N0); # valid, except in vicinity of buoyancy turning frequency (p. 285)
+ ? $pi / $GM_b * sqrt(($N**2 - $omega**2) / ($GM_N0**2 - $omega**2)) * $j
+ : $pi * $j * $N / ($GM_b * $GM_N0); # valid, except in vicinity of buoyancy turning frequency (p. 285)
}
sub B($) # structure function (omega dependence)
@@ -67,18 +74,35 @@
sub Sw($$$$)
{
- my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[1/m]> <lat[deg]> <N[rad/s]>',@_);
+ my($omega,$m,$lat,$N) = &antsFunUsage(4,'fff','<frequency[1/s]> <vertical wavenumber[rad/m]> <lat[deg]> <N[rad/s]>',@_);
local($f) = abs(&f($lat));
return nan if ($omega < $f || $omega > $N);
- my($E0) = 6.3e-5; # dimensionless spectral level
- my($j_star) = 3; # peak mode number
- my($b) = 1300; #m # pycnocline lengthscale
+ my($GM_b) = 1300; #m # pycnocline lengthscale
+
+ my($mstar) = &m($GM_jstar,$N,$omega);
+
+ return $GM_E0 * $GM_b * 2 * $f**2/$omega**2/B($omega) * $GM_jstar / ($m+$mstar)**2;
+}
+
+#----------------------------------------------------------------------
+# GM76, as per Gregg and Kunze (JGR 1991)
+# - beta is vertical wavenumber (m above)
+#----------------------------------------------------------------------
- my($mstar) = &m($j_star,$N,$omega);
+sub Su($$)
+{
+ my($beta,$N) = @_;
- return $E0 * $b * 2 * $f**2/$omega**2/B($omega) * $j_star / ($m+$mstar)**2;
+ my($beta_star) = &m($GM_jstar,$N); # A3
+ return 3*$GM_E0*$GM_b**3*$GM_N0**2 / (2*$GM_jstar*$pi) / (1+$beta/$beta_star)**2; # A2
+}
+
+sub Su_z($$)
+{
+ my($beta,$N) = &antsFunUsage(2,'ff','<vertical wavenumber[rad/m]> <N[rad/s]>',@_);
+ return $beta**2 * &Su($beta,$N);
}
1;
--- 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 Jan 18 18:46:33 2012
+# dlm: Thu Oct 25 21:24:59 2012
# (c) 2011 A.M. Thurnherr
-# uE-Info: 103 19 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 22 41 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -16,101 +16,170 @@
# Jan 5, 2012: - removed S(), which is just pwrdens/N^2 (rather than
# pwrdens/N^2/(2pi) as I erroneously thought)
# Jan 18, 2012: - added T_VI_alt() to allow assessment of tilt correction extrema
+# Aug 22, 2012: - added documentation
+# - added T_w()
+# Sep 24, 2012: - made "k" argument default in T_w()
+# Oct 25, 2012: - renamed T_SM to T_ASM
require "$ANTS/libvec.pl";
require "$ANTS/libfuns.pl";
+#------------------------------------------------------------------------------
+# Spectral corrections for LADCP data
+#
+# NOTES:
+# - see Polzin et al. (JAOT 2002), Thurnherr (JAOT 2012)
+# - to correct, multiply power densities (or power, I think) with corrections
+# - apply to down-/up-cast data only
+#------------------------------------------------------------------------------
+
#----------------------------------------------------------------------
-# Polzin et al., JAOT 2002 LADCP shear corrections
+# 1. Corrections for individual data acquisition and processing steps
#----------------------------------------------------------------------
-# NOTES:
-# - apply to downcast data only
-
-#----------------------------------------------------------------------
-# individual corrections
-#----------------------------------------------------------------------
-
-# NB: Dzb = (Dzt == Dzr) assumed
+#------------------------------------------------------------------------------
+# T_ravg(k,blen)
+# - correct for range averaging due to finite pulse and finite receive window
+# - this version assumes bin-length == pulse-length
+#------------------------------------------------------------------------------
sub T_ravg($$)
{
- my($kz,$Dzb) =
+ my($k,$blen) =
&antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <pulse/bin-length[m]>',@_);
- return 1 / sinc($kz*$Dzb/2/$PI)**4;
+ return 1 / sinc($k*$blen/2/$PI)**4;
}
+#-------------------------------------------------------------
+# T_fdiff(k,dz)
+# - correct for first differencing on a grid with dz spacing
+#-------------------------------------------------------------
sub T_fdiff($$)
{
- my($kz,$Dzd) =
+ my($k,$dz) =
&antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_);
- return 1 / sinc($kz*$Dzd/2/$PI)**2;
+ return 1 / sinc($k*$dz/2/$PI)**2;
}
+#------------------------------------------------------------
+# T_interp(k,blen,dz)
+# - correct for CODAS gridding-with-interpolation algorithm
+# - ONLY USED IN UH SOFTWARE
+#------------------------------------------------------------
sub T_interp($$$)
{
- my($kz,$Dzb,$Dzg) =
+ my($k,$blen,$dz) =
&antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_);
- return 1 / sinc($kz*$Dzb/2/$PI)**4 / sinc($kz*$Dzg/2/$PI)**2;
+ return 1 / sinc($k*$blen/2/$PI)**4 / sinc($k*$dz/2/$PI)**2;
}
+#-------------------------------------------------------------------------
+# T_binavg(k,dz)
+# - correct for simple bin averaging
+# - Polzin et al. suggest to use blen instead of dz; this must be a typo
+#-------------------------------------------------------------------------
-# NB: Polzin et al claim that Dz should be ADCP bin size, which does not seem to make sense
sub T_binavg($$)
{
- my($kz,$Dzg) =
+ my($k,$dz) =
&antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_);
- return 1 / sinc($kz*$Dzg/2/$PI)**2;
+ return 1 / sinc($k*$dz/2/$PI)**2;
}
+#--------------------------------------------------------------------------------
+# T_tilt(k,d')
+# - d' is a length scale that depends on tilt stats and range max
+# - on d' == 0, T_tilt() == 1, i.e. the correction is disabled
+# - d' = dprime(range_max)
+# - is from a quadratic fit to three data points given by Polzin et al.
+# - see Thurnherr (J. Tech. 2012) for notes
+# - on range_max == 0, d' == 0, i.e. the correction is disabled
+#--------------------------------------------------------------------------------
sub T_tilt($$)
{
- my($kz,$dprime) =
+ my($k,$dprime) =
&antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_);
- return 1 / sinc($kz*$dprime/2/$PI)**2;
+ return $dprime ? 1 / sinc($k*$dprime/2/$PI)**2 : 1;
+}
+
+sub dprime($)
+{
+ return $_[0] ? -1.2 + 0.0857 * $_[0] - 0.000136 * $_[0]**2 : 0;
+}
+
+#======================================================================
+# 2. Processing-Specific Corrections
+#======================================================================
+
+#----------------------------------------------------------------------
+# T_UH(k,blen,dz,range_max)
+# - UH implementation of the shear method (WOCE standard)
+# - range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+sub T_UH($$$$)
+{
+ my($k,$blen,$dz,$range_max) =
+ &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+ return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max));
}
#----------------------------------------------------------------------
-# convenient combinations
+# T_ASM(k,blen,dz,range_max)
+# - re-implemented shear method with simple depth binning
+# - range_max == 0 disables tilt correction
#----------------------------------------------------------------------
-sub LADCP_tilt_dprime($)
+sub T_ASM($$$$)
{
- return -1.2 + 0.0857 * $_[0] - 0.000136 * $_[0]**2;
-}
-
-sub T_UH($$$$)
-{
- my($kz,$blen,$grez,$maxrange) =
- &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <max range[m]>',@_);
- return T_ravg($kz,$blen) * T_fdiff($kz,$blen) * T_interp($kz,$blen,$grez) * T_tilt($kz,LADCP_tilt_dprime($maxrange));
+ my($k,$blen,$dz,$range_max) =
+ &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+ return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
}
-sub T_SM($$$$)
-{
- my($kz,$blen,$grez,$maxrange) =
- &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <max range[m]>',@_);
- return T_ravg($kz,$blen) * T_fdiff($kz,$blen) * T_binavg($kz,$grez) * T_tilt($kz,LADCP_tilt_dprime($maxrange));
-}
+#------------------------------------------------------------
+# T_VI(k,blen,preavg_dz,grid_dz,range_max)
+# T_VI_alt(k,blen,preavg_dz,grid_dz,dprime)
+# - velocity inversion method of Visbeck (J. Tech., 2002)
+# - only valid if pre-averaging into superensembles is used
+# - range_max == 0 disables tilt correction
+#------------------------------------------------------------
sub T_VI($$$$$)
{
- my($kz,$blen,$sel,$grez,$maxrange) =
- &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <max range[m]>',@_);
- return T_VI_alt($kz,$blen,$sel,$grez,LADCP_tilt_dprime($maxrange));
+ my($k,$blen,$sel,$dz,$range_max) =
+ &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <range max[m]>',@_);
+ return T_VI_alt($k,$blen,$sel,$dz,dprime($range_max));
}
sub T_VI_alt($$$$$)
{
- my($kz,$blen,$sel,$grez,$dprime) =
+ my($k,$blen,$sel,$dz,$dprime) =
&antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <tilt d-prime[m]>',@_);
croak("$0: tilt-dprime outside range [0..$blen]\n")
unless ($dprime>=0 && $dprime<=$blen);
- return ($sel>0) ? T_ravg($kz,$blen) * T_binavg($kz,$sel) * T_binavg($kz,$grez) * T_tilt($kz,$dprime)
- : T_ravg($kz,$blen) * T_binavg($kz,$grez) * T_tilt($kz,$dprime);
+ return ($sel>0) ? T_ravg($k,$blen) * T_binavg($k,$sel) * T_binavg($k,$dz) * T_tilt($k,$dprime)
+ : T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,$dprime);
+}
+
+#----------------------------------------------------------------------
+# T_w(k,blen,dz,range_max)
+# - vertical-velocity method of Thurnherr (IEEE 2011)
+# - range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+{ 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,@_);
+ return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
+ }
}
1;
--- 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 Apr 16 15:58:47 2010
+# dlm: Fri Sep 7 11:11:09 2012
# (c) 1999 A.M. Thurnherr
-# uE-Info: 269 45 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 286 38 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -13,6 +13,7 @@
# 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()
@@ -276,5 +277,15 @@
}
#----------------------------------------------------------------------
+# inverse hyperbolic cosine; mathworld
+# - requires argument >= 1
+#----------------------------------------------------------------------
+
+sub acosh($)
+{
+ return log($_[0] + sqrt($_[0]**2-1));
+}
+
+#----------------------------------------------------------------------
1;
--- a/libstats.pl
+++ b/libstats.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B S T A T S . P L
# doc: Wed Mar 24 13:59:27 1999
-# dlm: Thu Apr 26 09:49:47 2012
+# dlm: Mon Oct 15 10:34:21 2012
# (c) 1999 A.M. Thurnherr
-# uE-Info: 33 64 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 90 30 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -31,6 +31,7 @@
# Mar 10, 2012: - medianF() -> medianAnts_(); mad2F() -> mad2Ants_()
# - added sum()
# Apr 26, 2012: - BUG: std() did not allow nan as stddev input
+# Oct 15, 2012: - added max_i(), min_i()
require "$ANTS/libfuns.pl";
@@ -60,6 +61,17 @@
return $min<9e99 ? $min : nan;
}
+sub min_i(@)
+{
+ my($min) = 9e99;
+ my($min_i);
+
+ for (my($i)=0; $i<=$#_; $i++) {
+ $min_i=$i,$min=$_[$i] if (numberp($_[$i]) && $_[$i] < $min);
+ }
+ return $min<9e99 ? $min_i : nan;
+}
+
sub max(@)
{
my($max) = -9e99;
@@ -69,6 +81,17 @@
return $max>-9e99 ? $max : nan;
}
+sub max_i(@)
+{
+ my($max) = -9e99;
+ my($max_i);
+
+ for (my($i)=0; $i<=$#_; $i++) {
+ $max_i=$i,$max=$_[$i] if (numberp($_[$i]) && $_[$i] > $max);
+ }
+ return $max>-9e99 ? $max_i : nan;
+}
+
sub N(@)
{
my($N) = 0;
--- a/libvec.pl
+++ b/libvec.pl
@@ -1,9 +1,9 @@
#======================================================================
-# . . / L I B / L I B V E C . P L
+# . . / A N T S L I B / L I B V E C . P L
# doc: Sat Mar 20 12:50:32 1999
-# dlm: Tue Jun 5 08:38:52 2012
+# dlm: Mon Jun 11 11:12:06 2012
# (c) 1999 A.M. Thurnherr
-# uE-Info: 318 0 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 35 69 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -32,6 +32,7 @@
# Nov 5, 2009: - added angle(); vel_bias() => angle_diff()
# Apr 22, 2010: - added angle_ts()
# Jun 5, 2012: - added &closestPointOnStraightLine()
+# Jun 11, 2012: - addeed $t output to &closestPointOnStraightLine()
require "$ANTS/libPOSIX.pl"; # acos()
@@ -299,9 +300,11 @@
}
#----------------------------------------------------------------------
-# &closestPointOnStraightLine(lat,lon,lat1,lonA,lat2,lon2)
+# ($lat,$lon,$t) = &closestPointOnStraightLine(lat,lon,lat1,lonA,lat2,lon2)
# - determine point on line segment from <lat1,lonA> to <lat2,lon2> that
# is closest to target point <lat,lon>
+# - $t [0-1] output indicates where along the line segment the closest
+# point lies
# - http://stackoverflow.com/questions/3120357/get-closest-point-to-a-line
# - NOT DONE IN PLANAR GEOMETRY => USE ONLY IN SMALL DOMAINS
#----------------------------------------------------------------------
@@ -315,7 +318,7 @@
my($APlon) = $lonP - $lonA; my($APlat) = $latP - $latA;
my($t) = ($APlon*$ABlon + $APlat*$ABlat) / ($ABlon**2 + $ABlat**2);
return (undef,undef) unless ($t>=0 && $t<=1);
- return ($latA + $t*$ABlat,$lonA + $t*$ABlon);
+ return ($latA + $t*$ABlat,$lonA + $t*$ABlon, $t);
}
1;