--- a/LADCPproc.bestLag
+++ b/LADCPproc.bestLag
@@ -1,14 +1,22 @@
#======================================================================
# L A D C P P R O C . B E S T L A G
# doc: Tue Sep 28 21:58:48 2010
-# dlm: Tue Oct 26 13:37:19 2010
+# dlm: Wed Jan 5 19:44:09 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 58 94 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 67 28 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
# - first lag is always(?) nan unless CTD is turned on with LADCP in water
+# HISTORY:
+# Sep 28, 2010: - created
+# Dec 9, 2010: - adapted to %CTD
+# Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags
+# Jan 5, 2011: - changed first guess from 80% down to 10% down
+# - added LADCP time lag to %PARAMs
+# - added support of -i
+
sub interp_LADCP_w($$)
{
my($elapsed,$ens) = @_;
@@ -34,8 +42,8 @@
for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
my($Li) = $Ci + $Llag;
next if ($Li<0 || $Li>=$opt_w);
- next unless numberp($CTD_w[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
- $sad += abs($CTD_w[$ws+$Ci] - $LADCP_w[$ws+$Li]);
+ next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
+ $sad += abs($CTD{w}[$ws+$Ci] - $LADCP_w[$ws+$Li]);
$nad++;
}
next unless ($nad > 0);
@@ -50,29 +58,34 @@
sub lagLADCP2CTD()
{
#------------------------------------------------------------------------
- # find 1st rec & ensemble >=80% down to max depth & make 1st guess at lag
+ # find 1st rec & ensemble >=10% down to max depth & make 1st guess at lag
#------------------------------------------------------------------------
- my($CTD_80pct_down) = 0;
- $CTD_80pct_down++ # "until" formulation allows for missing pressures
- until ($CTD_press[$CTD_80pct_down]-$CTD_press[0] >= 0.8*($CTD_maxpress-$CTD_press[0]));
-
- my($LADCP_80pct_down) = 0;
- $LADCP_80pct_down++
- while ($LADCP{ENSEMBLE}[$LADCP_80pct_down]->{DEPTH} < 0.8*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
-
+ my($first_guess_lag); # in units of CTD records
- my($first_guess_lag) = $LADCP{ENSEMBLE}[$LADCP_80pct_down]->{ELAPSED_TIME} -
- $CTD_80pct_down*$CTD_sampint;
-
- printf(STDERR "\n\t1st guess offset [CTD pressure, LADCP estimated depth] = %ds [%ddbar, %dm]\n",
- $first_guess_lag*$CTD_sampint,$CTD_press[$CTD_80pct_down],$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH})
- if ($opt_d);
- croak("\n$0: cannot determine valid 1st guess offset ($LADCP_80pct_down,$CTD_80pct_down,$CTD_press[$CTD_80pct_down])\n")
- if ($first_guess_lag*$CTD_sampint > 1200);
+ if (defined($opt_i)) {
+ $first_guess_lag = -$opt_i / $CTD{sampint};
+ } else {
+ my($CTD_10pct_down) = 0;
+ $CTD_10pct_down++ # "until" formulation allows for missing pressures
+ until ($CTD{press}[$CTD_10pct_down]-$CTD{press}[0] >= 0.1*($CTD{maxpress}-$CTD{press}[0]));
+
+ my($LADCP_10pct_down) = 0;
+ $LADCP_10pct_down++
+ while ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} < 0.1*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+
+ $first_guess_lag = ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED_TIME} -
+ $CTD_10pct_down*$CTD{sampint}) / $CTD{sampint};
+
+ printf(STDERR "\n\t1st guess offset [CTD pressure, LADCP estimated depth] = %ds [%ddbar, %dm]\n",
+ $first_guess_lag*$CTD{sampint},$CTD{press}[$CTD_10pct_down],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH})
+ if ($opt_d);
+ croak("\n$0: cannot determine valid 1st guess offset ($LADCP_10pct_down,$CTD_10pct_down,$CTD{press}[$CTD_10pct_down])\n")
+ if ($first_guess_lag*$CTD{sampint} > 1200);
+ }
#------------------------------------------------------------------------------------
- # Linearly interpolate LADCP time series onto a new grid with $CTD_sampint resolution
+ # Linearly interpolate LADCP time series onto a new grid with $CTD{sampint} resolution
# ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data
# chunks large
#------------------------------------------------------------------------------------
@@ -80,8 +93,8 @@
$nGaps = 0;
for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {
- while ($r*$CTD_sampint < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
- $LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD_sampint,$ens)
+ while ($r*$CTD{sampint} < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
+ $LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD{sampint},$ens)
unless ($first_guess_lag > $r);
$r++;
}
@@ -96,8 +109,8 @@
# Calculate lags
#----------------------------------------------------------------------
- printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]: ",$opt_w);
- $opt_w = int($opt_w / $CTD_sampint);
+ printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]:",$opt_w);
+ $opt_w = int($opt_w / $CTD{sampint});
#---------------------------------------------------------------
# carry out opt_n lag correlations and keep tally of the results
@@ -105,20 +118,37 @@
my(%nBest);
my($nLags) = 0;
my($lags) = '';
+ my($lastLag) = 9e99; my($nSame) = 1;
for (my($window)=0; $window<$opt_n; $window++) {
my($ws) = $window * int(@LADCP_w/$opt_n); # window start
$ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);
-
$bestLag = bestLag($ws);
if (defined($bestLag)) {
- printf(STDERR "%d ",$bestLag*$CTD_sampint);
- $lags .= sprintf(" %s",$bestLag*$CTD_sampint);
+ if (defined($lastLag) && $bestLag == $lastLag) {
+ $nSame++;
+ } else {
+ printf(STDERR "(x%d)",$nSame)
+ if ($nSame > 1);
+ printf(STDERR " %d",$bestLag*$CTD{sampint});
+ $nSame = 1;
+ $lastLag = $bestLag;
+ }
+ $lags .= sprintf(" %s",$bestLag*$CTD{sampint});
$nBest{$bestLag}++;
$nLags++;
} else {
- printf(STDERR "nan ");
+ if (!defined($lastLag)) {
+ $nSame++;
+ } else {
+ printf(STDERR "(x%d)",$nSame)
+ if ($nSame > 1);
+ printf(STDERR " nan");
+ $nSame = 1;
+ $lastLag = $bestLag;
+ }
}
}
+ printf(STDERR "(x%d)",$nSame) if ($nSame > 1);
&antsAddParams('time_lags',$lags);
#----------------------
@@ -129,17 +159,18 @@
$best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
}
croak("\n$0: cannot determine a valid lag\n")
- unless ($nBest{$best_lag} > 1);
+ unless ($nBest{$best_lag} > $opt_n/3);
print(STDERR "\n\n\t\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
if ($nBest{$best_lag} < $opt_n/2);
if ($nBest{$best_lag} == $nLags) {
- printf(STDERR "\n\t\tunanimous lag = %ds\n",$best_lag*$CTD_sampint);
+ printf(STDERR "\n\t\tunanimous lag = %ds\n",$best_lag*$CTD{sampint});
} else {
- printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD_sampint);
+ printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD{sampint});
}
- return ($first_guess_lag + $best_lag) * $CTD_sampint;
+ &antsAddParams('LADCP_time_lag',($first_guess_lag + $best_lag) * $CTD{sampint});
+ return ($first_guess_lag + $best_lag) * $CTD{sampint};
}
1;