0
|
1 |
#======================================================================
|
|
2 |
# C A L C _ T I M E L A G S . P L
|
|
3 |
# doc: Fri Dec 17 21:59:07 2010
|
|
4 |
# dlm: Tue Dec 21 00:40:33 2010
|
|
5 |
# (c) 2010 A.M. Thurnherr
|
|
6 |
# uE-Info: 117 20 NIL 0 0 72 2 2 4 NIL ofnI
|
|
7 |
#======================================================================
|
|
8 |
|
|
9 |
# HISTORY:
|
|
10 |
# Dec 17, 2010: - created
|
|
11 |
# Dec 18, 2010: - adapted for multi-pass lagging
|
|
12 |
# Dec 20: 2010: - added code to adjust start and end of profile ens
|
|
13 |
# based on extent of CTD profile and guestimated time
|
|
14 |
# ofset
|
|
15 |
|
|
16 |
# TODO:
|
|
17 |
# - better seabed code (from LADCPproc)
|
|
18 |
# - intermediate-step timelagging guess
|
|
19 |
# - flip aliased ensembles
|
|
20 |
|
|
21 |
sub mad_w($$$) # mean absolute deviation
|
|
22 |
{
|
|
23 |
my($fe,$le,$so) = @_; # first/last LADCP ens, CTD scan offset
|
|
24 |
my($sad) = my($n) = 0;
|
|
25 |
|
|
26 |
for (my($e)=$fe; $e<=$le; $e++) {
|
|
27 |
my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5);
|
|
28 |
die("assertion failed\n" .
|
|
29 |
"\ttest: abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2\n" .
|
|
30 |
"\te = $e, s = $s, ensemble = $LADCP{ENSEMBLE}[$e]->{NUMBER}"
|
|
31 |
) unless (abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2);
|
|
32 |
|
|
33 |
my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W} - $CTD{W}[$s+$so];
|
|
34 |
next unless (abs($dw) <= $opt_m);
|
|
35 |
|
|
36 |
$sad += abs($dw);
|
|
37 |
$n++;
|
|
38 |
}
|
|
39 |
return ($n>0) ? $sad/$n : 9e99; # n == 0, e.g. in bottom gap
|
|
40 |
}
|
|
41 |
|
|
42 |
|
|
43 |
sub bestLag($$$$) # find best lag in window
|
|
44 |
{
|
|
45 |
my($fe,$le,$ww,$soi) = @_; # first/last LADCP ens, window width, scan-offset increment
|
|
46 |
|
|
47 |
my($bestso) = 0; # error at first-guess offset
|
|
48 |
my($bestmad) = mad_w($fe,$le,0);
|
|
49 |
|
|
50 |
for (my($dso) = 1; $dso <= int($ww/2/$CTD{DT} + 0.5); $dso+=$soi) {
|
|
51 |
my($mad) = mad_w($fe,$le,-$dso);
|
|
52 |
$bestmad=$mad,$bestso=-$dso if ($mad < $bestmad);
|
|
53 |
$mad = mad_w($fe,$le,$dso);
|
|
54 |
$bestmad=$mad,$bestso=$dso if ($mad < $bestmad);
|
|
55 |
}
|
|
56 |
return ($bestso,$bestmad);
|
|
57 |
}
|
|
58 |
|
|
59 |
#----------------------------------------------------------------------
|
|
60 |
# carry out lag correlations and keep tally of the results
|
|
61 |
# - fist and last 10% of LADCP profile ignored
|
|
62 |
#----------------------------------------------------------------------
|
|
63 |
|
|
64 |
sub calc_lag($$$)
|
|
65 |
{
|
|
66 |
my($n_windows,$w_size,$scan_increment) = @_;
|
|
67 |
|
|
68 |
RETRY:
|
|
69 |
progress("Calculating $n_windows time-lags from ${w_size}s-long windows at %dHz resolution...\n",
|
|
70 |
int(1/$scan_increment/$CTD{DT}+0.5));
|
|
71 |
|
|
72 |
my($approx_CTD_profile_start_ens) =
|
|
73 |
$firstGoodEns + int(($CTD{ELAPSED}[0] - $CTD{TIME_LAG}) / $LADCP{MEAN_DT});
|
|
74 |
my($approx_CTD_profile_end_ens) =
|
|
75 |
$firstGoodEns + int(($CTD{ELAPSED}[$#{$CTD{ELAPSED}}] + $CTD{ELAPSED}[0] - $CTD{TIME_LAG}) / $LADCP{MEAN_DT});
|
|
76 |
|
|
77 |
my($approx_joint_profile_start_ens) = max($firstGoodEns,$approx_CTD_profile_start_ens);
|
|
78 |
my($approx_joint_profile_end_ens) = min($lastGoodEns,$approx_CTD_profile_end_ens);
|
|
79 |
debugmsg("profile start: $firstGoodEns -> $approx_joint_profile_start_ens\n");
|
|
80 |
debugmsg("profile end : $lastGoodEns -> $approx_joint_profile_end_ens\n");
|
|
81 |
|
|
82 |
my($skip_ens) = int(($approx_joint_profile_end_ens - $approx_joint_profile_start_ens) / 10 + 0.5);
|
|
83 |
|
|
84 |
my(%nBest);
|
|
85 |
for (my($wi)=0; $wi<$n_windows; $wi++) {
|
|
86 |
my($fe) = $approx_joint_profile_start_ens + $skip_ens + $wi*int(($approx_joint_profile_end_ens-$approx_joint_profile_start_ens-2*$skip_ens)/$n_windows+0.5);
|
|
87 |
my($so,$mad) = bestLag($fe,$fe+int($w_size/$LADCP{MEAN_DT}+0.5),$w_size,$scan_increment);
|
|
88 |
debugmsg("%.1f cm/s mad(w) at %3d scans offset\n",100*$mad,$so);
|
|
89 |
$nBest{$so}++;
|
|
90 |
}
|
|
91 |
|
|
92 |
my(@best_lag);
|
|
93 |
foreach my $i (keys(%nBest)) {
|
|
94 |
$best_lag[0] = $i if ($nBest{$i} > $nBest{$best_lag[0]});
|
|
95 |
}
|
|
96 |
foreach my $i (keys(%nBest)) {
|
|
97 |
next if ($i == $best_lag[0]);
|
|
98 |
$best_lag[1] = $i if ($nBest{$i} > $nBest{$best_lag[1]});
|
|
99 |
}
|
|
100 |
foreach my $i (keys(%nBest)) {
|
|
101 |
next if ($i == $best_lag[0] || $i == $best_lag[1]);
|
|
102 |
$best_lag[2] = $i if ($nBest{$i} > $nBest{$best_lag[2]});
|
|
103 |
}
|
|
104 |
progress("\t3 most popular offsets: %d (%d%%), %d (%d%%), %d (%d%%)\n",
|
|
105 |
$best_lag[0],int(($nBest{$best_lag[0]}/$n_windows)*100+0.5),
|
|
106 |
$best_lag[1],int(($nBest{$best_lag[1]}/$n_windows)*100+0.5),
|
|
107 |
$best_lag[2],int(($nBest{$best_lag[2]}/$n_windows)*100+0.5));
|
|
108 |
|
|
109 |
if ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} <= 6) {
|
|
110 |
warning(0,"cannot determine valid lag => trying again with doubled window size\n");
|
|
111 |
undef(%nBest);
|
|
112 |
$w_size *= 2;
|
|
113 |
goto RETRY;
|
|
114 |
}
|
|
115 |
|
|
116 |
croak(sprintf("$0: cannot determine a valid lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
|
|
117 |
int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_windows+0.5)))
|
|
118 |
unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_windows);
|
|
119 |
|
|
120 |
my($wmo) = ($nBest{$best_lag[0]}*$best_lag[0] +
|
|
121 |
$nBest{$best_lag[1]}*$best_lag[1] +
|
|
122 |
$nBest{$best_lag[2]}*$best_lag[2]) / ($nBest{$best_lag[0]} +
|
|
123 |
$nBest{$best_lag[1]} +
|
|
124 |
$nBest{$best_lag[2]});
|
|
125 |
progress("\tweighted-mean offset = %.1f scans\n",$wmo);
|
|
126 |
|
|
127 |
if ($wmo > 0.9*$w_size/2/$CTD{DT}) {
|
|
128 |
warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
|
|
129 |
$CTD{TIME_LAG} += $w_size/2;
|
|
130 |
undef(%nBest);
|
|
131 |
goto RETRY;
|
|
132 |
}
|
|
133 |
if (-$wmo > 0.9*$w_size/2/$CTD{DT}) {
|
|
134 |
warning(0,"lag too close to the edge of the window --- trying again after adjusting the guestimated offset\n");
|
|
135 |
$CTD{TIME_LAG} -= $w_size/2;
|
|
136 |
undef(%nBest);
|
|
137 |
goto RETRY;
|
|
138 |
}
|
|
139 |
|
|
140 |
return $CTD{TIME_LAG}+$wmo*$CTD{DT};
|
|
141 |
}
|
|
142 |
|
|
143 |
|
|
144 |
1;
|