.
new file mode 100644
--- /dev/null
+++ b/.interp.ADCP
@@ -0,0 +1,90 @@
+#======================================================================
+# . I N T E R P . A D C P
+# doc: Fri Apr 16 16:07:48 2010
+# dlm: Tue Aug 9 23:07:35 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 81 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# interpolation scheme mimicking RDI ADCP response as described in
+# Broadband primer, p.17
+
+# HISTORY:
+# Apr 16, 2010: - created
+# Aug 9, 2011: - added -u
+
+# NOTES:
+# - interface is described in [.interp.linear]
+
+$IS_opts = 'b:u';
+$IS_optsUsage = '[pass -u)nfiltered] -b)in/pulse <length>';
+
+sub IS_usage()
+{
+ &antsUsageError()
+ unless defined($opt_b);
+}
+
+sub IS_init($$$$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_interpolate(br,idr,xf,xv,xi,f) interpolate field f
+# br data buffer reference
+# idr init-data reference
+# xf x field
+# xv x value
+# xi index of last record with x-value <= x
+# f field number to interpolate
+# <ret val> interpolated value
+#
+# NB:
+# - handle f == xf
+# - return NaN if any of the y values required is NaN
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+ return $xv if ($xf == $f); # return target x
+
+ my($tow) = $xi; # top of triangular sampling window
+ while ($tow>=0 && (!numberp($bR->[$tow][$xf]) || $bR->[$tow][$xf] > $xv-$opt_b)) {
+ $tow--;
+ }
+ if ($tow < 0) { # incomplete window
+ return nan unless ($opt_u);
+ $tow = 0;
+ } else {
+ $tow++;
+ }
+
+ my($bow) = $xi+1; # bottom of triangular sampling window
+ while ($bow<=$#{$bR} && (!numberp($bR->[$bow][$xf]) || $bR->[$bow][$xf] < $xv+$opt_b)) {
+ $bow++;
+ }
+
+ if ($bow > $#{$bR}) { # incomplete window
+ return nan unless ($opt_u);
+ $bow = $#{$bR};
+ } else {
+ $bow--;
+ }
+
+ my($sweight) = 0; # calculate weighted average
+ my($sum) = 0;
+ for (my($i)=$tow; $i<=$bow; $i++) {
+ next unless (numberp($bR->[$i][$xf]) && numberp($bR->[$i][$f]));
+ my($weight) = 1 - abs($bR->[$i][$xf]-$xv)/$opt_b;
+ $sum += $weight * $bR->[$i][$f];
+ $sweight += $weight;
+ }
+
+ return ($sweight>0) ? $sum/$sweight : nan;
+}
+
+#----------------------------------------------------------------------
+
+1;
new file mode 100644
--- /dev/null
+++ b/.interp.fillnan
@@ -0,0 +1,37 @@
+#======================================================================
+# . I N T E R P . F I L L N A N
+# doc: Tue Sep 18 16:36:28 2012
+# dlm: Tue Sep 18 16:43:36 2012
+# (c) 2012 A.M. Thurnherr
+# uE-Info: 34 5 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# fill missing with nans
+
+# HISTORY:
+# Sep 18, 2012: - adapted from [.interp.nnbr]
+
+# see [.interp.linear] for documentation of interface
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+sub IS_usage() {}
+sub IS_init($$$$) {}
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+ if ($f == $xf) {
+ return $xv;
+ } elsif ($bR->[$xi][$xf] == $xv) {
+ return $bR->[$xi][$f];
+ } elsif ($bR->[$xi+1][$xf] == $xv) {
+ return $bR->[$xi+1][$f];
+ } else {
+ return nan;
+ }
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.interp.linear
@@ -0,0 +1,87 @@
+#======================================================================
+# . I N T E R P . L I N E A R
+# doc: Wed Nov 22 21:01:09 2000
+# dlm: Tue Aug 5 14:13:14 2008
+# (c) 2000 A.M. Thurnherr
+# uE-Info: 23 47 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# linear interpolation
+
+# HISTORY:
+# Nov 22, 2000: - created
+# Jul 28, 2006: - Version 3.3 [HISTORY]
+# - added xf to ISInit() args
+# Aug 22, 2006: - modified to allow use with [match]
+# Aug 5, 2008: - added idr param to IS_init()
+
+# NOTES:
+# - the [.interp.*] routines assume that x increases strictly monotonically;
+# for utilities dealing with non-monotonic data, use [.nminterp.*] instead
+# - the [.interp.*] routines are written to work on multiple buffers,
+# rather than just @ants_; this implies that data created by IS_init()
+# must be stored separately for each buffer
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $IS_opts string of allowed options
+# $IS_optsUsage usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_usage() mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_usage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_init(br,f,xf) init interpolation of field f
+# br data buffer reference
+# idr init-data reference
+# f field number
+# xf x field number
+# <ret val> none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_init($$$$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_interpolate(br,idr,xf,xv,xi,f) interpolate field f
+# br data buffer reference
+# idr init-data reference
+# xf x field
+# xv x value
+# xi index of last record with x-value <= x
+# f field number to interpolate
+# <ret val> interpolated value
+#
+# NB:
+# - handle f == xf
+# - return NaN if any of the y values required is NaN
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+ return $xv if ($xf == $f);
+ return nan unless (numberp($bR->[$xi][$f]) && numberp($bR->[$xi+1][$f]));
+
+ my($sc) = ($xv - $bR->[$xi][$xf]) / ($bR->[$xi+1][$xf] - $bR->[$xi][$xf]);
+ return $bR->[$xi][$f] + $sc * ($bR->[$xi+1][$f] - $bR->[$xi][$f]);
+}
+
+#----------------------------------------------------------------------
+
+1;
new file mode 100644
--- /dev/null
+++ b/.interp.linear.orig
@@ -0,0 +1,90 @@
+#======================================================================
+# . I N T E R P . L I N E A R
+# doc: Wed Nov 22 21:01:09 2000
+# dlm: Fri Sep 23 16:32:20 2011
+# (c) 2000 A.M. Thurnherr
+# uE-Info: 62 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# linear interpolation
+
+# HISTORY:
+# Nov 22, 2000: - created
+# Jul 28, 2006: - Version 3.3 [HISTORY]
+# - added xf to ISInit() args
+# Aug 22, 2006: - modified to allow use with [match]
+# Aug 5, 2008: - added idr param to IS_init()
+# Sep 23, 2011: - added support for xfnr==-1 (%RECNO)
+
+# NOTES:
+# - the [.interp.*] routines assume that x increases strictly monotonically;
+# for utilities dealing with non-monotonic data, use [.nminterp.*] instead
+# - the [.interp.*] routines are written to work on multiple buffers,
+# rather than just @ants_; this implies that data created by IS_init()
+# must be stored separately for each buffer
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $IS_opts string of allowed options
+# $IS_optsUsage usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_usage() mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_usage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_init(br,f,xf) init interpolation of field f
+# br data buffer reference
+# idr init-data reference
+# f field number
+# xf x field number or -1 for %RECNO
+# <ret val> none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_init($$$$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &IS_interpolate(bR,idR,xf,xv,xi,f) interpolate field f
+# bR data buffer reference
+# idR init-data reference
+# xf x field or -1 for %RECNO
+# xv x value
+# xi index of last record with x-value <= x
+# f field number to interpolate
+# <ret val> interpolated value
+#
+# NB:
+# - handle f == xf
+# - return NaN if any of the y values required is NaN
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+ return $xv if ($xf == $f);
+ return nan unless (numberp($bR->[$xi][$f]) && numberp($bR->[$xi+1][$f]));
+
+ my($sc) = ($xf < 0)
+ ? ($xv - $bR->[$xi][$xf]) / ($bR->[$xi+1][$xf] - $bR->[$xi][$xf]);
+ : ($xv - $bR->[$xi][$xf]) / ($bR->[$xi+1][$xf] - $bR->[$xi][$xf]);
+ return $bR->[$xi][$f] + $sc * ($bR->[$xi+1][$f] - $bR->[$xi][$f]);
+}
+
+#----------------------------------------------------------------------
+
+1;
new file mode 100644
--- /dev/null
+++ b/.interp.nnbr
@@ -0,0 +1,37 @@
+#======================================================================
+# . I N T E R P . N N B R
+# doc: Wed Nov 22 21:01:09 2000
+# dlm: Wed Aug 3 11:48:59 2011
+# (c) 2000 A.M. Thurnherr
+# uE-Info: 32 0 NIL 0 0 72 8 2 4 NIL ofnI
+#======================================================================
+
+# nearest neighbor resampling
+
+# HISTORY:
+# Nov 22, 2000: - adapted from [.interp.linear]
+# Apr 3, 2004: - added nominal -x
+# Jan 12, 2006: - renamed from [.interp.subsample]
+# Jul 28, 2006: - added xf to ISInit() args
+# Aug 22, 2006: - adapted to work with [match]
+# Aug 5, 2008: - added idr param to IS_init()
+# Aug 3, 2011: - removed -x, which is now handled by [resample]
+
+# see [.interp.linear] for documentation of interface
+
+$IS_opts = "";
+$IS_optsUsage = "";
+
+sub IS_usage() {}
+sub IS_init($$$$) {}
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+ return $bR->[$xi+1][$f]
+ if ($bR->[$xi+1][$xf] - $xv < $xv - $bR->[$xi][$xf]);
+ return $bR->[$xi][$f];
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.interp.poly
@@ -0,0 +1,94 @@
+#======================================================================
+# . I N T E R P . P O L Y
+# doc: Thu Nov 23 21:30:25 2000
+# dlm: Tue Aug 5 14:20:19 2008
+# (c) 2000 A.M. Thurnherr
+# uE-Info: 58 0 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# polynomial interpolation
+
+# HISTORY:
+# Nov 23, 2000: - created
+# Nov 24, 2000: - added -c)enter x-val in window
+# Dec 13, 2002: - replaced -n by -o
+# Oct 22, 2003: - changed -e to -r (conflicts with -e option of [resample])
+# Jul 28, 2006: - removed debugging code
+# - added xf to ISInit() args
+# - BUG: sufficient-input data test was buggy
+# Aug 22, 2006: - adapted to work with [match]
+# Sep 11, 2006: - BUG: sufficient-input ermesg was wrong
+# Aug 5, 2008: - added idr param to IS_init()
+
+# see [.interp.linear] for documentation of interface
+
+# NOTES:
+# - $opt_o <poly order> is required
+# - $opt_r is used to return the error estimates instead of the
+# interpolated values
+# - -o 0 is almost like using -s subsample, the only difference
+# being that subsample returns the real x values while -o 0
+# returns the subsampled x values
+# - on -c the window for each interpolation estimate is centered
+# around the interpolated x value; this implies window-sliding
+# between tabulated x values resulting in functional disconti-
+# nuities
+
+require "$ANTS/polint.pl";
+
+$IS_opts = "co:r";
+$IS_optsUsage = "[e-r)ror estimates] [-c)enter x-val] -o <poly order>";
+
+sub IS_usage()
+{
+ die("$0 (.interp.poly): ERROR! -o required\n")
+ unless (defined($opt_o) && $opt_o >= 0);
+ &antsInfo("WARNING: -c generates discontinuous function")
+ if ($opt_c);
+ &antsInfo("WARNING: high-order polynomial interpolation deprecated")
+ if ($opt_o > 5);
+ &antsInfo("order-1 polynomial; -s linear is more efficient")
+ if ($opt_o == 1);
+ &antsInfo("giving error estimates") if ($opt_r);
+}
+
+sub IS_init($$$$)
+{
+ my($bR,$idR,$f,$xf) = @_;
+ die(sprintf("$0 (.interp.poly): ERROR! " .
+ "not enough data points (%d) for order-$opt_o poly\n",
+ $#{$bR}+1))
+ if ($opt_o > $#{$bR});
+}
+
+sub findWindow($$$$$) # find window of size n=opt_o+1
+{
+ my($bR,$f,$v,$i,$n) = @_;
+ my($cn); # current win size
+
+ return ($bR->[$i+1][$f] - $v < $v - $bR->[$i][$f]) ? $i+1 : $i
+ if ($n == 1); # nearest neighbor
+
+ if ($opt_c) { # center window around x value
+ for ($cn=2; $cn<$n; $cn++) { # grow window
+ $i-- if (($i>0 && $bR->[$i+$cn-1][$f] - $v >= $v - $bR->[$i][$f]) ||
+ ($i+$cn-1 == $#{$bR}));
+ }
+ } else {
+ $i -= int(($n-1)/2); # ideally
+ $i = 0 if ($i < 0); # off leftmost column
+ $i = $#{$bR}+1 - $n if ($i+$n > $#{$bR}+1);# off rightmost column
+ }
+
+ return $i;
+}
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+ return $xv if ($xf == $f);
+
+ $xi = &findWindow($bR,$xf,$xv,$xi,$opt_o+1);
+ my($y,$dy) = &polint($bR,$xf,$xv,$xi,$opt_o+1,$f);
+ return $opt_r ? $dy : $y;
+}
new file mode 100644
--- /dev/null
+++ b/.interp.spline
@@ -0,0 +1,125 @@
+#======================================================================
+# . I N T E R P . S P L I N E
+# doc: Wed Nov 22 21:01:09 2000
+# dlm: Tue Aug 5 14:20:39 2008
+# (c) 2000 A.M. Thurnherr
+# uE-Info: 31 0 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# spline interpolation
+
+# HISTORY:
+# May 27, 2003: - adapted from [.interp.linear]
+# Jul 19, 2004: - BUG: i instead of $i --- 'orrible!
+# May 15, 2006: - BUG: assertion did not work on -e
+# Jul 28, 2006: - Version 3.3 [HISTORY]
+# - made 100% compatible with [.interp.spline]
+# - added xf to ISInit() args
+# Jul 30, 2006: - continued debugging
+# Aug 22, 2006: - adapted to work with [match]
+# Aug 5, 2008: - added idr param to IS_init()
+
+# see [.interp.linear] for documentation of interface
+
+# EXAMPLE:
+# Matlab: x = 0:10; y = sin(x); xx = 0:.25:10;
+# yy = spline(x,y,xx); plot(x,y,'o',xx,yy)
+# gnuplot/ANTS:
+# set style function points;
+# set style data lines;
+# plot [0:10] sin(x), \
+# '<Cat -f x:0,1,10 | list x y=sin($x) | resample -s spline x \#0-10:0.25'
+
+$IS_opts = "B:";
+$IS_optsUsage = "[-B)oundary cond <1st/last>]";
+
+#----------------------------------------------------------------------
+
+my($yp1,$ypn);
+
+sub IS_usage()
+{
+ if (defined($opt_B)) { # NR p.115
+ ($yp1,$ypn) = split('/',$opt_B);
+ croak("$0: can't decode -B $opt_B\n")
+ unless (numberp($yp1) && numberp($ypn));
+ }
+}
+
+#----------------------------------------------------------------------
+
+sub IS_init($$$$) {
+ my($bR,$idR,$f,$xf) = @_;
+ my($i,$k,$p,$qn,$sig,$un,@u);
+ my($n) = scalar(@{$bR});
+
+ if (defined($opt_B)) { # handle boundary cond
+ $idR->[1][$f] = -0.5;
+ $u[1] = (3/($bR->[1][$xf]-$bR->[0][$xf])) *
+ (($bR->[1][$f]-$bR->[0][$f]) /
+ ($bR->[1][$xf]-$bR->[0][$xf]) - $yp1);
+ } else {
+ $idR->[1][$f] = $u[1] = 0;
+ }
+
+ for ($i=2; $i<=$n-1; $i++) {
+ $sig = ($bR->[$i-1][$xf]-$bR->[$i-2][$xf]) /
+ ($bR->[$i][$xf]-$bR->[$i-2][$xf]);
+ $p = $sig*$idR->[$i-1][$f] + 2;
+ $idR->[$i][$f] = ($sig-1)/$p;
+ $u[$i] = ($bR->[$i][$f]-$bR->[$i-1][$f]) /
+ ($bR->[$i][$xf]-$bR->[$i-1][$xf])
+ - ($bR->[$i-1][$f]-$bR->[$i-2][$f]) /
+ ($bR->[$i-1][$xf]-$bR->[$i-2][$xf]);
+ $u[$i] = (6*$u[$i]/($bR->[$i][$xf]-$bR->[$i-2][$xf]) -
+ $sig * $u[$i-1]) / $p;
+ }
+
+ if (defined($opt_B)) {
+ $qn = 0.5;
+ $un = (3/($bR->[$n-1][$xf]-$bR->[$n-2][$xf])) *
+ ($ypn-($bR->[$n-1][$f]-$bR->[$n-2][$f]) /
+ ($bR->[$n-1][$xf]-$bR->[$n-2][$xf]));
+ } else {
+ $qn = $un = 0;
+ }
+
+ $idR->[$n][$f] = ($un-$qn*$u[$n-1]) / ($qn*$idR->[$n-1][$f]+1);
+ for ($k=$n-1;$k>=1;$k--) {
+ $idR->[$k][$f] = $idR->[$k][$f] * $idR->[$k+1][$f] + $u[$k];
+ }
+}
+
+#----------------------------------------------------------------------
+
+sub IS_interpolate($$$$$$)
+{
+ my($bR,$idR,$xf,$xv,$xi,$f) = @_;
+
+ return $xv if ($xf == $f);
+
+ return $bR->[$xi][$f] # edge values are ok
+ if equal($xv,$bR->[$xi][$xf]);
+ return $bR->[$xi+1][$f]
+ if equal($xv,$bR->[$xi+1][$xf]);
+
+ return nan unless (numberp($bR->[$xi][$f]) &&
+ numberp($bR->[$xi+1][$f]));
+
+ croak("$0: assertion $bR->[$xi+1][$xf] >= $xv >= $bR->[$xi][$xf] failed")
+ unless (defined($opt_e) || $bR->[$xi+1][$xf] >= $xv && $xv >= $bR->[$xi][$xf]);
+
+ my($h) = $bR->[$xi+1][$xf] - $bR->[$xi][$xf];
+ croak("$0: assertion #2 failed") unless ($h > 0);
+ my($a) = ($bR->[$xi+1][$xf] - $xv) / $h;
+ my($b) = ($xv - $bR->[$xi][$xf]) / $h;
+
+ return $a*$bR->[$xi][$f] + $b*$bR->[$xi+1][$f] +
+ (($a*$a*$a-$a) * $idR->[$xi+1][$f] +
+ ($b*$b*$b-$b) * $idR->[$xi+2][$f]) *
+ ($h*$h)/6;
+}
+
+#----------------------------------------------------------------------
+
+1;
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.gamma_n
@@ -0,0 +1,36 @@
+#======================================================================
+# . I S O P Y C N A L _ T S . G A M M A _ N
+# doc: Tue Dec 13 21:50:18 2005
+# dlm: Mon Dec 19 13:04:15 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 27 30 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 14, 2005: - created
+# Dec 19, 2005: - finalized
+
+# NOTES:
+# - requires %lat/%lon PARAMs
+
+require "$ANTS/libgamma.pl"; # load equation of state
+
+unless (defined($P{ITS})) {
+ &antsInfo("using default %ITS=90");
+ &antsAddParams(ITS,90);
+}
+$gamma::temp_scale = $P{ITS};
+
+croak("$0: need %lat/%lon\n")
+ unless defined($P{lat}) && defined($P{lon});
+croak("$0: need %press\n")
+ unless defined($P{press});
+
+sub density($$)
+{
+ my($S,$T) = @_;
+ my($gamma) = gamma::gamma_n($S,$T,$P{press},$P{lat},$P{lon});
+ return $gamma > 0 ? $gamma : undef;
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.sigma0
@@ -0,0 +1,29 @@
+#======================================================================
+# . I S O P Y C N A L _ T S . S I G M A 0
+# doc: Mon Dec 19 12:28:01 2005
+# dlm: Mon Dec 19 13:27:24 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 19 56 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+unless (defined($P{ITS})) {
+ &antsInfo("using default %ITS=90");
+ &antsAddParams(ITS,90);
+}
+
+&antsInfo("Warning: ignoring non-zero %press = $P{press}")
+ if ($P{press} != 0);
+$P{press} = 0;
+
+sub density($$)
+{
+ my($S,$T) = @_;
+ return sigma($S,$T,$P{press},0);
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.sigma2
@@ -0,0 +1,29 @@
+#======================================================================
+# . I S O P Y C N A L _ T S . S I G M A 2
+# doc: Mon Dec 19 12:28:01 2005
+# dlm: Mon Dec 19 13:29:23 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 26 32 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+unless (defined($P{ITS})) {
+ &antsInfo("using default %ITS=90");
+ &antsAddParams(ITS,90);
+}
+
+&antsInfo("Warning: ignoring non-zero %press = $P{press}")
+ if ($P{press} != 0);
+$P{press} = 0;
+
+sub density($$)
+{
+ my($S,$T) = @_;
+ return sigma($S,$T,$P{press},2000);
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.sigma4
@@ -0,0 +1,29 @@
+#======================================================================
+# . I S O P Y C N A L _ T S . S I G M A 4
+# doc: Mon Dec 19 12:28:01 2005
+# dlm: Wed Sep 16 10:19:10 2009
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 26 34 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Sep 16, 2009: created
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+unless (defined($P{ITS})) {
+ &antsInfo("using default %ITS=90");
+ &antsAddParams(ITS,90);
+}
+
+&antsInfo("Warning: ignoring non-zero %press = $P{press}")
+ if ($P{press} != 0);
+$P{press} = 0;
+
+sub density($$)
+{
+ my($S,$T) = @_;
+ return sigma($S,$T,$P{press},4000);
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.theta0
@@ -0,0 +1,22 @@
+#======================================================================
+# . I S O P Y C N A L _ T S . T H E T A 0
+# doc: Mon Dec 19 09:43:02 2005
+# dlm: Mon Dec 19 13:38:50 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 21 42 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# calculate in-situ temperature at surface given theta2
+
+# HISTORY:
+# Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+$temperature_fname = 'theta0';
+
+sub temperature($$)
+{
+ my($salin,$potemp) = @_;
+ return &temp($salin,$potemp,$P{press},0);
+}
new file mode 100644
--- /dev/null
+++ b/.isopycnal_TS.theta2
@@ -0,0 +1,22 @@
+#======================================================================
+# . I S O P Y C N A L _ T S . T H E T A 2
+# doc: Mon Dec 19 09:43:02 2005
+# dlm: Mon Dec 19 13:43:18 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 21 32 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# calculate in-situ temperature at surface given theta2
+
+# HISTORY:
+# Dec 19, 2005: - created
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+$temperature_fname = 'theta2';
+
+sub temperature($$)
+{
+ my($salin,$potemp) = @_;
+ return &temp($salin,$potemp,$P{press},2000);
+}
new file mode 100644
--- /dev/null
+++ b/.lmfit.exp
@@ -0,0 +1,117 @@
+#======================================================================
+# . L M F I T . E X P
+# doc: Wed Feb 24 09:40:06 1999
+# dlm: Fri Jul 28 13:40:56 2006
+# (c) 1999 A.M. Thurnherr
+# uE-Info: 30 41 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+# - a number of global variables to be set during loading
+# - a number of subs to perform admin tasks (usage, init, ...)
+# - a sub to evaluate the model function which is to be fitted using
+# a number of pararams which are all stored in @A (beginning at
+# A[1]!!!). You also need to return the partial derivatives of
+# the model function wrt all params.
+# - the interface is documented between +++++++ lines
+
+# fit exponential A[3]+A[2]*exp(A[1]*x) to data
+#
+# NOTES:
+# - initial parameter estimates are crucial
+# - there is currently no heuristics
+
+# HISTORY:
+# Mar 11, 1999: - created from [./.mfit.poly] & [./.mfit.gauss]
+# Jul 31, 1999: - typecheck usage
+# Mar 17, 2001: - param->arg
+# Jan 16, 2006: - added notes
+# Jul 28, 2006: - Version 3.3 [HISTORY]
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $modelOpts string of allowed options
+# $modelOptsUsage usage information string for options
+# $modelMinArgs min # of arguments of model
+# $modelArgsUsage usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+# $modelNFit number of params to fit in model
+# @nameA symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "";
+$modelOptsUsage = "";
+$modelMinArgs = 0;
+$modelArgsUsage = "[exp [mul [add guess]]]";
+$modelNFit = 3;
+$nameA[1] = "exp";
+$nameA[2] = "mul";
+$nameA[3] = "add";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage() mangle parameters; NB: there may be `infinite' # of
+# filenames after model arguments; this usually sets
+# @A (the model parameters) but these can later be
+# calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+ $A[1] = nan; $A[2] = nan; $A[3] = nan; # usage
+ $A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ $A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ $A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit() initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+ $A[1] = 1 unless (numberp($A[1]));
+ $A[2] = 1 unless (numberp($A[2]));
+ $A[3] = 0 unless (numberp($A[3]));
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda) evaluate polynom and derivatives
+# x x value (NOT xfnr)
+# A reference to @A
+# dyda reference to array for partial derivatives
+# (wrt individaul params in @A)
+# <ret val> y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$) # y = A[3]+A[2]*exp(A[1]*x)
+{
+ my($x,$AR,$dydaR) = @_;
+ my($e) = exp($AR->[1]*$x);
+
+ $dydaR->[1] = $AR->[2]*$x*$e;
+ $dydaR->[2] = $e;
+ $dydaR->[3] = 1; # partial derivatives
+
+ return $AR->[3] + $AR->[2]*$e;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup() cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+}
new file mode 100644
--- /dev/null
+++ b/.lmfit.gauss
@@ -0,0 +1,176 @@
+#======================================================================
+# . L M F I T . G A U S S
+# doc: Wed Feb 24 09:40:06 1999
+# dlm: Fri Jul 28 13:32:35 2006
+# (c) 1999 A.M. Thurnherr
+# uE-Info: 35 51 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+# - a number of global variables to be set during loading
+# - a number of subs to perform admin tasks (usage, init, ...)
+# - a sub to evaluate the model function which is to be fitted using
+# a number of pararams which are all stored in @A (beginning at
+# A[1]!!!). You also need to return the partial derivatives of
+# the model function wrt all params.
+# - the interface is documented between +++++++ lines
+
+# Gauss data model (i.e. fit Gaussian curve)
+# NB: - fitting is rather sensitive to the input parameters, thus
+# a heuristic has been added to guess them (by setting them
+# to NaN)
+# - another fickle parameter is the y-offset (zero line); thus
+# a heuristics has been added for this one as well
+# - the parameters are peak, mean, standard deviation
+
+# HISTORY:
+# Feb 24, 1999: - created together with [./cfit]
+# Feb 25, 1999: - cosmetic changes
+# Jul 31, 1999: - parameter typecheck
+# Oct 04, 1999: - changed param names
+# Oct 05, 1999: - improved heuristics
+# - changed e-scale to sigma
+# Mar 17, 2001: - param->arg
+# Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan()
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $modelOpts string of allowed options
+# $modelOptsUsage usage information string for options
+# $modelMinArgs min # of arguments of model
+# $modelArgsUsage usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+# $modelNFit number of params to fit in model
+# @nameA symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "y";
+$modelOptsUsage = "[-y)shift]";
+$modelMinArgs = 0;
+$modelArgsUsage = "[peak guess [mean guess [sigma guess]]]";
+$modelNFit = 3;
+$nameA[1] = "peak";
+$nameA[2] = "mean";
+$nameA[3] = "sigma";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage() mangle parameters; NB: there may be `infinite' # of
+# filenames after model arguments; this usually sets
+# @A (the model parameters) but these can later be
+# calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+ $A[1] = nan; $A[2] = nan; $A[3] = nan; # usage
+ $A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ $A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ $A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit() initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+ my($i,$j,$ymin,$ymax,$xatymax);
+
+# --------------------------------------------------
+# heuristics for initial model param values
+# --------------------------------------------------
+
+ $ymin = 1e33, $ymax = -1e33, $xatymax = 0;
+ for ($i=0; $i<=$#ants_; $i++) {
+ next if ($antsFlagged[$i]);
+ $ymin = $ants_[$i][$yfnr]
+ if ($ants_[$i][$yfnr] < $ymin);
+ $ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr]
+ if ($ants_[$i][$yfnr] > $ymax);
+ }
+ $A[1] = $ymax - $ymin if isnan($A[1]); # peak guess
+ $A[2] = $xatymax if isnan($A[2]); # mean guess
+ if (isnan($A[3])) { # sigma guess
+ for ($i=1;
+ $i<=$#ants_ && !$antsFlagged[$i]
+ && $ants_[$i][$yfnr]-$ymin<0.36*$A[1];
+ $i++) {}
+ for ($j=$#ants_;
+ $j>=1 && !$antsFlagged[$i]
+ && $ants_[$j][$yfnr]-$ymin < 0.36*$A[1];
+ $j--) {}
+ $A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0;
+ $A[3] *= 0.71; # scale by 1/sqrt(2)
+ if ($A[3] == 0) {
+ &antsInfo("$model: sigma heuristic failed (set to 1)!");
+ $A[3] = 1;
+ }
+ }
+
+# --------------------------------------------------
+# y shift (-y option)
+# --------------------------------------------------
+
+ if ($opt_y) {
+ for ($i=1; $i<=$#ants_; $i++) {
+ $ants_[$i][$yfnr] -= $ymin;
+ }
+ }
+
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda) evaluate sum of Gaussians (p.528) at x
+# x x value (NOT xfnr)
+# A reference to @A
+# dyda reference to array for partial derivatives
+# (wrt individaul params in @A)
+# <ret val> y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+ my($x,$AR,$dydaR) = @_;
+ my($i,$fac,$ex,$arg,$sqrt2sig);
+ my($y) = 0;
+
+ for ($i=1; $i < $#{$AR}; $i+=3) {
+ $sqrt2sig = (1.4142135623731*$AR->[$i+2]);
+ $arg = ($x - $AR->[$i+1]) / $sqrt2sig;
+ $ex = exp(-$arg*$arg);
+ $fac = $AR->[$i] * $ex * 2*$arg;
+ $y += $AR->[$i] * $ex;
+
+ $dydaR->[$i] = $ex;
+ $dydaR->[$i+1] = $fac / $sqrt2sig;
+ $dydaR->[$i+2] = $fac * $arg / $sqrt2sig;
+ }
+ return $y;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup() cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+ if ($opt_y) {
+ $A[1] += $ymin;
+ for ($i=1; $i<=$#ants_; $i++) {
+ $ants_[$i][$yfnr] += $ymin;
+ }
+ }
+}
new file mode 100644
--- /dev/null
+++ b/.lmfit.normal
@@ -0,0 +1,175 @@
+#======================================================================
+# . L M F I T . N O R M A L
+# doc: Wed Feb 24 09:40:06 1999
+# dlm: Fri Jul 28 13:35:24 2006
+# (c) 1999 A.M. Thurnherr
+# uE-Info: 34 51 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+# - a number of global variables to be set during loading
+# - a number of subs to perform admin tasks (usage, init, ...)
+# - a sub to evaluate the model function which is to be fitted using
+# a number of pararams which are all stored in @A (beginning at
+# A[1]!!!). You also need to return the partial derivatives of
+# the model function wrt all params.
+# - the interface is documented between +++++++ lines
+
+# check if a given distribution is normal
+# NB:
+# - fitting is based on gauss curve fitting [.lmfit.gauss]
+# - heuristics are taken from there and scaled for the normal
+# parameter choices
+# - simplified, e.g. y-shift is removed (does not make sense for
+# distribution)
+# - added chi^2 significance testing to &modelCleanup() on -x
+
+# HISTORY:
+# Oct 04, 1999: - created from [.lmfit.gauss]
+# Oct 05, 1999: - added chi^2 significance test
+# - removed -y
+# - improved heuristics
+# Mar 17, 2001: - param->arg
+# Jul 28, 2006: - Version 3.3 [HISTORY]; &isnan()
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $modelOpts string of allowed options
+# $modelOptsUsage usage information string for options
+# $modelMinArgs min # of arguments of model
+# $modelArgsUsage usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+# $modelNFit number of params to fit in model
+# @nameA symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "x";
+$modelOptsUsage = "[-x chi^2 test]";
+$modelMinArgs = 0;
+$modelArgsUsage = "[area guess [mean guess [sigma guess]]]";
+$modelNFit = 3;
+$nameA[1] = "area";
+$nameA[2] = "mean";
+$nameA[3] = "sigma";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage() mangle parameters; NB: there may be `infinite' # of
+# filenames after model arguments; this usually sets
+# @A (the model parameters) but these can later be
+# calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+ $A[1] = nan; $A[2] = nan; $A[3] = nan; # usage
+ $A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ $A[2] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ $A[3] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit() initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+ my($i,$j,$ymin,$ymax,$xatymax);
+
+# --------------------------------------------------
+# heuristics for initial model param values
+# --------------------------------------------------
+
+ $ymin = 1e33, $ymax = -1e33, $xatymax = 0;
+ for ($i=0; $i<=$#ants_; $i++) {
+ next if ($antsFlagged[$i]);
+ $ymin = $ants_[$i][$yfnr]
+ if ($ants_[$i][$yfnr] < $ymin);
+ $ymax = $ants_[$i][$yfnr], $xatymax = $ants_[$i][$xfnr]
+ if ($ants_[$i][$yfnr] > $ymax);
+ }
+ $A[1] = $ymax - $ymin if isnan($A[1]); # peak guess
+ $A[2] = $xatymax if isnan($A[2]); # mean guess
+ if (isnan($A[3])) { # e-scale guess
+ for ($i=1;
+ $i<=$#ants_ && !$antsFlagged[$i]
+ && $ants_[$i][$yfnr]-$ymin<0.36*$A[1];
+ $i++) {}
+ for ($j=$#ants_;
+ $j>=1 && !$antsFlagged[$i]
+ && $ants_[$j][$yfnr]-$ymin < 0.36*$A[1];
+ $j--) {}
+ $A[3] = abs($ants_[$i][$xfnr]-$ants_[$j][$xfnr]) / 2.0;
+ if ($A[3] == 0.0) {
+ &antsInfo("$model: sigma heuristic failed (set to 1)!");
+ $A[3] = 1.0;
+ }
+ }
+
+ $A[1] *= 1.77 * $A[3] # gauss -> normal
+ unless (isnan($A[1]) || isnan($A[3]));
+ $A[3] *= 0.71 unless isnan($A[3]);
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda) evaluate Normal distribution curve at x
+# x x value (NOT xfnr)
+# A reference to @A
+# dyda reference to array for partial derivatives
+# (wrt individual params in @A)
+# <ret val> y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+ my($x,$AR,$dydaR) = @_;
+
+ my($peak) = $AR->[1] / (2.506628274631 * $AR->[3]);
+ my($dx ) = $x - $AR->[2];
+ my($sig2) = $AR->[3] * $AR->[3];
+ my($expo) = exp(-$dx*$dx/(2*$sig2));
+ my($norm) = $peak * $expo;
+
+ if (defined($dydaR)) {
+ $dydaR->[1] = $norm / $AR->[1];
+ $dydaR->[2] = $norm * $dx / $sig2;
+ $dydaR->[3] = $norm/$AR->[3] * ($dx*$dx/$sig2 - 1);
+ }
+
+ return $norm;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup() cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+ return unless ($opt_x);
+
+ require "$ANTS/libfuns.pl";
+ my($chisq) = 0;
+ my($nval,$prob,$sign);
+
+ for ($i=0; $i<=$#ants_; $i++) {
+ next if ($antsFlagged[$i]);
+# next if ($ants_[$i][$yfnr] <= 1); # IGNORE TAIL HEURISTICS
+ $nval = &modelEvaluate($ants_[$i][$xfnr],\@A);
+ $chisq += ($ants_[$i][$yfnr] - $nval)**2 / $nval;
+ }
+ $prob = &gammq(($ndata-3)/2,$chisq/2);
+ $sign = int($prob*100);
+ &antsInfo("$model: normal-distr. hypothesis disproved at $sign%% sign. level");
+}
new file mode 100644
--- /dev/null
+++ b/.lmfit.poly
@@ -0,0 +1,129 @@
+#======================================================================
+# . L M F I T . P O L Y
+# doc: Wed Feb 24 09:40:06 1999
+# dlm: Fri Jul 28 13:35:50 2006
+# (c) 1999 A.M. Thurnherr
+# uE-Info: 28 41 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+# - a number of global variables to be set during loading
+# - a number of subs to perform admin tasks (usage, init, ...)
+# - a sub to evaluate the model function which is to be fitted using
+# a number of pararams which are all stored in @A (beginning at
+# A[1]!!!). You also need to return the partial derivatives of
+# the model function wrt all params.
+# - the interface is documented between +++++++ lines
+
+# fit polynomial (sum of A_i x^i) to data
+# NB:
+
+# HISTORY:
+# Feb 25, 1999: - created
+# Mar 14, 1999: - cosmetic changes
+# Jul 31, 1999: - argument typechecking
+# Mar 17, 2001: - param->arg
+# Jan 12, 2006: - specify order with -o as in [.lsfit.poly]
+# Jul 28, 2006: - Version 3.3 [HISTORY]
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $modelOpts string of allowed options
+# $modelOptsUsage usage information string for options
+# $modelMinArgs min # of arguments of model
+# $modelArgsUsage usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+# $modelNFit number of params to fit in model
+# @nameA symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "o:";
+$modelOptsUsage = "-o)rder <n>";
+$modelMinArgs = 0;
+$modelArgsUsage = "[c0 [c1 [...]]]";
+
+&antsInfo("non-linear method deprecated; use `lsfit' instead");
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage() mangle parameters; NB: there may be `infinite' # of
+# filenames after model arguments; this usually sets
+# @A (the model parameters) but these can later be
+# calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+ my($c);
+
+ die("$0 (.lmfit.poly): ERROR! -o required\n") # order of polynomial
+ unless (defined($opt_o) && $opt_o >= 0);
+ $modelNFit = &antsCardOpt($opt_o)+1;
+
+ for ($c=0; $c<$modelNFit; $c++) { # init coefficients
+ if ($#ARGV >= 0 && ! -r $ARGV[0]) {
+ $A[$c+1] = &antsFloatArg();
+ } else {
+ $A[$c+1] = nan;
+ }
+ $nameA[$c+1] = "c$c"; # and names
+ }
+ &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit() initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+ my($c);
+
+ for ($c=0; $c<$modelNFit; $c++) { # init coefficients
+ $A[$c+1] = 10**-$c unless (numberp($A[$c+1]));
+ }
+
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda) evaluate polynomial and derivatives
+# x x value (NOT xfnr)
+# A reference to @A
+# dyda reference to array for partial derivatives
+# (wrt individaul params in @A)
+# <ret val> y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+ my($x,$AR,$dydaR) = @_;
+ my($i);
+ my($pow) = 1;
+ my($y) = 0;
+
+ for ($i=1; $i<=$modelNFit; $i++) {
+ $y += $AR->[$i]*$pow;
+ $dydaR->[$i] = $pow;
+ $pow *= $x;
+ }
+ return $y;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup() cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+}
new file mode 100644
--- /dev/null
+++ b/.lmfit.sqrt
@@ -0,0 +1,105 @@
+#======================================================================
+# . L M F I T . S Q R T
+# doc: Fri Oct 10 15:50:42 2014
+# dlm: Sat Oct 11 09:30:48 2014
+# (c) 2014 A.M. Thurnherr
+# uE-Info: 27 32 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# model function to your data:
+# - a number of global variables to be set during loading
+# - a number of subs to perform admin tasks (usage, init, ...)
+# - a sub to evaluate the model function which is to be fitted using
+# a number of pararams which are all stored in @A (beginning at
+# A[1]!!!). You also need to return the partial derivatives of
+# the model function wrt all params.
+# - the interface is documented between +++++++ lines
+
+# fit square root A[1]*sqrt(x) to data
+#
+# NOTES:
+# - initial parameter estimates may be important
+# - there is currently no heuristics
+
+# HISTORY:
+# Oct 10, 2014: - created from [.lmfit.exp]
+# Oct 11, 2014: - made it work
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $modelOpts string of allowed options
+# $modelOptsUsage usage information string for options
+# $modelMinArgs min # of arguments of model
+# $modelArgsUsage usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+# $modelNFit number of params to fit in model
+# @nameA symbolic names of model parameters
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "";
+$modelOptsUsage = "";
+$modelMinArgs = 0;
+$modelArgsUsage = "[scale guess]";
+$modelNFit = 1;
+$nameA[1] = "scale";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage() mangle parameters; NB: there may be `infinite' # of
+# filenames after model arguments; this usually sets
+# @A (the model parameters) but these can later be
+# calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+ $A[1] = nan;
+ $A[1] = &antsFloatArg() if ($#ARGV >= 0 && ! -r $ARGV[0]);
+ &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit() initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit()
+{
+ $A[1] = 1 unless (numberp($A[1])); # scale
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(x,A,dyda) evaluate polynom and derivatives
+# x x value (NOT xfnr)
+# A reference to @A
+# dyda reference to array for partial derivatives
+# (wrt individaul params in @A)
+# <ret val> y value
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$) # y = A[1]*sqrt(*x)
+{
+ my($x,$AR,$dydaR) = @_;
+ my($v) = sqrt($x);
+
+ $dydaR->[1] = $v; # dy/dA[1] = sqrt(x)
+ return $AR->[1]*$v;
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup() cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+}
new file mode 100644
--- /dev/null
+++ b/.lsfit.bilin
@@ -0,0 +1,141 @@
+#======================================================================
+# . L S F I T . B I L I N
+# doc: Wed Feb 24 09:40:06 1999
+# dlm: Fri Jul 28 13:36:36 2006
+# (c) 1999 A.M. Thurnherr
+# uE-Info: 32 41 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# What you need to provide if you wanna fit a different
+# linear model function to your data:
+# - a number of global variables to be set during loading
+# - a number of subs to perform admin tasks (usage, init, ...)
+# - a sub to evaluate the basis funs at a given x value; each
+# y value must be stored in @A (beginning with
+# A[1]!!!).
+# - the interface is documented between +++++++ lines
+
+# fit bi-linear function to data, i.e. y = A + B*x1 + C*x2
+
+# HISTORY:
+# Aug 01, 1999: - adapted from [.lsfit.poly]
+# Aug 02, 1999: - added &antsDescription()
+# Sep 26, 1999: - cosmetics
+# - added vars & covars
+# Sep 27, 1999: - changed from covar to sigmas
+# Oct 01, 1999: - cosmetics
+# Oct 06, 1999: - added -l
+# Mar 17, 2001: - param->arg
+# May 2, 2001: - updated doc
+# Nov 17, 2005: - commented out antsDescription()
+# - updated stats on -p
+# Jul 28, 2006: - Version 3.3 [HISTORY]
+
+# NOTES:
+# - could be easily extended to multidimensional linear fit but
+# what's the use?
+# - -p zeroes param C for bilinear spice method (T = A + B*sig + C*neph)
+# so that residual field becomes spice anomaly (linearly correlated
+# with neph)
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $modelOpts string of allowed options
+# $modelOptsUsage usage information string for options
+# $modelMinArgs min # of arguments of model
+# $modelArgsUsage usage information string for arguments
+#
+# The following variables may be set later but not after &modelInit()
+#
+# $modelNFit number of params to fit in model
+# @nameA symbolic names of model parameters
+#
+# You should call &antsDescription() for the -ct options here
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$modelOpts = "pl:";
+$modelOptsUsage = "[-p)artial f/r] [-l)imit r_BC <min,max>]";
+$modelMinArgs = 1;
+$modelArgsUsage = "<2nd x-field>";
+$modelNFit = 3;
+$nameA[1] = "A"; $nameA[2] = "B"; $nameA[3] = "C";
+$A[1] = nan; $A[2] = nan; $A[3] = nan;
+#&antsDescription("c","bilin_$nameA[1]",
+# "c","bilin_$nameA[2]",
+# "c","bilin_$nameA[3]");
+#&antsDescription("t","bilin_sigma_$nameA[1]",
+# "t","bilin_sigma_$nameA[2]",
+# "t","bilin_sigma_$nameA[3]",
+# "t","bilin_ccc_$nameA[1]_$nameA[2]",
+# "t","bilin_ccc_$nameA[1]_$nameA[3]",
+# "t","bilin_ccc_$nameA[2]_$nameA[3]");
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelUsage() mangle parameters; NB: there may be `infinite' # of
+# filenames after model arguments; this usually sets
+# @A (the model parameters) but these can later be
+# calculated heuristically during &modelInit()
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelUsage()
+{
+ if (defined($opt_l)) {
+ ($minBC,$maxBC) = split(',',$opt_l);
+ &antsUsageError("\n>>> error with -l")
+ unless (defined($maxBC) && $maxBC > $minBC);
+ }
+ $x2fnr = &antsFieldArg();
+ &antsUsageError() unless ($#ARGV < 0 || -r $ARGV[0]);
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelInit() initializes model after reading of data
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelInit() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &modelEvaluate(idx,xfnr,vals) evaluate basis funs at x (NB: x1, x2)
+# idx current index in @ants_
+# xfnr field number of x field
+# vals reference to return values (1-relative!)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelEvaluate($$$)
+{
+ my($idx,$xfnr,$valsR) = @_;
+ my($i);
+
+ $valsR->[1] = 1;
+ $valsR->[2] = $ants_[$idx][$xfnr];
+ $valsR->[3] = $ants_[$idx][$x2fnr];
+}
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+# &modelCleanup() cleans up after fitting but before output
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub modelCleanup()
+{
+ if (defined($opt_l)) {
+ my($ccc) = $covar[2][3]/(sqrt($covar[2][2])*sqrt($covar[3][3]));
+ if ($ccc < $minBC || $ccc > $maxBC) {
+ &antsInfo("CCC B/C = %.3g out of range, fit discarded",$ccc);
+ $suppressFit = 1;
+ }
+ }
+ if ($opt_p) {
+ &antsInfo("parameter $nameA[3] = %.3g discarded",$A[3]);
+ $A[3] = 'discarded';
+ $covar[3][3] = $RMS = $sig = nan;
+ }
+}
new file mode 100644
--- /dev/null
+++ b/.match_minimize.MAD
@@ -0,0 +1,54 @@
+#======================================================================
+# . M A T C H _ M I N I M I Z E . M A D
+# doc: Tue Aug 22 18:31:22 2006
+# dlm: Wed Aug 23 23:15:13 2006
+# (c) 2006 A.M. Thurnherr
+# uE-Info: 39 26 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# sum of absolute deviations
+
+#======================================================================
+
+sub MM_usage(@)
+{
+ @MM_fname = @_;
+
+ croak("match-minimize usage: -m SAD(<field>[,...])\n")
+ unless (@MM_fname > 0);
+ for (my($fi)=0; $fi<@MM_fname; $fi++) {
+ $MM_fnr[$fi] = &fnr($MM_fname[$fi]);
+ $MM_wfnr[$fi] = &wfnr($MM_fname[$fi]);
+ &IS_init($MM_wfnr[$fi],$MM_cwfnr); # prepare interpolation
+ }
+ for (my($r)=0; $r<@ants_; $r++) {
+ for (my($fi)=0; $fi<@MM_fname; $fi++) {
+ croak("$0: $MM_fname[$fi] must be numeric\n")
+ unless numberp($ants_[$r][$MM_fnr[$fi]]) &&
+ numberp($ants_[$r][$MM_wfnr[$fi]]);
+ }
+ }
+}
+
+#======================================================================
+
+sub MM_eval(@)
+{
+ &MW_warp(@_);
+
+ my($sad) = my($nad) = 0;
+ for (my($r)=0; $r<@ants_; $r++) {
+ my($xi) = &bSearch($MW_cwfnr,$ants_[$r][$MW_cfnr]);
+ next unless defined($xi);
+ $nad++;
+ for (my($fi)=0; $fi<@MM_fnr; $fi++) {
+ my($wv) = &IS_interpolate(\@wf_,$MW_cwfnr,$ants_[$r][$MW_cfnr],$xi,$MM_wfnr[$fi]);
+ $sad += abs($ants_[$r][$MM_fnr[$fi]] - $wv);
+ }
+ }
+ return $sad/$nad;
+}
+
+#======================================================================
+
+1;
new file mode 100644
--- /dev/null
+++ b/.match_warp.stretch
@@ -0,0 +1,59 @@
+#======================================================================
+# . M A T C H _ W A R P . S T R E T C H
+# doc: Tue Aug 22 18:31:22 2006
+# dlm: Wed Aug 23 23:51:25 2006
+# (c) 2006 A.M. Thurnherr
+# uE-Info: 52 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# stretch (and shift!) monotonically increasing coordinate
+
+#======================================================================
+
+sub MW_usage(@)
+{
+ croak("match-warp usage: -w stretch(<field>)\n")
+ unless (@_ == 1);
+ ($MW_cfname) = @_;
+
+ $MW_cfnr = fnr($MW_cfname);
+ $MW_cwfnr = wfnr($MW_cfname);
+
+ croak("$0: field $MW_cfname must be numeric\n")
+ unless numberp($ants_[0][$MW_cfnr]);
+ for (my($r)=1; $r<@ants_; $r++) {
+ croak("$0: field $MW_cfname of rec $r must be numeric ($ants_[$r][$MW_cfnr])\n")
+ unless numberp($ants_[$r][$MW_cfnr]);
+ croak("$0: field $MW_cfname must be monotonically increasing (" .
+ "$ants_[$r][$MW_cfnr]=>$ants_[$r-1][$MW_cfnr]" .
+ ")\n")
+ unless ($ants_[$r][$MW_cfnr] >= $ants_[$r-1][$MW_cfnr]);
+ }
+
+ croak("$0: warp-file field $MW_cfname must be numeric\n")
+ unless numberp($wf_[0][$MW_cwfnr]);
+ for (my($r)=1; $r<@wf_; $r++) {
+ croak("$0: warp-file field $MW_cfname of rec $r must be numeric ($wf_[$r][$MW_cwfnr])\n")
+ unless numberp($wf_[$r][$MW_cwfnr]);
+ croak("$0: warp-file field $MW_cfname must be monotonically increasing (" .
+ "$wf_[$r][$MW_cwfnr]=>$wf_[$r-1][$MW_cwfnr]" .
+ ")\n")
+ unless ($wf_[$r][$MW_cwfnr] >= $wf_[$r-1][$MW_cwfnr]);
+ $MW_unwarped[$r] = $wf_[$r][$MW_cwfnr];
+ }
+}
+
+#======================================================================
+
+sub MW_warp(@)
+{
+ my($offset,$scale) = @_;
+
+ for (my($r)=0; $r<@wf_; $r++) {
+ $wf_[$r][$MW_cwfnr] = $MW_unwarped[$r]*$scale + $offset;
+ }
+}
+
+#======================================================================
+
+1;
--- a/.nminterp.linear
+++ b/.nminterp.linear
@@ -1,9 +1,9 @@
#======================================================================
# . N M I N T E R P . L I N E A R
# doc: Wed Nov 22 21:01:09 2000
-# dlm: Thu Jan 23 13:39:13 2014
+# dlm: Fri Aug 30 13:26:58 2019
# (c) 2000 A.M. Thurnherr
-# uE-Info: 19 61 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 101 9 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# linear interpolation for [fillgaps]
@@ -96,6 +96,9 @@
if ($xf>=0 && equal($ants_[$cr][$xf],$ants_[$nvr[$f]][$xf])) ||
($xf<0 && $cr==$nvr[$f]);
+ if ($xf < 0) {
+ die("$nvr[$f] - $lvr") unless ($nvr[$f] - $lvr > 0);
+ }
my($sc) = ($xf >= 0)
? ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / ($ants_[$nvr[$f]][$xf] - $ants_[$lvr][$xf])
: ($cr - $lvr) / ($nvr[$f] - $lvr);
new file mode 100644
--- /dev/null
+++ b/.nminterp.nan
@@ -0,0 +1,71 @@
+#======================================================================
+# . N M I N T E R P . N A N
+# doc: Sun Jul 6 20:37:31 2014
+# dlm: Sun Jul 6 20:39:26 2014
+# (c) 2014 A.M. Thurnherr
+# uE-Info: 67 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# filling gaps with nan [fillgaps]
+
+# HISTORY:
+# Jul 6, 2014: - created from [.nminterp.linear]
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $ISOpts string of allowed options
+# $ISOptsUsage usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$ISOpts = "";
+$ISOptsUsage = "";
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISUsage() mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub ISUsage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISInit(f,xf) init interpolation of field f
+# f field number
+# xf x field number or -1 for %RECNO
+# <ret val> none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub ISInit($$) {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &interpolate(xf,mg,f,lvr,cr) interpolate field f
+# xf x field or -1 for %RECNO
+# mg max allowed x gap
+# f field number to interpolate
+# lvr last record with valid y val
+# cr current record
+# <ret val> interpolated value
+#
+# NB:
+# - handle f == xf
+# - return undef if interpolation cannot be carried out
+# - x VALUES MAY NOT BE MONOTONIC
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub interpolate($$$$$$)
+{
+ my($xf,$mg,$f,$lvr,$cr) = @_;
+
+ return ($xf == $f) ? $ants_[$cr][$xf] : nan;
+} # static scope
+
+#----------------------------------------------------------------------
+
+1;
new file mode 100644
--- /dev/null
+++ b/.nminterp.spline
@@ -0,0 +1,153 @@
+#======================================================================
+# . N M I N T E R P . S P L I N E
+# doc: Wed Nov 22 21:01:09 2000
+# dlm: Fri Sep 23 16:52:31 2011
+# (c) 2000 A.M. Thurnherr
+# uE-Info: 19 50 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# spline interpolation for [fillgaps]
+
+# HISTORY:
+# Apr 3, 2006: - created from [.interp.spline]
+# Jul 28, 2006: - made to work
+# Jul 30, 2006: - added max-gap support
+# - BUG: extended interpolation into last interval
+# - BUG: assertions made it incompatible with decreasing data
+# Aug 8, 2008: - renamed
+# Sep 23, 2011: - removed $xv from interpolate() args
+# - added error when xf<0 (%RECNO)
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# THE FOLLOWING VARIABLES MUST BE SET GLOBALLY (i.e. during loading)
+#
+# $ISOpts string of allowed options
+# $ISOptsUsage usage information string for options
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+$ISOpts = '';
+$ISOptsUsage = '';
+
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISUsage() mangle parameters (options, really)
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub ISUsage() {}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &ISInit(f,xf) init interpolation of field f
+# f field number
+# xf x field number
+# <ret val> none
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+my(@y2);
+
+sub ISInit($$) {
+ my($f,$xf) = @_;
+ my($k,$i,$p,$qn,$sig,$un,@u);
+
+ croak(".nminterp.spline: cannot use %RECNO as x-field (implementation restriction)\n")
+ unless ($xf >= 0);
+
+ my($pi,$ppi,$li); # find idx of first 3 & last valid recs
+ for (my($r)=0; $r<=$#ants_; $r++) {
+ if (numberp($ants_[$r][$f])) {
+ unless (defined($ppi)) {
+ $ppi = $pi; $pi = $i; $i = $r;
+ }
+ $li = $r;
+ }
+ }
+ unless (defined($ppi)) {
+ &antsInfo("WARNING: not enough <$antsLayout[$f]> data for spline interpolation");
+ return;
+ }
+
+
+ my($fpi) = $pi;
+ $y2[$pi][$f] = $u[$pi] = 0;
+
+ for (; $i<=$li; $i++) {
+ next unless numberp($ants_[$i][$f]);
+
+ $sig = ($ants_[$pi][$xf]-$ants_[$ppi][$xf]) /
+ ($ants_[$i][$xf]-$ants_[$ppi][$xf]);
+ $p = $sig*$y2[$pi][$f] + 2;
+ $y2[$i][$f] = ($sig-1)/$p;
+ $u[$i] = ($ants_[$i][$f]-$ants_[$pi][$f]) /
+ ($ants_[$i][$xf]-$ants_[$pi][$xf])
+ - ($ants_[$pi][$f]-$ants_[$ppi][$f]) /
+ ($ants_[$pi][$xf]-$ants_[$ppi][$xf]);
+ $u[$i] = (6*$u[$i]/($ants_[$i][$xf]-$ants_[$ppi][$xf]) -
+ $sig * $u[$pi]) / $p;
+
+ $ppi = $pi; $pi = $i;
+ }
+
+ $qn = $un = 0;
+
+ $y2[$li][$f] = ($un-$qn*$u[$pi]) / ($qn*$y2[$pi][$f]+1);
+ my($pk) = $li;
+ for ($k=$pi; $k>=$fpi; $k--) {
+ next unless numberp($ants_[$k][$f]);
+ $y2[$k][$f] = $y2[$k][$f] * $y2[$pk][$f] + $u[$k];
+ $pk = $k;
+ }
+}
+
+#+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+#
+# &interpolate(xf,mg,f,lvr,cr) interpolate field f
+# xf x field
+# mg max allowed x gap
+# f field number to interpolate
+# lvr last record with valid y val
+# cr current record
+# <ret val> interpolated value
+#
+# NB:
+# - handle f == xf
+# - return undef if interpolation cannot be carried out
+# - x VALUES MAY NOT BE MONOTONIC
+#
+#++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
+
+sub interpolate($$$$$$)
+{
+ my($xf,$mg,$f,$lvr,$cr) = @_;
+
+ return $ants_[$cr][$xf] if ($xf == $f);
+
+ return $ants_[$lvr][$f] # edge values are ok
+ if equal($ants_[$cr][$xf],$ants_[$lvr][$xf]);
+
+ my($nvr1) = $cr + 1; # find next
+ while ($nvr1 <= $#y2 && !numberp($y2[$nvr1][$f])) { $nvr1++; }
+ return undef # none or gap too large
+ if ($nvr1>$#ants_ || abs($ants_[$nvr1][$xf]-$ants_[$lvr][$xf])>$mg);
+ return $ants_[$nvr1][$f] # edge values are ok
+ if equal($ants_[$cr][$xf],$ants_[$nvr1][$xf]);
+
+ my($nvr2) = $nvr1 + 1;
+ while ($nvr2 <= $#y2 && !numberp($y2[$nvr2][$f])) { $nvr2++; }
+
+ my($h) = $ants_[$nvr1][$xf] - $ants_[$lvr][$xf];
+ my($a) = ($ants_[$nvr1][$xf] - $ants_[$cr][$xf]) / $h;
+ my($b) = ($ants_[$cr][$xf] - $ants_[$lvr][$xf]) / $h;
+
+ return $a*$ants_[$lvr][$f] + $b*$ants_[$nvr1][$f] +
+ (($a*$a*$a-$a) * $y2[$nvr1][$f] +
+ ($b*$b*$b-$b) * $y2[$nvr2][$f]) *
+ ($h*$h)/6;
+}
+
+#----------------------------------------------------------------------
+
+1;
new file mode 100644
--- /dev/null
+++ b/.xycontour.sigma0
@@ -0,0 +1,47 @@
+#======================================================================
+# . I S O P L E T H . S I G M A 0
+# doc: Tue Dec 13 21:50:18 2005
+# dlm: Wed Dec 14 09:40:50 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 44 43 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 14, 2005: - created
+
+# NOTES:
+# - assumes T/S fields to be called 'temp' and 'salin'
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+unless (defined($P{ITS})) {
+ &antsInfo("using default %ITS=90");
+ &antsAddParams(ITS,90);
+}
+
+sub yfield_isopleth($) # return y-field based on x-field
+{
+ my($xf) = @_;
+ return $xf eq 'salin' ? 'temp' : 'salin';
+}
+
+sub bracket_isopleth($$$$) # heuristically bracket y field
+{
+ my($xf,$xv,$iv,$prev_y) = @_;
+
+ if ($xf eq 'salin') {
+ return (-10,30);
+ } else {
+ return (34,37);
+ }
+}
+
+sub eval_isopleth($$$) # evaluate function for isopleth
+{
+ my($xf,$x,$y) = @_;
+
+ return $xf eq 'salin' ? sigma($x,$y,0,0)
+ : sigma($y,$x,0,0)
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/.xycontour.sigma2
@@ -0,0 +1,47 @@
+#======================================================================
+# . I S O P L E T H . S I G M A 2
+# doc: Tue Dec 13 21:50:18 2005
+# dlm: Wed Dec 14 09:11:01 2005
+# (c) 2005 A.M. Thurnherr
+# uE-Info: 19 26 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 13, 2005: - created
+
+# NOTES:
+# - assumes T/S fields to be called 'temp' and 'salin'
+
+require "$ANTS/libEOS83.pl"; # load equation of state
+
+unless (defined($P{ITS})) {
+ &antsInfo("using default %ITS=90");
+ &antsAddParams(ITS,90);
+}
+
+sub yfield_isopleth($) # return y-field based on x-field
+{
+ my($xf) = @_;
+ return $xf eq 'salin' ? 'temp' : 'salin';
+}
+
+sub bracket_isopleth($$$$) # heuristically bracket y field
+{
+ my($xf,$xv,$iv,$prev_y) = @_;
+
+ if ($xf eq 'salin') {
+ return (-10,30);
+ } else {
+ return (34,37);
+ }
+}
+
+sub eval_isopleth($$$) # evaluate function for isopleth
+{
+ my($xf,$x,$y) = @_;
+
+ return $xf eq 'salin' ? sigma($x,$y,0,2000)
+ : sigma($y,$x,0,2000)
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/OLD
@@ -0,0 +1,721 @@
+#======================================================================
+# O L D
+# doc: Tue Nov 26 21:59:40 2013
+# dlm: Mon Jul 1 15:15:19 2019
+# (c) 2017 A.M. Thurnherr
+# uE-Info: 698 42 NIL 0 0 70 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Nov 26, 2013: - created
+# Dec 1, 2013: - removed HDG_ROT
+# - added support for IMP data gaps
+# Mar 4, 2014: - added support for DATA_SOURCE_ID
+# May 22, 2014: - use +/-2deg to assess quality of heading offset
+# May 23, 2014: - added $dhist_binsize, $dhist_min_pirom
+# Jul 27, 2016: - updated variable names for consistency
+# Jul 28, 2016: - major re-write if merging routines
+# Jul 29, 2016: - cosmetics
+# - increased heading-offset resolution from 2 to 1 degrees
+# - BUG: inconsistent heading definition used (from old IMP with
+# confused coordinates)
+# - modified initial timelag guess (there was a bug and it is
+# likely more robust based on end time rather than start time)
+# Aug 5, 2016: - BUG: weird statement accessing LADCP_begin-1
+# - BUG: DSID of first ensemble was not left original
+# Aug 22, 2016: - major changes to timelagging
+# Aug 23, 2016: - changed semantics for removing ensembles with bad attitudes:
+# instead of setting attitude to undef (or large pitch/roll),
+# clearEns() is used
+# Aug 24, 2016: - overhauled time-lagging
+# Aug 25, 2016: - significant code cleanup
+# Aug 26, 2016: - added _hdg_err.ps output plot
+# Oct 13, 2016: - made hdg nan for invalid records (BUG with current versions of IMP+LADCP, IMPatchPD0)
+# Nov 22, 2016: - added heading-offset plot
+# - added sensor info to plots
+# Nov 29, 2016: - added stats to compass error plot
+# Dec 29, 2016: - improved histogram plot
+# Nov 16, 2017: - adapted rot_vecs() to KVM coordinates
+# - made sensor information optional in
+# Nov 20, 2017: - major code cleanup
+# Nov 22, 2017: - replaced "IMP" output in routines used by KVH by "IMU"
+# Dec 8, 2017: - replaced remaing "IMP" output (e.g. in plot) by "IMU"
+# May 22, 2018: - added data trimming to rot_vec
+# May 23, 2018: - added horizontal field strength to mag calib plot
+# - added support for S/N 8 board inside UL WH300 (neg_piro == 2)
+# May 24, 2018: - continued working (coord_trans == 2)
+# May 25, 2018: - continued working
+# May 30, 2018: - BUG: trimming did not re-calculate elapsed time
+# Jun 5, 2018: - BUG: on -c main magnetometer and accelerometer vectors were
+# modified twice
+# Jun 7, 2018: - relaxed definition of -e
+# Jun 9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does not matter in practice)
+
+#----------------------------------------------------------------------
+# gRef() library
+#----------------------------------------------------------------------
+
+sub pl_mag_calib_begin($$$) # initialize mag_calib plot
+{
+ my($pfn,$plotsize,$axlim) = @_;
+ my($xmin,$xmax) = (-$axlim,$axlim);
+ my($ymin,$ymax) = (-$axlim,$axlim);
+ GMT_begin($pfn,"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+ GMT_psxy('-Sc0.05');
+}
+
+sub pl_mag_calib_plot($$$) # plot data point
+{
+ my($valid,$magX,$magY) = @_;
+ if ($valid) { print(GMT "> -Wred -Gred\n$magX $magY\n"); }
+ else { print(GMT "> -Wgreen -Ggreen\n$magX $magY\n"); }
+}
+
+sub pl_mag_calib_end($$) # finish mag_calib plot
+{
+ my($axlim,$HF_mag,$sensor_info) = @_;
+
+ GMT_psxy('-Sc0.1 -Gblue'); # calibration circle
+ for (my($a)=0; $a<2*$pi; $a+=0.075) {
+ printf(GMT "%f %f\n",$HF_mag*sin($a),$HF_mag*cos($a));
+ }
+
+ if ($axlim <= 0.1) { # axes labels
+ GMT_psbasemap('-Bg1a.04f.001:"X Magnetic Field [Gauss]":/g1a0.02f0.001:"Y Magnetic Field [Gauss]":WeSn');
+ } else {
+ GMT_psbasemap('-Bg1a.1f.01:"X Magnetic Field [Gauss]":/g1a0.1f0.01:"Y Magnetic Field [Gauss]":WeSn');
+ }
+ GMT_unitcoords(); # horizontal field strength
+ GMT_pstext('-F+f12,Helvetica,blue+jTR -N');
+ printf(GMT "0.98 0.98 HF = %.2f Gauss\n",$HF_mag);
+
+ printf(GMT "0.98 0.94 $sensor_info\n") # sensor info
+ if ($sensor_info ne '');
+
+ GMT_pstext('-F+f14,Helvetica,blue+jTL -N'); # profile id
+ printf(GMT "0.01 1.06 $P{profile_id}\n");
+ GMT_end();
+}
+
+sub rot_vecs(@) # rotate & output IMU vector data
+{
+ my($coord_trans,$min_elapsed,$max_elapsed,$plot_milapsed,$plot_malapsed) = @_; # negate KVH pitch/roll data if first arg set to 1
+ $min_elapsed = 0 unless defined($min_elapsed);
+ $max_elapsed = 9e99 unless defined($max_elapsed);
+ $plot_milapsed = $min_elapsed unless defined($plot_milapsed);
+ $plot_malapsed = $max_elapsed unless defined($plot_malapsed);
+
+ while (&antsIn()) {
+ next if ($ants_[0][$elapsedF] < $min_elapsed); # trim data
+ last if ($ants_[0][$elapsedF] > $max_elapsed);
+
+ my($cpiro) = -1; # current pitch/roll accelerometer
+ my(@R); # rotation matrix
+ for (my($i)=0; $i<@vecs; $i++) { # rotate vector data
+ if ($piro[$i][0] != $cpiro) { # next sensor chip
+ $cpiro = $piro[$i][0];
+
+ my($accX) = $ants_[0][$vecs[$cpiro][0]];
+ my($accY) = $ants_[0][$vecs[$cpiro][1]];
+ my($accZ) = $ants_[0][$vecs[$cpiro][2]];
+
+ if ($coord_trans == 2) { # S/N 8 inside UL WH300
+ $accY *= -1; $accZ *= -1;
+ }
+
+ my($roll) = atan2($accY,$accZ); # eqn 25 from Freescale AN3461
+ my($pitch) = atan2($accX,sqrt($accY**2+$accZ**2)); # eqn 26 from <Freescale AN3461
+ if ($coord_trans == 1) { # KVH
+ $pitch *= -1;
+ $roll *= -1;
+ }
+ $ants_[0][$piro[$i][1]] = deg($pitch); # add pitch/roll to data
+ $ants_[0][$piro[$i][2]] = deg($roll);
+
+ my($sp) = sin($pitch); my($cp) = cos($pitch); # define rotation matrix
+ my($sr) = sin($roll); my($cr) = cos($roll);
+ @R = ([ $cp, 0, -$sp ],
+ [-$sp*$sr, $cr, -$cp*$sr],
+ [ $sp*$cr, $sr, $cp*$cr]);
+ }
+ my($xval) = $ants_[0][$vecs[$i][0]];
+ my($yval) = $ants_[0][$vecs[$i][1]];
+ my($zval) = $ants_[0][$vecs[$i][2]];
+
+ if ($coord_trans == 2) { # S/N 8 inside UL WH300
+ $yval *= -1; $zval *= -1;
+ }
+
+ $ants_[0][$vecs[$i][0]] = ($xval-$bias[$i][0]) * $R[0][0] +
+ ($yval-$bias[$i][1]) * $R[0][1] +
+ ($zval-$bias[$i][2]) * $R[0][2];
+ $ants_[0][$vecs[$i][1]] = ($xval-$bias[$i][0]) * $R[1][0] +
+ ($yval-$bias[$i][1]) * $R[1][1] +
+ ($zval-$bias[$i][2]) * $R[1][2];
+ $ants_[0][$vecs[$i][2]] = ($xval-$bias[$i][0]) * $R[2][0] +
+ ($yval-$bias[$i][1]) * $R[2][1] +
+ ($zval-$bias[$i][2]) * $R[2][2];
+ }
+
+ my($magX) = $ants_[0][$magXF];
+ my($magY) = $ants_[0][$magYF];
+ my($magZ) = $ants_[0][$magZF];
+ my($accX) = $ants_[0][$accXF];
+ my($accY) = $ants_[0][$accYF];
+ my($accZ) = $ants_[0][$accZF];
+
+ my($HF) = sqrt($magX**2+$magY**2);
+ my($valid)= ($HF >= $minfac*$HF_mag) && ($HF <= $maxfac*$HF_mag);
+ my($hdg) = $valid ? mag_heading($magX,$magY) : nan;
+
+ &antsOut($ants_[0][$elapsedF]-$min_elapsed,$ants_[0][$tempF],
+ RDI_pitch($ants_[0][$pitchF],$ants_[0][$rollF]),$ants_[0][$rollF],
+ $hdg,$accX,$accY,$accZ,$magX,$magY,$magZ,
+ vel_u($HF,$hdg),vel_v($HF,$hdg),$valid);
+
+ pl_mag_calib_plot($valid,$magX,$magY)
+ if defined($P{profile_id}) &&
+ ($ants_[0][$elapsedF] >= $plot_milapsed) &&
+ ($ants_[0][$elapsedF] <= $plot_malapsed);
+ }
+}
+
+#----------------------------------------------------------------------
+# LADCP merging library
+#----------------------------------------------------------------------
+
+#-----------------------------------------------------
+# Instrument Offset Estimation
+#
+# 1: resolution of histogram in deg
+# 1 deg okay for good sensors
+# 2 deg sometimes required (2016 P18 003 2nd sensor)
+# 2: min tilt anom to consider for offset estimation
+# 0.3 deg works even for calm casts
+# increased values improve histogram
+# 2.0 deg is too high for quiet casts (2016 P18 003)
+# 3: minimum fraction of hist mode required
+# 10% default
+# decreasing histogram resolution is better than
+# decreasing this value, I think
+#-----------------------------------------------------
+
+#----------------------------------------------------------------------
+# trim_out_of_water()
+# - attempt to remove out-of-water records from time-series of
+# horizontal acceleration
+#----------------------------------------------------------------------
+
+sub trim_out_of_water($)
+{
+ my($verbose) = @_;
+
+ #--------------------------------------------------------------------------
+ # first-difference horizontal acceleration at full resolution to pass-filter
+ # dangling motion
+ #--------------------------------------------------------------------------
+
+ $IMP{Ah}[0] = sqrt($ants_[0][$accXF]**2+$ants_[0][$accYF]**2);
+ $IMP{dAhdt}[0] = nan;
+ for (my($r)=1; $r<@ants_; $r++) {
+ $IMP{Ah}[$r] = sqrt($ants_[$r][$accXF]**2+$ants_[$r][$accYF]**2);
+ $IMP{dAhdt}[$r] = ($IMP{Ah}[$r]-$IMP{Ah}[$r-1]) / ($ants_[$r][$elapsedF]-$ants_[$r-1][$elapsedF]);
+ }
+
+ #--------------------------------------------------------------------------------------
+ # create 10-s binned time series to calculate rms values of this quantity (dAhdt), and
+ # scale this with cos(sqrt($$pitch**2+$$roll**2)) to dampen underwater peaks (when the
+ # instrument has a large tilt because it is being dragged)
+ # NB: dAhdt, pitch and roll are only set up to last full bin (not so, sum and n)
+ #--------------------------------------------------------------------------------------
+
+ my(@dAhdt,@pitch,@roll,@dAhdt_rms);
+ my(@sum) = my(@sume) = my(@sump) = my(@sumr) = my(@n) = (0);
+
+ my($bin_start) = $ants_[0][$elapsedF];
+ for (my($r)=1; $r<@ants_; $r++) {
+
+ if ($ants_[$r][$elapsedF] - $bin_start <= 10) { # within 10-s bin
+ $sum[$#sum] += $IMP{dAhdt}[$r]; # sums
+ $sume[$#sume] += $ants_[$r][$elapsedF];
+ $sump[$#sump] += $ants_[$r][$pitchF];
+ $sumr[$#sumr] += $ants_[$r][$rollF];
+ $n[$#n]++;
+ next;
+ }
+
+ $dAhdt[$#sum] = $sum[$#sum] / $n[$#n]; # bin done => means
+ $elapsed[$#sum] = $sume[$#sume] / $n[$#n];
+ $pitch[$#sum] = $sump[$#sump] / $n[$#n];
+ $roll[$#sum] = $sumr[$#sumr] / $n[$#n];
+
+ my($sumsq) = 0; # sum of squares for rms(accel)
+ for (my($rr)=$r-$n[$#n]; $rr<$r; $rr++) {
+ $sumsq += ($IMP{dAhdt}[$rr] - $dAhdt[$#sum])**2;
+ }
+ $dAhdt_rms[$#sum] = sqrt($sumsq / $n[$#n]);
+
+ push(@sum,$IMP{dAhdt}[$r]); # begin next bin
+ push(@sume,$ants_[$r][$elapsedF]);
+ push(@sump,$ants_[$r][$pitchF]);
+ push(@sumr,$ants_[$r][$rollF]);
+ push(@n,1);
+ $bin_start = $ants_[$r][$elapsedF];
+ }
+
+ #--------------------------------------------
+ # trim beginning/end when IMP is out of water
+ #--------------------------------------------
+
+ my($i,$si);
+ for ($i=int(@dAhdt_rms/2); $i>0; $i--) {
+ last if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0);
+ }
+ if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0) {
+ for ($si=0; $ants_[$si][$elapsedF]<=$elapsed[$i]; $si++) {}
+ splice(@ants_,0,$si);
+ printf(STDERR "\n\t\t%5d leading out-of-water IMU records removed",$si)
+ if ($si>0 && $verbose);
+ } else {
+ print(STDERR "\n\t\tWARNING: no leading out-of-water IMU records detected/removed") if $verbose;
+ }
+
+ for ($i=int(@dAhdt_rms/2); $i<@dAhdt_rms; $i++) {
+ last if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0);
+ }
+ if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0) {
+ for ($si=$#ants_; $ants_[$si][$elapsedF]>=$elapsed[$i]; $si--) {}
+ my($rem) = @ants_ - $si;
+ splice(@ants_,$si);
+ printf(STDERR "\n\t\t%5d trailing out-of-water IMU records removed",$rem)
+ if ($rem>0 && $verbose);
+ } else {
+ print(STDERR "\n\t\tWARNING: no trailing out-of-water IMU records detected/removed") if $verbose;
+ }
+
+ printf(STDERR "\n\t\tcast duration : %.1f min",
+ ($ants_[$#ants_][$elapsedF] - $ants_[0][$elapsedF]) / 60)
+ if $verbose;
+}
+
+#--------------------------------------------------------------------------------
+# prep_piro_IMP()
+# - calculate pitch/roll offsets & tilt azimuth for IMP
+# - during an attempt to improve time lagging for the 2015 IOPAS data set,
+# it was noticed that one particular instrument, WHM300#12973 maxes out
+# pitch at 27.36 degrees, whereas the roll may not be maxed out at 28.76 deg,
+# the max observed during the cruise.
+# - therefore, IMP{TILT_AZIM} and IMP{TILT_ANOM} are calculated here, first,
+# with a pitch/roll cutoff value of 29 degrees
+# - after the time lagging, when the LADCP start and end times are known,
+# the TILT values are re-calculated without the pitch/roll limit, and
+# using only the correct time range
+#---------------------------------------------------------------------------------
+
+sub prep_piro_IMP($)
+{
+ my($verbose) = @_;
+ my($RDI_max_tilt) = 29;
+ my($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
+
+ for (my($r)=0; $r<@ants_; $r++) {
+ next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
+ $nPR++;
+ $IMP_pitch_mean += min($ants_[$r][$pitchF],$RDI_max_tilt);
+ $IMP_roll_mean += min($ants_[$r][$rollF],$RDI_max_tilt);
+ }
+ $IMP_pitch_mean /= $nPR;
+ $IMP_roll_mean /= $nPR;
+ printf(STDERR "\n\t\tIMU mean pitch/roll : %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
+ if $verbose;
+
+ for (my($r)=0; $r<@ants_; $r++) {
+ next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
+ $IMP{TILT_AZIMUTH}[$r] = tilt_azimuth(min($ants_[$r][$pitchF],$RDI_max_tilt)-$IMP_pitch_mean,
+ min($ants_[$r][$rollF],$RDI_max_tilt) -$IMP_roll_mean);
+ $IMP{TILT_ANOM}[$r] = angle_from_vertical(min($ants_[$r][$pitchF],$RDI_max_tilt)-$IMP_pitch_mean,
+ min($ants_[$r][$rollF],$RDI_max_tilt) -$IMP_roll_mean);
+ }
+ return ($IMP_pitch_mean,$IMP_roll_mean);
+}
+
+#----------------------------------------------------------------------
+# prep_piro_LADCP()
+# - calculate pitch/roll offsets & tilt azimuth of LADCP
+#----------------------------------------------------------------------
+
+sub prep_piro_LADCP($)
+{
+ my($verbose) = @_;
+
+ my($LADCP_pitch_mean,$LADCP_roll_mean) = (0,0);
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+ gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+ $LADCP_pitch_mean += $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH};
+ $LADCP_roll_mean += $LADCP{ENSEMBLE}[$ens]->{ROLL};
+ }
+ $LADCP_pitch_mean /= ($LADCP_end-$LADCP_begin+1);
+ $LADCP_roll_mean /= ($LADCP_end-$LADCP_begin+1);
+ printf(STDERR "\n\t\tLADCP mean pitch/roll : %.1f/%.1f deg",$LADCP_pitch_mean,$LADCP_roll_mean)
+ if $verbose;
+
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH} =
+ tilt_azimuth($LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$LADCP_pitch_mean,
+ $LADCP{ENSEMBLE}[$ens]->{ROLL}-$LADCP_roll_mean);
+ $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} =
+ angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$LADCP_pitch_mean,
+ $LADCP{ENSEMBLE}[$ens]->{ROLL}-$LADCP_roll_mean);
+ }
+
+ print(STDERR "\n") if $verbose;
+ return ($LADCP_pitch_mean,$LADCP_roll_mean);
+}
+
+#------------------------------------------------------------------
+# sub calc_hdg_offset()
+# - estimate heading offset from tilt time series
+# - returns heading offset and updated IMP mean tilts
+# - also creates diagnostic plot with pl_hdg_offset()
+#------------------------------------------------------------------
+
+sub pl_hdg_offset($@)
+{
+ my($dhist_binsize,$modefrac,@dhist) = @_;
+
+ my($plotsize) = '13c';
+ my($xmin,$xmax) = (-180.5,180.5);
+ my($ymin) = 0;
+ my($ymax) = 1.05 * $dhist[$HDG_offset];
+
+ GMT_begin("$P{profile_id}${opt_a}_hdg_offset.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+ GMT_psxy("-Sb${dhist_binsize}u -GCornFlowerBlue");
+ for (my($i)=0; $i<@dhist; $i++) {
+ next unless $dhist[$i];
+ printf(GMT "%f $dhist[$i]\n",$i*$dhist_binsize>180 ? $i*$dhist_binsize-360 : $i*$dhist_binsize);
+ }
+ GMT_psbasemap('-Bg45a90f15:"IMU Heading Offset [\260]":/ga100f10:"Frequency":WeSn');
+ GMT_unitcoords();
+ GMT_pstext('-F+f14,Helvetica,CornFlowerBlue+jTR -N');
+ printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+ GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+ printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
+ GMT_end();
+}
+
+sub calc_hdg_offset($)
+{
+ my($verbose) = @_;
+
+ print(STDERR "\n\tRe-calculating IMU pitch/roll anomalies") if $verbose;
+ ($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+ if ($r < 0 && $ens == $LADCP_begin) {
+ $r = int(($LADCP{ENSEMBLE}[++$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT})
+ while ($r < 0);
+ printf(STDERR "\n\tIMU data begin with instrument already in water => skipping %ds of LADCP data",
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME})
+ if ($verbose);
+ $LADCP_begin = $ens;
+ }
+ if ($r > $#ants_) {
+ printf(STDERR "\n\tIMU data end while instrument is still in water => truncating %ds of LADCP data",
+ $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME})
+ if ($verbose);
+ $LADCP_end = $ens - 1;
+ last;
+ }
+ next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+ $nPR++;
+ $IMP_pitch_mean += $ants_[$r][$pitchF];
+ $IMP_roll_mean += $ants_[$r][$rollF];
+ }
+ $IMP_pitch_mean /= $nPR;
+ $IMP_roll_mean /= $nPR;
+
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+ next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+ $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} =
+ $IMP{TILT_AZIMUTH}[$r] = tilt_azimuth($ants_[$r][$pitchF]-$IMP_pitch_mean,
+ $ants_[$r][$rollF] -$IMP_roll_mean);
+ $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} =
+ $IMP{TILT_ANOM}[$r] = angle_from_vertical($ants_[$r][$pitchF]-$IMP_pitch_mean,
+ $ants_[$r][$rollF] -$IMP_roll_mean);
+ }
+
+ printf(STDERR "\n\t\tIMU mean pitch/roll: %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
+ if $verbose;
+
+ if (defined($opt_s)) {
+ $HDG_offset = $opt_s;
+ printf(STDERR "\n\tHEADING_OFFSET = %.1f (set by -s)\n",$HDG_offset) if $verbose;
+ } else {
+ my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
+ croak("$0: cannot decode -e $opt_e\n")
+ unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
+
+ my(@dhist); my($nhist) = my($modeFreq) = 0;
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+ next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+ next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
+ abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
+ $dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
+ $nhist++;
+ }
+ croak("$0: empty histogram\n")
+ unless ($nhist);
+
+ $HDG_offset = 0;
+ for (my($i)=1; $i<@dhist-1; $i++) { # make sure mode is not on edge
+ $HDG_offset = $i if ($dhist[$i] >= $dhist[$HDG_offset]);
+ }
+ my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
+
+ $HDG_offset *= $dhist_binsize;
+ pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
+
+ if ($opt_f) {
+ printf(STDERR "\n\nIGNORED WARNING (-f): Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100)
+ if ($modefrac < $dhist_min_mfrac);
+ } else {
+ croak(sprintf("\n$0: Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100))
+ if ($modefrac < $dhist_min_mfrac);
+ }
+
+ printf(STDERR "\n\t") if $verbose;
+ printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac) if $verbose;
+ }
+ return ($HDG_offset,$IMP_pitch_mean,$IMP_roll_mean);
+}
+
+#-----------------------------------------------------------
+# rot_IMP()
+# - rotate IMP Data Into LADCP Instrument Coords
+# - also replaced pitch/roll by corresponding anomalies!!!
+#-----------------------------------------------------------
+
+sub rot_IMP($)
+{
+ my($verbose) = @_;
+ my($crho) = cos(rad($HDG_offset));
+ my($srho) = sin(rad($HDG_offset));
+
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+
+ if (numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF])) { # pitch/roll
+ my($rot_p) = (($ants_[$r][$pitchF]-$IMP_pitch_mean)*$crho +
+ ($ants_[$r][$rollF]-$IMP_roll_mean)*$srho);
+ my($rot_r) = (-($ants_[$r][$pitchF]-$IMP_pitch_mean)*$srho +
+ ($ants_[$r][$rollF]-$IMP_roll_mean)*$crho);
+ $ants_[$r][$pitchF] = $rot_p;
+ $ants_[$r][$rollF] = $rot_r;
+ }
+
+ $ants_[$r][$hdgF] = angle_pos($ants_[$r][$hdgF] - $HDG_offset)
+ if numberp($ants_[$r][$hdgF]);
+ }
+
+ my($rot_p) = $IMP_pitch_mean * $crho + $IMP_roll_mean * $srho; # mean pitch roll
+ my($rot_r) = -$IMP_pitch_mean * $srho + $IMP_roll_mean * $crho;
+ $IMP_pitch_mean = $rot_p;
+ $IMP_roll_mean = $rot_r;
+
+ print(STDERR "\n") if $verbose;
+ return ($IMP_pitch_mean,$IMP_roll_mean);
+}
+
+#----------------------------------------------------------------------
+# create_merge_plots()
+# - tilt time series (*_time_lag.ps)
+#----------------------------------------------------------------------
+
+sub create_merge_plots($$$)
+{
+ my($basename,$plotsize,$verbose) = @_;
+
+ #---------------------------------
+ # Tilt Time Series (*_time_lag.ps)
+ #---------------------------------
+
+ my($mint,$maxt) = (99,-99);
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ if (($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 180) ||
+ ($LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 180)) {
+ $mint = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} < $mint);
+ $maxt = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} > $maxt);
+ $mint = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < $mint);
+ $maxt = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} > $maxt);
+ }
+ }
+
+ my($xmin,$xmax) = (-90,90);
+ my($ymin) = round($mint-0.5);
+ my($ymax) = round($maxt+0.5);
+
+ GMT_begin("${basename}_time_lag.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+
+ GMT_psxy('-W2,coral');
+ for (my($ens) = $LADCP_begin + 5;
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 90;
+ $ens++) {
+ printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM});
+ }
+ GMT_psxy('-W1');
+ for (my($ens) = $LADCP_begin + 5;
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 90;
+ $ens++) {
+ printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM});
+ }
+ GMT_psxy('-W2,SeaGreen');
+ for (my($ens) = $LADCP_end - 5;
+ $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 90;
+ $ens--) {
+ printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM});
+ }
+ GMT_psxy('-W1');
+ for (my($ens) = $LADCP_end - 5;
+ $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 90;
+ $ens--) {
+ printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM});
+ }
+
+ GMT_psbasemap('-Bg30a30f10:"Elapsed Time [sec]":/g5a5f1:"Tilt Magnitude [\260]":WeSn');
+ GMT_unitcoords();
+ GMT_pstext('-F+f14,Helvetica,Coral+jTL');
+ printf(GMT "0.52 0.98 downcast\n");
+ GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR');
+ printf(GMT "0.48 0.98 upcast\n");
+ GMT_psxy('-W4,LightSkyBlue');
+ printf(GMT "0.5 0\n0.5 1\n");
+ GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+ printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
+ GMT_end();
+
+ #------------------------------
+ # Heading Errors (*_hdg_err.ps)
+ #------------------------------
+
+ my(@err_binned,@err_nsamp);
+ my($sumErr) = 0; my($nErr) = $LADCP_end - $LADCP_begin + 1;
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ next unless ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < 10);
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+ my($bi) = int($ants_[$r][$hdgF]/5);
+ my($err) = angle_diff($ants_[$r][$hdgF],$LADCP{ENSEMBLE}[$ens]->{HEADING});
+ next unless numberp($err);
+ $err_binned[$bi] += $err; $sumErr += $err;
+ $err_nsamp[$bi]++;
+ }
+ for (my($bi)=0; $bi<@err_nsamp; $bi++) {
+ $err_binned[$bi] = ($err_nsamp[$bi] >= 5)
+ ? $err_binned[$bi]/$err_nsamp[$bi]
+ : undef;
+ }
+ my(@err_dssq);
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ next unless ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < 10);
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+ my($bi) = int($ants_[$r][$hdgF]/5);
+ my($err) = angle_diff($ants_[$r][$hdgF],$LADCP{ENSEMBLE}[$ens]->{HEADING});
+ next unless numberp($err);
+ $err_dssq[$bi] += ($err-$err_binned[$bi])**2;
+ }
+
+ my($xmin,$xmax) = (0,360);
+ my($ymin,$ymax) = (-45,45);
+ GMT_begin("${basename}_hdg_err.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+ GMT_psxy('-Ey/3,CornFlowerBlue');
+ my($sumSq,$sumBe) = my($nSq,$nBe) = (0,0);
+ for (my($bi)=0; $bi<@err_binned; $bi++) {
+ next unless ($err_nsamp[$bi] >= 2);
+ next unless numberp($err_binned[$bi]);
+ $sumSq += $err_binned[$bi]**2; $nSq++;
+ $sumBe += $err_binned[$bi]; $nBe++;
+# printf(GMT "%f %f\n",2.5+5*$bi,$err_binned[$bi]);
+ printf(GMT "%f %f %f\n",2.5+5*$bi,$err_binned[$bi],sqrt($err_dssq[$bi]/($err_nsamp[$bi]-1)));
+ }
+ GMT_psbasemap('-Bg90a45f5:"ADCP Heading [\260]":/g15a15f5:"ADCP Compass Error [\260]":WeSn');
+ GMT_unitcoords();
+ GMT_pstext('-F+f12,Helvetica,CornFlowerBlue+jTR -Gwhite -C25%');
+ printf(GMT "0.98 0.98 rms error = %7.1f \260\n",sqrt($sumSq/$nSq));
+ printf(GMT "0.98 0.94 time-averaged error = %7.1f \260\n",$sumErr/$nErr);
+ printf(GMT "0.98 0.90 heading-averaged error = %7.1f \260\n",$sumBe/$nBe);
+ GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
+ printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
+ GMT_end();
+
+ print(STDERR "\n") if $verbose;
+}
+
+#----------------------------------------------------------------------
+# output_merged()
+# - output merged data
+#----------------------------------------------------------------------
+
+sub output_merged($)
+{
+ my($verbose) = @_;
+
+ my($tazimF) = &antsNewField('tilt_azimuth');
+ my($tanomF) = &antsNewField('tilt_magnitude');
+ my($L_tazimF) = &antsNewField('LADCP_tilt_azimuth');
+ my($L_tanomF) = &antsNewField('LADCP_tilt_magnitude');
+ my($L_elapsedF) = &antsNewField('LADCP_elapsed');
+ my($L_ensF) = &antsNewField('LADCP_ens');
+ my($L_depthF) = &antsNewField('LADCP_depth_estimate');
+ my($L_pitchF) = &antsNewField('LADCP_pitch');
+ my($L_rollF) = &antsNewField('LADCP_roll');
+ my($L_hdgF) = &antsNewField('LADCP_hdg');
+ my($dcF) = &antsNewField('downcast');
+
+ for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+# print(STDERR "$r\[$ens\,$LADCP{ENSEMBLE}[$ens]->{NUMBER}] = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});\n");
+ if ($r<0 || $r>$#ants_) { # ensemble beyond limits of IMU data
+ my(@out);
+ if ($ens >= $LADCP_begin && $ens <= $LADCP_end) {
+ $out[$elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG};
+ $out[$L_elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME};
+ }
+ $out[$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+ $out[$dcF] = ($ens <= $LADCP_bottom);
+ &antsOut(@out);
+ } elsif ($ens < $LADCP_begin || $ens > $LADCP_end) { # pre deplyment or post recovery
+ $ants_[$r][$L_elapsedF] = undef;
+ $ants_[$r][$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+ $ants_[$r][$L_pitchF] = $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
+ $ants_[$r][$L_rollF] = $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
+ $ants_[$r][$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
+ $ants_[$r][$dcF] = nan;
+ &antsOut(@{$ants_[$r]});
+ } else {
+ $ants_[$r][$tazimF] = angle($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} + $HDG_offset);
+ $ants_[$r][$tanomF] = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM};
+ $ants_[$r][$L_tazimF] = $LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH};
+ $ants_[$r][$L_tanomF] = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM};
+ $ants_[$r][$L_elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME};
+ $ants_[$r][$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+ $ants_[$r][$L_depthF] = $LADCP{ENSEMBLE}[$ens]->{DEPTH};
+ $ants_[$r][$L_pitchF] = $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
+ $ants_[$r][$L_rollF] = $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
+ $ants_[$r][$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
+ $ants_[$r][$dcF] = ($ens <= $LADCP_bottom);
+ &antsOut(@{$ants_[$r]});
+ }
+ }
+
+ print(STDERR "\n") if $verbose;
+}
+
+#----------------------------------------------------------------------
+
+1; # return true for all the world to see
--- 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: Mon Apr 23 14:20:58 2018
+# dlm: Wed Apr 10 16:57:59 2019
# (c) 1998 A.M. Thurnherr
-# uE-Info: 216 82 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 217 48 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -214,6 +214,7 @@
# Mar 10, 2017: - BUG: antsCheckDeps() used ctime instead of mtime!!!
# Apr 5, 2017: - BUG: stale file mtime dependency info was not printed correctly
# Apr 23, 2018: - BUG: @antsLayout was not kept up-to-date when layout-changes are allowed
+# Apr 10, 2019: - disabled dependency warnings
# GENERAL NOTES:
# - %P was named without an ants-prefix because associative arrays
@@ -292,7 +293,7 @@
# - by default, tests current <> file
#----------------------------------------------------------------------
-{ my($warned);
+{ my($warned) = 1; # disable dependency warning
sub antsCheckDeps()
{
--- 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: Wed May 25 12:16:39 2016
+# dlm: Fri Aug 30 13:39:50 2019
# (c) 1998 A.M. Thurnherr
-# uE-Info: 104 31 NIL 0 0 70 10 2 4 NIL ofnI
+# uE-Info: 108 62 NIL 0 0 70 10 2 4 NIL ofnI
#======================================================================
# Miscellaneous auxillary functions
@@ -102,6 +102,10 @@
# May 18, 2015: - added antsFindParam()
# Jun 21, 2015: - added antsParam(), modified antsRequireParam()
# May 12, 2016: - added &div2() to prevent division by zero errors
+# Apr 5, 2019: - disabled weird line of code in antsFunUsage() (see comment)
+# - improved error messages in antsFunUsage()
+# - BUG: antsFunUsage did not work with -ve argc (variable argument funs)
+# Aug 30, 2019: - BUG: antsLoadModel() did not respect $ANTS
# fnr notes:
# - matches field names starting with the string given, i.e. "sig" is
@@ -446,8 +450,7 @@
require "$pref.$name";
return $name;
} else { # load from ANTSlib
- my($path) = ($0 =~ m{^(.*)/[^/]*$});
- require "$path/$pref.$name";
+ require "$ANTS/$pref.$name";
return $name;
}
}
@@ -467,8 +470,7 @@
require "$pref.$name";
return ($name,split(',',$args));
} else {
- my($path) = ($0 =~ m{^(.*)/[^/]*$});
- require "$path/$pref.$name";
+ require "$ANTS/$pref.$name";
return ($name,split(',',$args));
}
}
@@ -538,7 +540,7 @@
if (ref($params[0]) && @antsLayout>0 && @params<2*$argc+1) { # default params
my(@newparams); # 2nd test is for abc
- my($npi) = $argc+1;
+ my($npi) = abs($argc)+1;
$listAllRecs = 1; # special flag for list(1)
@@ -551,7 +553,7 @@
return(@newparams);
}
- for (my($i)=1; $i<=$argc; $i++) { # fill cache & do tests
+ for (my($i)=1; $i<=abs($argc); $i++) { # fill cache & do tests
if (defined($params[$i])) {
push(@{$params[0]},&fnr($params[$i]));
push(@newparams,&antsVal($params[0]->[$#{$params[0]}]));
@@ -564,24 +566,27 @@
croak("usage: $msg\n") unless ($npi > $#params);
@params = @newparams;
- } elsif (ref($params[0])) {
- splice(@params,0,$argc+1);
+ } elsif (ref($params[0])) { # remove array ref & list of field names
+ splice(@params,0,abs($argc)+1);
}
if ($argc >= 0) { # argument count
- croak("usage: $msg\n") unless (@params == $argc);
+ croak("usage: $msg [params = @params]\n") unless (@params == $argc);
} else {
- croak("usage: $msg\n") unless (@params >= -$argc);
+ croak("usage: $msg [params = @params])\n") unless (@params >= -$argc);
}
for (my($i)=0; $i<length($types); $i++) { # type checking
$_ = substr($types,$i,1);
SWITCH: {
- last unless defined($params[$i]);
- &antsNoCardErr("",$params[$i]),last SWITCH if (/c/);
- &antsNoIntErr("",$params[$i]),last SWITCH if (/i/);
- &antsNoFloatErr("",$params[$i]),last SWITCH if (/f/);
- &antsNoFileErr("",$params[$i]),last SWITCH if (/F/);
+# 4/5/19: The following line of code prevents proper type checking when one of the
+# arguments is undefined. I do not know under what circumstances the code
+# is required. Therfore I disabled it temporarily.
+# last unless defined($params[$i]);
+ &antsNoCardErr(sprintf("argument #%d in $msg (params = @params)",$i+1),$params[$i]),last SWITCH if (/c/);
+ &antsNoIntErr(sprintf("argument #%d in $msg",$i+1),$params[$i]),last SWITCH if (/i/);
+ &antsNoFloatErr(sprintf("argument #%d in $msg (params = @params)",$i+1),$params[$i]),last SWITCH if (/f/);
+ &antsNoFileErr(sprintf("argument #%d in $msg",$i+1),$params[$i]),last SWITCH if (/F/);
if (/\d/) {
croak("$0: $params[$i] is not a string of length $_\n")
unless ($_ == length($params[$i]));
--- a/libGM76.pl
+++ b/libGM76.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B G M 7 6 . P L
# doc: Sun Feb 20 14:43:47 2011
-# dlm: Sat Feb 2 02:33:11 2019
+# dlm: Sat Apr 6 19:33:05 2019
# (c) 2011 A.M. Thurnherr
-# uE-Info: 144 43 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 34 48 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -27,7 +27,11 @@
# - added GM_strain
# - BUG: Sw usage message had wrong parameter order
# - renamed Sw => GM_VKE; Su_z => GM_shear
-
+# Mar 31, 2019: - updated doc for shear and strain spectra
+# Apr 1, 2019: - added GM_vdiv
+# Apr 5, 2019: - BUG: GM_VKE was erroneous (argument shifting was wrong)
+# - adapted to improved antsFunUsage()
+# Apr 6, 2019: - fiddled during debugging of fs_finestructure issue
require "$ANTS/libEOS83.pl";
@@ -100,7 +104,7 @@
sub GM_VKE(@)
{
my($omega,$m,$lat,$b,$jstar,$N) =
- &antsFunUsage(-5,'fff','[frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+ &antsFunUsage(-5,'fffff','GM_VKE([frequency[1/s]] <vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
if (defined($N)) { # Sw(omega,m)
local($f) = abs(&f($lat));
@@ -110,8 +114,8 @@
my($mstar) = &m($jstar,$N,$omega);
return $GM_E0 * $b * 2 * $f**2/$omega**2/B($omega) * $jstar / ($m+$mstar)**2;
} else { # Sw(m), i.e. integrated over all omega; as in Thurnherr et al., GRL 2015
- $N = $lat; # shift arguments to account for missing omega
- $lat = $m;
+ $N = $jstar; $jstar = $b; # shift arguments to account for missing omega
+ $b = $lat; $lat = $m;
local($f) = abs(&f($lat));
$m = $omega;
undef($omega);
@@ -121,15 +125,38 @@
}
#----------------------------------------------------------------------
+# Vertical Divergence (normalized by N)
+# - implemented to investigate shear-to-vertical divergence ratio
+# - S[w_z/N] = S[w_z] / N^2 = S[w] * m^2 / N^2
+# - GM shear-to-vd variance ratio = 3N/2f
+# - GM strain-to-vd variance ratio = N/2f
+#----------------------------------------------------------------------
+
+sub GM_vdiv(@)
+{
+ my($m,$lat,$b,$jstar,$N) =
+ &antsFunUsage(5,'fffff','GM_vdiv(<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
+ return GM_VKE($m,$lat,$b,$jstar,$N) * $m**2 / $N**2;
+}
+
+#----------------------------------------------------------------------
# Shear and Strain m-spectra (i.e. integrated over f)
+# - shear is buoyancy-frequency normalized
# - spectral density
-# - from Thurnherr et al. (GRL 2015)
+# - from Thurnherr et al. (GRL 2015), which is from email info
+# Eric Kunze, 04/22/2015: First, the standard GHP parameterization that I
+# implement uses either the N-normalized shear or strain spectra. The GM
+# versions of these are
+# S[V_z/N](m) = (3*PI/2)*E0*b*jstar*m^2/(m+mstar)^2 and
+# S[Z_z](m) = (PI/2)*E0*b*jstar*m^2/(m+mstar)^2
+# though in practice mstar is not particularly important because the
+# variance of both these quantities is dominated by high m.
#----------------------------------------------------------------------
sub GM_shear(@)
{
my($m,$lat,$b,$jstar,$N) =
- &antsFunUsage(5,'fffff','<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+ &antsFunUsage(5,'fffff','GM_shear(<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
local($f) = abs(&f($lat));
my($mstar) = &m($jstar,$N);
return 3 * $pi/2 * $GM_E0 * $b * $jstar * $m**2 / ($m+$mstar)**2;
@@ -138,7 +165,7 @@
sub GM_strain(@)
{
my($m,$lat,$b,$jstar,$N) =
- &antsFunUsage(5,'fffff','<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>',@_);
+ &antsFunUsage(5,'fffff','GM_strain(<vertical wavenumber[rad/m]> <lat[deg]> <b[m]> <j*> <N[rad/s]>)',@_);
local($f) = abs(&f($lat));
my($mstar) = &m($jstar,$N);
return $pi/2 * $GM_E0 * $b * $jstar * $m**2 / ($m+$mstar)**2;
--- a/libIMP.pl
+++ b/libIMP.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B I M P . P L
# doc: Tue Nov 26 21:59:40 2013
-# dlm: Sat Jun 9 12:07:15 2018
+# dlm: Mon Jul 1 15:47:57 2019
# (c) 2017 A.M. Thurnherr
-# uE-Info: 52 109 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 689 0 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -49,7 +49,16 @@
# Jun 5, 2018: - BUG: on -c main magnetometer and accelerometer vectors were
# modified twice
# Jun 7, 2018: - relaxed definition of -e
-# Jun 9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does not matter in practice)
+# Jun 9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does
+# not matter in practice)
+# Jun 30, 2019: - renamed -s in IMP+LADCP to -o (option not used in KVH+LADCP)
+# - changed semantics so that offset histogram is plotted even on -o
+# - made create_merge_plots return errors
+# Jul 1, 2019: - BUG: pre-deployment output cleared IMU pitch/roll fields (1 record
+# affected)
+# - modified output of pre-/post-cast records to include all LADCP
+# but no IMU data
+# - BUG: GIMBAL_PITCH had not only been defined for in-water records
#----------------------------------------------------------------------
# gRef() library
@@ -349,6 +358,10 @@
my($verbose) = @_;
my($LADCP_pitch_mean,$LADCP_roll_mean) = (0,0);
+ for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+ gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+ }
for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
@@ -398,7 +411,11 @@
GMT_psbasemap('-Bg45a90f15:"IMU Heading Offset [\260]":/ga100f10:"Frequency":WeSn');
GMT_unitcoords();
GMT_pstext('-F+f14,Helvetica,CornFlowerBlue+jTR -N');
- printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+ if (defined($opt_o)) {
+ printf(GMT "0.99 1.06 -o %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+ } else {
+ printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+ }
GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
GMT_end();
@@ -449,35 +466,36 @@
printf(STDERR "\n\t\tIMU mean pitch/roll: %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
if $verbose;
- if (defined($opt_s)) {
- $HDG_offset = $opt_s;
- printf(STDERR "\n\tHEADING_OFFSET = %.1f (set by -s)\n",$HDG_offset) if $verbose;
- } else {
- my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
- croak("$0: cannot decode -e $opt_e\n")
- unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
+ my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
+ croak("$0: cannot decode -e $opt_e\n")
+ unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
- my(@dhist); my($nhist) = my($modeFreq) = 0;
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
- next unless numberp($IMP{TILT_AZIMUTH}[$r]);
- next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
- abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
- $dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
- $nhist++;
- }
- croak("$0: empty histogram\n")
- unless ($nhist);
+ my(@dhist); my($nhist) = my($modeFreq) = 0;
+ for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
+ next unless numberp($IMP{TILT_AZIMUTH}[$r]);
+ next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
+ abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
+ $dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
+ $nhist++;
+ }
+ croak("$0: empty histogram\n")
+ unless ($nhist);
+ if (defined($opt_o)) { # set heading offset, either with -o
+ $HDG_offset = $opt_o;
+ } else { # or from the data
$HDG_offset = 0;
for (my($i)=1; $i<@dhist-1; $i++) { # make sure mode is not on edge
$HDG_offset = $i if ($dhist[$i] >= $dhist[$HDG_offset]);
}
- my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
-
$HDG_offset *= $dhist_binsize;
- pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
-
+ }
+
+ my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
+ pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
+
+ unless (defined($opt_o)) { # make sure data are consistent (unless -o is used)
if ($opt_f) {
printf(STDERR "\n\nIGNORED WARNING (-f): Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100)
if ($modefrac < $dhist_min_mfrac);
@@ -485,10 +503,17 @@
croak(sprintf("\n$0: Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100))
if ($modefrac < $dhist_min_mfrac);
}
+ }
- printf(STDERR "\n\t") if $verbose;
- printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac) if $verbose;
+ printf(STDERR "\n\t") if $verbose;
+ if ($opt_o) {
+ printf(STDERR "IMU heading offset = -o %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac)
+ if $verbose;
+ } else {
+ printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac)
+ if $verbose;
}
+
return ($HDG_offset,$IMP_pitch_mean,$IMP_roll_mean);
}
@@ -532,6 +557,7 @@
#----------------------------------------------------------------------
# create_merge_plots()
# - tilt time series (*_time_lag.ps)
+# - heading errors (*_hdg_err.ps)
#----------------------------------------------------------------------
sub create_merge_plots($$$)
@@ -654,11 +680,13 @@
GMT_end();
print(STDERR "\n") if $verbose;
+ return (sqrt($sumSq/$nSq),$sumErr/$nErr,$sumBe/$nBe);
}
#----------------------------------------------------------------------
# output_merged()
# - output merged data
+# - there must be exactly one record for each PD0 ensemble
#----------------------------------------------------------------------
sub output_merged($)
@@ -690,16 +718,13 @@
$out[$dcF] = ($ens <= $LADCP_bottom);
&antsOut(@out);
} elsif ($ens < $LADCP_begin || $ens > $LADCP_end) { # pre deplyment or post recovery
- $ants_[$r][$L_elapsedF] = undef;
- $ants_[$r][$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
- $ants_[$r][$L_pitchF] = undef;
- $ants_[$r][$L_rollF] = undef;
- $ants_[$r][$L_hdgF] = undef;
- $ants_[$r][$pitchF] = undef;
- $ants_[$r][$rollF] = undef;
- $ants_[$r][$hdgF] = undef;
- $ants_[$r][$dcF] = ($ens <= $LADCP_bottom);
- &antsOut(@{$ants_[$r]});
+ my(@out); # correct IMP record NOT known
+ $out[$L_elapsedF] = undef;
+ $out[$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
+ $out[$L_pitchF] = $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
+ $out[$L_rollF] = $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
+ $out[$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
+ &antsOut(@out);
} else {
$ants_[$r][$tazimF] = angle($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} + $HDG_offset);
$ants_[$r][$tanomF] = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM};
@@ -710,8 +735,8 @@
$ants_[$r][$L_depthF] = $LADCP{ENSEMBLE}[$ens]->{DEPTH};
$ants_[$r][$L_pitchF] = $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
$ants_[$r][$L_rollF] = $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
- $ants_[$r][$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
- $ants_[$r][$dcF] = ($ens <= $LADCP_bottom);
+ $ants_[$r][$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
+ $ants_[$r][$dcF] = ($ens <= $LADCP_bottom) ? 1 : 0;
&antsOut(@{$ants_[$r]});
}
}
--- 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 Apr 25 17:41:36 2018
+# dlm: Sat Apr 6 19:42:07 2019
# (c) 2011 A.M. Thurnherr
-# uE-Info: 46 27 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 211 34 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -24,6 +24,8 @@
# - added parameter checks to processing-specific corrections
# May 18, 2015: - added pulse length to T_w() and T_w_z()
# Apr 25, 2018: - added eps_VKE() parameterization
+# Apr 5, 2018: - adapted to improved antsFunUsage()
+# - BUG eps_VKE() had erroneous string
require "$ANTS/libvec.pl";
require "$ANTS/libfuns.pl";
@@ -43,7 +45,7 @@
sub eps_VKE(@)
{
my($p0,$c) =
- &antsFunUsage(-1,'.','<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]]',@_);
+ &antsFunUsage(-1,'.','eps_VKE(<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]])',@_);
$c = 0.021 unless defined($c);
return numberp($p0) ? ($p0 / $c) ** 2 : nan;
}
@@ -69,7 +71,7 @@
sub T_ravg(@)
{
my($k,$blen,$plen) =
- &antsFunUsage(-2,'ff','<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]]',@_);
+ &antsFunUsage(-2,'ff','T_ravg(<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;
}
@@ -82,7 +84,7 @@
sub T_fdiff($$)
{
my($k,$dz) =
- &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_);
+ &antsFunUsage(2,'ff','T_fdiff(<vertical wavenumber[rad/s]> <differencing interval[m]>)',@_);
return 1 / sinc($k*$dz/2/$PI)**2;
}
@@ -95,7 +97,7 @@
sub T_interp($$$)
{
my($k,$blen,$dz) =
- &antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_);
+ &antsFunUsage(3,'fff','T_interp(<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>)',@_);
return 1 / sinc($k*$blen/2/$PI)**4 / sinc($k*$dz/2/$PI)**2;
}
@@ -108,7 +110,7 @@
sub T_binavg($$)
{
my($k,$dz) =
- &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_);
+ &antsFunUsage(2,'ff','T_binavg(<vertical wavenumber[rad/s]> <grid resolution[m]>)',@_);
return 1 / sinc($k*$dz/2/$PI)**2;
}
@@ -125,7 +127,7 @@
sub T_tilt($$)
{
my($k,$dprime) =
- &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_);
+ &antsFunUsage(2,'ff','T_tilt(<vertical wavenumber[rad/s]> <d-prime[m]>)',@_);
return $dprime ? 1 / sinc($k*$dprime/2/$PI)**2 : 1;
}
@@ -147,7 +149,7 @@
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]>',@_);
+ &antsFunUsage(-3,'fff','T_UH(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>)',@_);
croak("T_UH($k,$blen,$dz,$range_max): bad parameters\n")
unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);
return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max));
@@ -162,7 +164,7 @@
sub T_ASM($$$$)
{
my($k,$blen,$dz,$range_max) =
- &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+ &antsFunUsage(-3,'fff','T_ASM(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>)',@_);
croak("T_ASM($k,$blen,$dz,$range_max): bad parameters\n")
unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);
return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
@@ -180,14 +182,14 @@
sub T_VI($$$$$)
{
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]>',@_);
+ &antsFunUsage(-4,'ff.f','T_VI(<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($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]>',@_);
+ &antsFunUsage(-4,'ff.ff','T_VI_alt(<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <tilt d-prime[m]>)',@_);
croak("T_VI_alt($k,$blen,$sel,$dz,$range_max): bad parameters\n")
unless ($k>=0 && $blen>0 && $sel ne '' && $dz>0 && $range_max>=0);
croak("$0: tilt-dprime outside range [0..$blen]\n")
@@ -206,9 +208,9 @@
sub T_w(@)
{
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,@_);
+ &antsFunUsage(-4,'ffff',
+ 'T_w([vertical wavenumber[rad/s]] <ADCP bin length[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>)',
+ \@fc,'k',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));
@@ -227,9 +229,9 @@
sub T_w_z(@)
{
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,@_);
+ &antsFunUsage(-4,'ffff',
+ 'T_w_z([vertical wavenumber[rad/s]] <ADCP bin size[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>)',
+ \@fc,'k',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);
new file mode 100644
--- /dev/null
+++ b/libLADCP.pl.orig
@@ -0,0 +1,240 @@
+#======================================================================
+# L I B L A D C P . P L
+# doc: Wed Jun 1 20:38:19 2011
+# dlm: Mon Apr 1 16:13:46 2019
+# (c) 2011 A.M. Thurnherr
+# uE-Info: 27 65 NIL 0 0 70 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Jun 1, 2011: - created
+# Jul 29, 2011: - improved
+# Aug 10, 2011: - made correct
+# Aug 11, 2011: - added "convenient combinations"
+# Aug 18, 2011: - made buoyancy frequency non-constant in S()
+# Jan 4, 2012: - improved T_VI to allow correcting even without superensembles
+# 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()
+# 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()
+# Apr 25, 2018: - added eps_VKE() parameterization
+# Apr 1, 2019: - removed plen parameter from T_w() and T_w_z()
+
+require "$ANTS/libvec.pl";
+require "$ANTS/libfuns.pl";
+
+#------------------------------------------------------------------------------
+# VKE parameterization for epsilon
+#
+# NOTES:
+# - see Thurnherr et al. (GRL 2015)
+# - calculate eps from p0
+# - optional second argument allows free choice of parameterization constant
+# - default value is from paper, which is slightly lower than the current value
+# used in [LADCP_VKE], which applies the parameterization only to spectra
+# passing a few tests
+#------------------------------------------------------------------------------
+
+sub eps_VKE(@)
+{
+ my($p0,$c) =
+ &antsFunUsage(-1,'.','<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]]',@_);
+ $c = 0.021 unless defined($c);
+ return numberp($p0) ? ($p0 / $c) ** 2 : nan;
+}
+
+#------------------------------------------------------------------------------
+# 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
+#------------------------------------------------------------------------------
+#----------------------------------------------------------------------
+# 1. Corrections for individual data acquisition and processing steps
+#----------------------------------------------------------------------
+
+#------------------------------------------------------------------------------
+# T_ravg(k,blen[,plen])
+# - correct for range averaging due to finite pulse and finite receive window
+# - when called with 2 arguments, bin-length == pulse-length assumed
+#------------------------------------------------------------------------------
+
+sub T_ravg(@)
+{
+ 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;
+}
+
+#-------------------------------------------------------------
+# T_fdiff(k,dz)
+# - correct for first differencing on a grid with dz spacing
+#-------------------------------------------------------------
+
+sub T_fdiff($$)
+{
+ my($k,$dz) =
+ &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_);
+ 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($k,$blen,$dz) =
+ &antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_);
+ 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
+#-------------------------------------------------------------------------
+
+sub T_binavg($$)
+{
+ my($k,$dz) =
+ &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_);
+ 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($k,$dprime) =
+ &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_);
+ 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]>',@_);
+ croak("T_UH($k,$blen,$dz,$range_max): bad parameters\n")
+ unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);
+ return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max));
+}
+
+#----------------------------------------------------------------------
+# T_ASM(k,blen,dz,range_max)
+# - re-implemented shear method with simple depth binning
+# - range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+sub T_ASM($$$$)
+{
+ my($k,$blen,$dz,$range_max) =
+ &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_);
+ croak("T_ASM($k,$blen,$dz,$range_max): bad parameters\n")
+ unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0);
+ return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max));
+}
+
+#------------------------------------------------------------
+# 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
+# - sel == nan disables pre-averaging correction
+#------------------------------------------------------------
+
+sub T_VI($$$$$)
+{
+ 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($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("T_VI_alt($k,$blen,$sel,$dz,$range_max): bad parameters\n")
+ unless ($k>=0 && $blen>0 && $sel ne '' && $dz>0 && $range_max>=0);
+ croak("$0: tilt-dprime outside range [0..$blen]\n")
+ unless ($dprime>=0 && $dprime<=$blen);
+ 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 length[m]> <output grid resolution[m]> <range max[m]>',
+ \@fc,'k',undef,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));
+ }
+}
+
+#----------------------------------------------------------------------
+# T_w_z(k,blen,dz,range_max)
+# - vertical-velocity method of Thurnherr (IEEE 2011)
+# - first differencing of gridded shear to calculate dw/dz
+# - NB: grid-scale differentiation assumed
+# - range_max == 0 disables tilt correction
+#----------------------------------------------------------------------
+
+{ 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,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);
+ }
+}
+
+1;
--- a/libstats.pl
+++ b/libstats.pl
@@ -1,9 +1,9 @@
#======================================================================
-# L I B S T A T S . P L
+# . . / L I B / L I B S T A T S . P L
# doc: Wed Mar 24 13:59:27 1999
-# dlm: Sat Jan 31 15:51:14 2015
+# dlm: Thu Mar 12 15:11:24 2020
# (c) 1999 A.M. Thurnherr
-# uE-Info: 150 0 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 45 44 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -36,10 +36,37 @@
# - added fdiff()
# Jan 30, 2015: - added log_avg()
# - added noise_avg()
+# Mar 26, 2019: - added regress()
require "$ANTS/libfuns.pl";
#----------------------------------------------------------------------
+# simple linear regression
+# - returns regression coefficient (slope)
+#----------------------------------------------------------------------
+
+sub regress(@)
+{
+ my($n) = @_/2;
+ return nan unless ($n > 1);
+
+ my(@X) = @_[0..$n-1];
+ my(@Y) = @_[$n..2*$n-1];
+
+ my($sumx) = 0;
+ for (my($i)=0; $i<$n; $i++) { $sumx += $X[$i]; }
+ my($midx) = $sumx / @Y;
+ my($rcoeff) = my($sumxdsq) = 0;
+ for (my($i)=0; $i<$n; $i++) {
+ my($xd) = $X[$i] - $midx;
+ $sumxdsq += $xd**2;
+ $rcoeff += $xd * $Y[$i];
+ }
+ $rcoeff /= $sumxdsq;
+ return $rcoeff;
+}
+
+#----------------------------------------------------------------------
# estimate stderr given stddev & degrees of freedom
# - return nan for dof <= 0
#----------------------------------------------------------------------