0
|
1 |
#!/usr/bin/perl
|
|
2 |
#======================================================================
|
|
3 |
# L I S T B T
|
|
4 |
# doc: Sat Jan 18 18:41:49 2003
|
|
5 |
# dlm: Sun May 23 16:33:33 2010
|
|
6 |
# (c) 2003 A.M. Thurnherr
|
|
7 |
# uE-Info: 527 42 NIL 0 0 72 11 2 4 NIL ofnI
|
|
8 |
#======================================================================
|
|
9 |
|
|
10 |
# Extract Bottom-Track Data
|
|
11 |
|
|
12 |
# NOTE: NO SOUND-SPEED CORRECTION APPLIED YET!!!
|
|
13 |
|
|
14 |
# HISTORY:
|
|
15 |
# Jan 18, 2003: - created
|
|
16 |
# Jan 23, 2003: - added magnetic declination
|
|
17 |
# Jan 25, 2003: - continued construction
|
|
18 |
# Feb 11, 2003: - finally made it work
|
|
19 |
# Feb 12, 2003: - added default profile output
|
|
20 |
# Feb 13, 2003: - corrected raw output
|
|
21 |
# Feb 14, 2003: - added errors if instrument-BT filters are more strict
|
|
22 |
# than command-line values
|
|
23 |
# Feb 18, 2003: - removed -d dependency on -W
|
|
24 |
# Mar 3, 2003: - added -C)ompass correction
|
|
25 |
# Mar 10, 2003: - added -f)orce to allow visbeck-style post processing
|
|
26 |
# Mar 16, 2003: - added range comment
|
|
27 |
# Feb 26, 2004: - added Earth-coordinate support
|
|
28 |
# Feb 27, 2004: - made water-track calculation conditional (-E || -B)
|
|
29 |
# Mar 9, 2004: - added magnetic_declination to %PARAMs
|
|
30 |
# Apr 1, 2004: - added CTD u/v stats to %PARAMs
|
|
31 |
# Apr 2, 2004: - added CTD_msf (mean square fluctuation) stat
|
|
32 |
# Apr 3, 2004: - BUG: CTD vels were repeated for stats
|
|
33 |
# - removed non-ANTS option
|
|
34 |
# Nov 8, 2005: - UNIXTIME => UNIX_TIME
|
|
35 |
# - adapted to new binary read library
|
|
36 |
# - output editing statistics
|
|
37 |
# Aug 15, 2006: - added -b
|
|
38 |
# Aug 25, 2006: - fiddled
|
|
39 |
# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
|
|
40 |
# Nov 1, 2008: - BUG: sig(u) was reported instead of sig(v)
|
|
41 |
# Jul 30, 2009: - NaN => nan
|
|
42 |
|
|
43 |
# NOTES:
|
|
44 |
# - the RDI BT data contains ranges that are greater than the
|
|
45 |
# WT ping ranges. I don't know if those data are valid!
|
|
46 |
# - there is a fair bit of heuristic used, especially in the
|
|
47 |
# reference-layer calculation
|
|
48 |
# - depth-correction (-m) is highly recommended because it allows
|
|
49 |
# much better bad-BT detection and it is required for a valid
|
|
50 |
# comparison with LADCP profiles
|
|
51 |
# - the criterion for bottom-interference of the water-track data
|
|
52 |
# is derived from Firing's [merge.c] (adding 1.5 bin lengths to
|
|
53 |
# the calculated range), modified by taking the real beam angle
|
|
54 |
# into account.
|
|
55 |
# - from the RDI manuals it is not entirely clear whether the BT range
|
|
56 |
# is given in vertical or in along-beam meters; comparison with the
|
|
57 |
# WT range (calculated from the bin with the maximum echo amplitude)
|
|
58 |
# shows that vertical meters are used
|
|
59 |
|
|
60 |
# NOTES on quality checks:
|
|
61 |
# -a minimum BT amplitude; setting this to 50 (RDI default is 30)
|
|
62 |
# reduces the vertical range over which the bottom is detected but
|
|
63 |
# not the quality of the bottom track data; therefore, this should
|
|
64 |
# probably not be used.
|
|
65 |
# -c minimum BT correlation; the RDI default for this parameter is 220,
|
|
66 |
# which seems to work fine.
|
|
67 |
# -e max error velocity (BT & WT); this is primarily used for detecting
|
|
68 |
# good BT data, i.e. it should be set to a small value (Firing uses
|
|
69 |
# 0.1m/s in merge); if too small a value is chosen too many good
|
|
70 |
# data are discarded; note that the same error-velocity criterion
|
|
71 |
# is used to separate good from bad data when mean profiles are
|
|
72 |
# constructed.
|
|
73 |
# -w max difference between reference-layer w and BT w; this is a
|
|
74 |
# powerful criterion for determining good BT data; I like a value of
|
|
75 |
# 0.03 m/s.
|
|
76 |
# -d when the depth is corrected (-m) the...
|
|
77 |
|
|
78 |
$0 =~ m{(.*)/[^/]+};
|
|
79 |
require "$1/RDI_BB_Read.pl";
|
|
80 |
require "$1/RDI_Coords.pl";
|
|
81 |
require "$1/RDI_Utils.pl";
|
|
82 |
use Getopt::Std;
|
|
83 |
|
|
84 |
$USAGE = "$0 @ARGV";
|
|
85 |
die("Usage: $0 " .
|
|
86 |
"[use -b)ins <1st,last>] " .
|
|
87 |
"[write -R)aw data] [write -B)T data] " .
|
|
88 |
"[write -E)nsembles <pref>] [-F)ilter ensembles <script>] " .
|
|
89 |
"[-C)ompass correction <amp/phase/bias>] " .
|
|
90 |
"[-w) <max-diff|0.03>] [-a)mp <min|30>] [-e)rr-vel <max|0.05>] " .
|
|
91 |
"[-c)orrelation <min|220>] " .
|
|
92 |
"[-W)ater <depth> [allowed -d)epth-diff <maxdiff|20>]] " .
|
|
93 |
"[-f)orce (no setup tests)] " .
|
|
94 |
"[-M)agnetic <declination>] " .
|
|
95 |
"<RDI file>\n")
|
|
96 |
unless (&getopts("BC:E:F:M:RW:a:b:c:d:e:fw:") && @ARGV == 1);
|
|
97 |
|
|
98 |
print(STDERR "WARNING: magnetic declination not set!\n")
|
|
99 |
unless defined($opt_M);
|
|
100 |
|
|
101 |
$opt_c = 220 unless defined($opt_c); # defaults
|
|
102 |
$opt_a = 30 unless defined($opt_a);
|
|
103 |
$opt_e = 0.05 unless defined($opt_e);
|
|
104 |
$opt_w = 0.03 unless defined($opt_w);
|
|
105 |
$opt_d = 20 unless defined($opt_d);
|
|
106 |
|
|
107 |
if (defined($opt_C)) { # compass correction
|
|
108 |
($CC_amp,$CC_phase,$CC_bias) = split('/',$opt_C);
|
|
109 |
die("$0: can't decode -C$opt_C\n")
|
|
110 |
unless defined($CC_bias);
|
|
111 |
}
|
|
112 |
|
|
113 |
unless ($opt_f) { # check BT setup
|
|
114 |
readHeader($ARGV[0],\%dta);
|
|
115 |
die("$0: minimum instrument BT correlation ($dta{BT_MIN_CORRELATION}) " .
|
|
116 |
"too large for selected criterion (-c $opt_c) --- use -f to override\n")
|
|
117 |
if ($dta{BT_MIN_CORRELATION} > $opt_c);
|
|
118 |
die("$0: minimum instrument BT amplitude ($dta{BT_MIN_EVAL_AMPLITUDE}) " .
|
|
119 |
"too large for selected criterion (-a $opt_a) --- use -f to override\n")
|
|
120 |
if ($dta{BT_MIN_EVAL_AMPLITUDE} > $opt_a);
|
|
121 |
die("$0: maximum instrument BT error velocity ($dta{BT_MAX_ERROR_VELOCITY}) " .
|
|
122 |
"too small for selected criterion (-e $opt_e) --- use -f to override\n")
|
|
123 |
if ($dta{BT_MAX_ERROR_VELOCITY} < $opt_e);
|
|
124 |
}
|
|
125 |
|
|
126 |
require $opt_F if defined($opt_F); # load filter
|
|
127 |
|
|
128 |
print(STDERR "reading $ARGV[0]...");
|
|
129 |
readData($ARGV[0],\%dta); # read data
|
|
130 |
print(STDERR "done\n");
|
|
131 |
|
|
132 |
$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
|
|
133 |
ensure_BT_RANGE(\%dta); # make sure they're there
|
|
134 |
|
|
135 |
$firstBin = $lastBin = '*'; # bins to use
|
|
136 |
($firstBin,$lastBin) = split(',',$opt_b)
|
|
137 |
if defined($opt_b);
|
|
138 |
$firstBin = 1 if ($firstBin eq '*');
|
|
139 |
$lastBin = $dta{N_BINS} if ($lastBin eq '*');
|
|
140 |
$firstBin--; $lastBin--;
|
|
141 |
die("$ARGV[0]: not enough bins for ref layer\n")
|
|
142 |
unless ($lastBin-$firstBin >= 6);
|
|
143 |
|
|
144 |
if ($dta{BEAM_COORDINATES}) { # coords used
|
|
145 |
$beamCoords = 1;
|
|
146 |
} elsif (!$dta{EARTH_COORDINATES}) {
|
|
147 |
die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
|
|
148 |
}
|
|
149 |
|
|
150 |
#======================================================================
|
|
151 |
# Calculate reference-layer w
|
|
152 |
#======================================================================
|
|
153 |
|
|
154 |
sub w($)
|
|
155 |
{
|
|
156 |
my($ens) = @_;
|
|
157 |
my($i,$n,@v,$w);
|
|
158 |
|
|
159 |
for (my($b)=$firstBin; $b<=$lastBin; $b++) {
|
|
160 |
if ($beamCoords) {
|
|
161 |
next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][0] < 100 ||
|
|
162 |
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][1] < 100 ||
|
|
163 |
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][2] < 100 ||
|
|
164 |
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][3] < 100);
|
|
165 |
@v = velInstrumentToEarth(\%dta,$ens,
|
|
166 |
velBeamToInstrument(\%dta,
|
|
167 |
@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$b]}));
|
|
168 |
} else {
|
|
169 |
next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][0] > 0 ||
|
|
170 |
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][1] > 0 ||
|
|
171 |
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][2] > 0 ||
|
|
172 |
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$b][3] < 100);
|
|
173 |
@v = velApplyHdgBias(\%dta,$ens,
|
|
174 |
@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$b]});
|
|
175 |
}
|
|
176 |
next unless (defined($v[3]) && abs($v[3]) <= $opt_e);
|
|
177 |
$w += $v[2]; $n++;
|
|
178 |
}
|
|
179 |
# printf(STDERR "$ens $n %.3f\n",$n>=1?$w/$n:-999);
|
|
180 |
return $n>=2 ? $w/$n : undef;
|
|
181 |
}
|
|
182 |
|
|
183 |
#======================================================================
|
|
184 |
# Dump raw BT data from one ensemble
|
|
185 |
#======================================================================
|
|
186 |
|
|
187 |
sub dumpRaw($)
|
|
188 |
{
|
|
189 |
my($e) = @_;
|
|
190 |
|
|
191 |
unless ($headerDone) {
|
|
192 |
print("#ANTS# [] $USAGE\n");
|
|
193 |
print("#ANTS#PARAMS# RDI_file{$ARGV[0]}\n");
|
|
194 |
print("#ANTS#FIELDS# {ens} {range1} {range2} {range3} {range4} " .
|
|
195 |
"{beamvel1} {beamvel2} {beamvel3} {beamvel4} {cor1} " .
|
|
196 |
"{cor2} {cor3} {cor4} {amp1} {amp2} {amp3} {amp4}\n");
|
|
197 |
$headerDone = 1;
|
|
198 |
}
|
|
199 |
|
|
200 |
printf("%d %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d\n",
|
|
201 |
$dta{ENSEMBLE}[$e]->{NUMBER},
|
|
202 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[0],
|
|
203 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[1],
|
|
204 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[2],
|
|
205 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[3],
|
|
206 |
$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0],
|
|
207 |
$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1],
|
|
208 |
$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2],
|
|
209 |
$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3],
|
|
210 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0],
|
|
211 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1],
|
|
212 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2],
|
|
213 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3],
|
|
214 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[0],
|
|
215 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[1],
|
|
216 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[2],
|
|
217 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[3]
|
|
218 |
);
|
|
219 |
}
|
|
220 |
|
|
221 |
#======================================================================
|
|
222 |
# Dump processed BT data from one ensemble
|
|
223 |
#======================================================================
|
|
224 |
|
|
225 |
sub dumpBT($)
|
|
226 |
{
|
|
227 |
my($e) = @_;
|
|
228 |
|
|
229 |
unless ($headerDone) {
|
|
230 |
print("#ANTS# [] $USAGE\n");
|
|
231 |
printf("#ANTS#PARAMS# RDI_file{$ARGV[0]} bottom_time{%.1f}\n",
|
|
232 |
$dta{ENSEMBLE}[$maxz_e]->{ELAPSED});
|
|
233 |
print("#ANTS#FIELDS# {ens} {unix_time} {time} {depth} {BT_range} " .
|
|
234 |
"{WT_range} {u} {v} {w} {e} {w_ref} {corr} {amp}\n");
|
|
235 |
$headerDone = 1;
|
|
236 |
}
|
|
237 |
|
|
238 |
printf("%d %.2f %.2f %.1f %.1f %.1f %.4f %.4f %.4f %.4f %.4f %.1f %.1f\n",
|
|
239 |
$dta{ENSEMBLE}[$e]->{NUMBER},
|
|
240 |
$dta{ENSEMBLE}[$e]->{UNIX_TIME},
|
|
241 |
$dta{ENSEMBLE}[$e]->{ELAPSED},
|
|
242 |
$dta{ENSEMBLE}[$e]->{DEPTH},
|
|
243 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE},
|
|
244 |
$dta{ENSEMBLE}[$e]->{WT_MEAN_RANGE},
|
|
245 |
@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}},
|
|
246 |
$dta{ENSEMBLE}[$e]->{W_REF},
|
|
247 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_CORRELATION},
|
|
248 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_EVAL_AMPLITUDE}
|
|
249 |
);
|
|
250 |
}
|
|
251 |
|
|
252 |
#======================================================================
|
|
253 |
# Dump a single ensemble with valid BT data to separate file
|
|
254 |
#======================================================================
|
|
255 |
|
|
256 |
sub dumpEns(@) # write profile
|
|
257 |
{
|
|
258 |
my($e) = @_;
|
|
259 |
my($b,$i);
|
|
260 |
|
|
261 |
open(P,">$opt_E.$e") || die("$opt_E.$e: $!\n");
|
|
262 |
print(P "#ANTS#PARAMS# " .
|
|
263 |
"depth{$dta{ENSEMBLE}[$e]->{DEPTH}} " .
|
|
264 |
"range{$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}} " .
|
|
265 |
"wt_range{$dta{ENSEMBLE}[$e]->{WT_MEAN_RANGE}} " .
|
|
266 |
"w_ref{$dta{ENSEMBLE}[$e]->{W_REF}} " .
|
|
267 |
"BT_u{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0]} " .
|
|
268 |
"BT_v{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1]} " .
|
|
269 |
"BT_w{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]} " .
|
|
270 |
"BT_e{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3]} " .
|
|
271 |
"BT_cor1{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0]} " .
|
|
272 |
"BT_cor2{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1]} " .
|
|
273 |
"BT_cor3{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2]} " .
|
|
274 |
"BT_cor4{$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3]} " .
|
|
275 |
"BT_amp1{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[0]} " .
|
|
276 |
"BT_amp2{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[1]} " .
|
|
277 |
"BT_amp3{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[2]} " .
|
|
278 |
"BT_amp4{$dta{ENSEMBLE}[$e]->{BT_AMPLITUDE}[3]} " .
|
|
279 |
"BTFWT_u{$dta{ENSEMBLE}[$e]->{BTFWT_VELOCITY}[0]} " .
|
|
280 |
"BTFWT_v{$dta{ENSEMBLE}[$e]->{BTFWT_VELOCITY}[1]} " .
|
|
281 |
"BTFWT_w{$dta{ENSEMBLE}[$e]->{BTFWT_VELOCITY}[2]} " .
|
|
282 |
"\n"
|
|
283 |
);
|
|
284 |
print(P "#ANTS#FIELDS# " .
|
|
285 |
"{depth} {hab} {u} {v} {w} {e} {cor1} {cor2} {cor3} {cor4} " .
|
|
286 |
"{amp1} {amp2} {amp3} {amp4} {pcg1} {pcg2} {pcg3} {pcg4}\n"
|
|
287 |
);
|
|
288 |
|
|
289 |
my($slc) = (1-cos(rad($dta{BEAM_ANGLE})))*$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}
|
|
290 |
+ 1.5*$dta{BIN_LENGTH}; # side-lobe contamination
|
|
291 |
for ($b=$firstBin; $b<=$lastBin; $b++) {
|
|
292 |
next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]) &&
|
|
293 |
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]) &&
|
|
294 |
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]) &&
|
|
295 |
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]));
|
|
296 |
my($dz) = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
|
|
297 |
last if ($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}-$dz <= $slc);
|
|
298 |
my(@v) = $beamCoords
|
|
299 |
? velInstrumentToEarth(\%dta,$e,
|
|
300 |
velBeamToInstrument(\%dta,
|
|
301 |
@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}))
|
|
302 |
: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
|
|
303 |
next unless defined($v[0]);
|
|
304 |
next if (abs($v[3]) > $opt_e ||
|
|
305 |
abs($v[2]-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]) > 0.1);
|
|
306 |
$v[0] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0];
|
|
307 |
$v[1] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1];
|
|
308 |
my(@out) = (
|
|
309 |
$dta{ENSEMBLE}[$e]->{DEPTH}+$dz,
|
|
310 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}-$dz,
|
|
311 |
$v[0],$v[1],$v[2],$v[3],
|
|
312 |
@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]},
|
|
313 |
@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]},
|
|
314 |
@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}
|
|
315 |
);
|
|
316 |
for ($i=0; $i<17; $i++) { $out[$i] = nan unless defined($out[$i]); }
|
|
317 |
print(P "@out\n");
|
|
318 |
}
|
|
319 |
close(P);
|
|
320 |
}
|
|
321 |
|
|
322 |
#======================================================================
|
|
323 |
# Add Ensemble With Valid BT Data to Profile
|
|
324 |
#======================================================================
|
|
325 |
|
|
326 |
sub binEns($)
|
|
327 |
{
|
|
328 |
my($e) = @_;
|
|
329 |
my($slc) = (1-cos(rad($dta{BEAM_ANGLE})))*$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}
|
|
330 |
+ 1.5*$dta{BIN_LENGTH}; # side-lobe contamination
|
|
331 |
for (my($b)=$firstBin; $b<=$lastBin; $b++) {
|
|
332 |
next unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]) &&
|
|
333 |
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]) &&
|
|
334 |
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]) &&
|
|
335 |
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]));
|
|
336 |
my($dz) = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
|
|
337 |
last if ($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE}-$dz <= $slc);
|
|
338 |
my(@v) = $beamCoords
|
|
339 |
? velInstrumentToEarth(\%dta,$e,
|
|
340 |
velBeamToInstrument(\%dta,
|
|
341 |
@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}))
|
|
342 |
: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});
|
|
343 |
next unless defined($v[0]);
|
|
344 |
next if (abs($v[3]) > $opt_e ||
|
|
345 |
abs($v[2]-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2]) > 0.1);
|
|
346 |
|
|
347 |
$v[0] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0];
|
|
348 |
$v[1] -= $dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1];
|
|
349 |
|
|
350 |
my($bin) = int(($dta{ENSEMBLE}[$e]->{DEPTH}+$dz) / $dta{BIN_LENGTH});
|
|
351 |
$minBin = $bin unless ($bin >= $minBin);
|
|
352 |
$maxBin = $bin unless ($bin <= $maxBin);
|
|
353 |
if (defined($BTn[$bin])) {
|
|
354 |
my($f1) = $BTn[$bin] / ($BTn[$bin]+1);
|
|
355 |
my($f2) = ($BTn[$bin]-1) / $BTn[$bin];
|
|
356 |
$BTsigu[$bin] =
|
|
357 |
$f2*$BTsigu[$bin] + $f1*($v[0]-$BTu[$bin])**2/$BTn[$bin];
|
|
358 |
$BTsigv[$bin] =
|
|
359 |
$f2*$BTsigv[$bin] + $f1*($v[1]-$BTv[$bin])**2/$BTn[$bin];
|
|
360 |
$BTn[$bin]++;
|
|
361 |
$BTu[$bin] = $f1*$BTu[$bin] + $v[0]/$BTn[$bin];
|
|
362 |
$BTv[$bin] = $f1*$BTv[$bin] + $v[1]/$BTn[$bin];
|
|
363 |
} else {
|
|
364 |
$BTu[$bin] = $v[0];
|
|
365 |
$BTv[$bin] = $v[1];
|
|
366 |
$BTn[$bin] = 1;
|
|
367 |
}
|
|
368 |
}
|
|
369 |
}
|
|
370 |
|
|
371 |
#======================================================================
|
|
372 |
# Output Bottom-Referenced Profile
|
|
373 |
#======================================================================
|
|
374 |
|
|
375 |
sub dumpProf($$$)
|
|
376 |
{
|
|
377 |
my($db,$wd,$md) = @_;
|
|
378 |
my(@sum,@mean);
|
|
379 |
|
|
380 |
for (my($i)=0; $i<=$#listCTDu; $i++) { # CTD vel mean
|
|
381 |
$sum[0] += $listCTDu[$i];
|
|
382 |
$sum[1] += $listCTDv[$i];
|
|
383 |
}
|
|
384 |
@mean = ($sum[0]/@listCTDu,$sum[1]/@listCTDv);
|
|
385 |
@sum = (0,0); # stddev
|
|
386 |
for (my($i)=0; $i<=$#listCTDu; $i++) {
|
|
387 |
$sum[0] += ($listCTDu[$i]-$mean[0])**2;
|
|
388 |
$sum[1] += ($listCTDv[$i]-$mean[1])**2;
|
|
389 |
}
|
|
390 |
@sigma = ($sum[0]/sqrt($#listCTDu),$sum[1]/sqrt($#listCTDv));
|
|
391 |
@sum = (0,0); # mean speed fluct
|
|
392 |
for (my($i)=1; $i<=$#listCTDu; $i++) { # also: list for median
|
|
393 |
push(@cfluc,sqrt(($listCTDu[$i]-$listCTDu[$i-1])**2 +
|
|
394 |
($listCTDv[$i]-$listCTDv[$i-1])**2));
|
|
395 |
$sum[0] += $cfluc[$#cfluc];
|
|
396 |
}
|
|
397 |
|
|
398 |
printf("#ANTS#PARAMS# LADCP_depth_bias{%.1f} water_depth{%.1f} magnetic_declination{%.1f}\n",
|
|
399 |
$db,$wd,$md);
|
|
400 |
printf("#ANTS#PARAMS# CTD_u{%.3f} CTD_v{%.3f} CTD_sig_u{%.3f} CTD_sig_v{%.3f} CTD_mean_cfluc{%.4f} CTD_median_cfluc{%.4f}\n",
|
|
401 |
@mean,@sigma,$sum[0]/$#listCTDu,(sort{$a<=>$b}@cfluc)[@cfluc/2]);
|
|
402 |
printf("#ANTS#PARAMS# good_ens{$good}\n");
|
|
403 |
print("#ANTS#FIELDS# {depth} {u} {v} {sig_u} {sig_v} {n_data}\n");
|
|
404 |
|
|
405 |
for (my($bin)=$minBin; $bin<=$maxBin; $bin++) {
|
|
406 |
next unless defined($BTu[$bin]);
|
|
407 |
printf("%d %.3f %.3f %.3f %.3f %d\n",
|
|
408 |
($bin+0.5)*$dta{BIN_LENGTH},
|
|
409 |
$BTu[$bin],$BTv[$bin],
|
|
410 |
sqrt($BTsigu[$bin]),sqrt($BTsigv[$bin]),
|
|
411 |
$BTn[$bin]);
|
|
412 |
}
|
|
413 |
}
|
|
414 |
|
|
415 |
#======================================================================
|
|
416 |
# STEP 1: Calculate Depth (integrate w)
|
|
417 |
#======================================================================
|
|
418 |
|
|
419 |
for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
|
|
420 |
checkEnsemble(\%dta,$e);
|
|
421 |
$dta{ENSEMBLE}[$e]->{W_REF} = w($e);
|
|
422 |
next unless (defined($start_e) ||
|
|
423 |
defined($dta{ENSEMBLE}[$e]->{W_REF}));
|
|
424 |
$start_e = $e unless defined($start_e);
|
|
425 |
$end_e = $e if defined($dta{ENSEMBLE}[$e]->{W_REF});
|
|
426 |
$lasttime = $curtime;
|
|
427 |
$curtime = $dta{ENSEMBLE}[$e]->{UNIX_TIME};
|
|
428 |
$dta{ENSEMBLE}[$e]->{ELAPSED} =
|
|
429 |
$curtime - $dta{ENSEMBLE}[$start_e]->{UNIX_TIME};
|
|
430 |
filterEnsemble(\%dta,$e)
|
|
431 |
if (defined($opt_F) &&
|
|
432 |
$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][3] > 0);
|
|
433 |
$z += $dta{ENSEMBLE}[$e]->{W_REF} *
|
|
434 |
(defined($lasttime) ? ($curtime - $lasttime) : 0);
|
|
435 |
$maxz_e=$e,$maxz = $z unless ($z < $maxz);
|
|
436 |
$dta{ENSEMBLE}[$e]->{DEPTH} = $z;
|
|
437 |
}
|
|
438 |
|
|
439 |
unless ($opt_R) {
|
|
440 |
($w_depth,$swd) = find_seabed(\%dta,$maxz_e,$beamCoords);
|
|
441 |
die("$0: can't determine water depth (sigma = $swd)\n")
|
|
442 |
unless (defined($w_depth) && $swd < 10);
|
|
443 |
|
|
444 |
if (defined($opt_W)) { # adjust depth
|
|
445 |
$zbias = $w_depth - $opt_W; $w_depth = $opt_W;
|
|
446 |
for ($e=$start_e; $e<=$end_e; $e++) {
|
|
447 |
$dta{ENSEMBLE}[$e]->{DEPTH} -= $zbias;
|
|
448 |
}
|
|
449 |
}
|
|
450 |
}
|
|
451 |
|
|
452 |
# print(STDERR "maxz = $maxz, w_depth = $w_depth\n");
|
|
453 |
|
|
454 |
#======================================================================
|
|
455 |
# STEP 2: Process BT Data
|
|
456 |
#======================================================================
|
|
457 |
|
|
458 |
for ($e=$start_e; $e<=$end_e; $e++) {
|
|
459 |
next unless ($dta{ENSEMBLE}[$e]->{BT_RANGE}[0] && # BT data available
|
|
460 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[1] &&
|
|
461 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[2] &&
|
|
462 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[3]);
|
|
463 |
# die("$0: don't know what to do with non-zero %-good and " .
|
|
464 |
# "signal-strength values at ensemble " .
|
|
465 |
# "#$dta{ENSEMBLE}[$e]->{NUMBER}\n")
|
|
466 |
# if ($dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[0] ||
|
|
467 |
# $dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[1] ||
|
|
468 |
# $dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[2] ||
|
|
469 |
# $dta{ENSEMBLE}[$e]->{BT_PERCENT_GOOD}[3] ||
|
|
470 |
# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[0] ||
|
|
471 |
# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[1] ||
|
|
472 |
# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[2] ||
|
|
473 |
# $dta{ENSEMBLE}[$e]->{BT_SIGNAL_STRENGHT}[3]);
|
|
474 |
|
|
475 |
if ($opt_R) { # dump raw data
|
|
476 |
dumpRaw($e);
|
|
477 |
next;
|
|
478 |
}
|
|
479 |
|
|
480 |
$dta{ENSEMBLE}[$e]->{HEADING} -= # compass correction
|
|
481 |
$CC_amp * sin(rad($dta{ENSEMBLE}[$e]->{HEADING} - $CC_phase))
|
|
482 |
+ $CC_bias
|
|
483 |
if defined($opt_C);
|
|
484 |
|
|
485 |
@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}} = $beamCoords # xform BT vel
|
|
486 |
? velInstrumentToEarth(\%dta,$e,
|
|
487 |
velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}}))
|
|
488 |
: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}});
|
|
489 |
|
|
490 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE} = # mean vals
|
|
491 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[0]/4 +
|
|
492 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[1]/4 +
|
|
493 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[2]/4 +
|
|
494 |
$dta{ENSEMBLE}[$e]->{BT_RANGE}[3]/4;
|
|
495 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_CORRELATION} =
|
|
496 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0]/4 +
|
|
497 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1]/4 +
|
|
498 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2]/4 +
|
|
499 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3]/4;
|
|
500 |
$dta{ENSEMBLE}[$e]->{BT_MEAN_EVAL_AMPLITUDE} =
|
|
501 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[0]/4 +
|
|
502 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[1]/4 +
|
|
503 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[2]/4 +
|
|
504 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[3]/4;
|
|
505 |
|
|
506 |
# next # could add this
|
|
507 |
# if ($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE} < 50);
|
|
508 |
|
|
509 |
$bad_amp++,next
|
|
510 |
if ($dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[0] < $opt_a ||
|
|
511 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[1] < $opt_a ||
|
|
512 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[2] < $opt_a ||
|
|
513 |
$dta{ENSEMBLE}[$e]->{BT_EVAL_AMPLITUDE}[3] < $opt_a);
|
|
514 |
$bad_corr++,next
|
|
515 |
if ($dta{ENSEMBLE}[$e]->{BT_CORRELATION}[0] < $opt_c ||
|
|
516 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[1] < $opt_c ||
|
|
517 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[2] < $opt_c ||
|
|
518 |
$dta{ENSEMBLE}[$e]->{BT_CORRELATION}[3] < $opt_c);
|
|
519 |
|
|
520 |
$bad_w_ref++,next # quality checks
|
|
521 |
if (abs($dta{ENSEMBLE}[$e]->{BT_VELOCITY}[2] -
|
|
522 |
$dta{ENSEMBLE}[$e]->{W_REF}) > $opt_w);
|
|
523 |
$bad_e_vel++,next
|
|
524 |
if (abs($dta{ENSEMBLE}[$e]->{BT_VELOCITY}[3]) > $opt_e);
|
|
525 |
$bad_depth++,next
|
|
526 |
if (abs($dta{ENSEMBLE}[$e]->{BT_MEAN_RANGE} +
|
|
527 |
$dta{ENSEMBLE}[$e]->{DEPTH} - $w_depth) > $opt_d);
|
|
528 |
|
|
529 |
$good++;
|
|
530 |
push(@listCTDu,-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[0]);
|
|
531 |
push(@listCTDv,-$dta{ENSEMBLE}[$e]->{BT_VELOCITY}[1]);
|
|
532 |
|
|
533 |
if ($opt_E || $opt_B) {
|
|
534 |
my(@maxamp) = (0,0,0,0); # water-track range
|
|
535 |
my(@btm_e) = (0,0,0,0);
|
|
536 |
for ($b=$firstBin; $b<=$lastBin; $b++) {
|
|
537 |
for ($i=0; $i<4; $i++) {
|
|
538 |
if ($dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][$i] > $maxamp[$i]) {
|
|
539 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[$i] = $b;
|
|
540 |
$maxamp[$i] = $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][$i];
|
|
541 |
$btm_e[$i] = $e;
|
|
542 |
}
|
|
543 |
}
|
|
544 |
}
|
|
545 |
for ($i=0; $i<4; $i++) {
|
|
546 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[$i] *= $dta{BIN_LENGTH};
|
|
547 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[$i] += $dta{DISTANCE_TO_BIN1_CENTER};
|
|
548 |
}
|
|
549 |
$dta{ENSEMBLE}[$e]->{WT_MEAN_RANGE} =
|
|
550 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[0]/4 +
|
|
551 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[1]/4 +
|
|
552 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[2]/4 +
|
|
553 |
$dta{ENSEMBLE}[$e]->{WT_RANGE}[3]/4;
|
|
554 |
|
|
555 |
my($btm_e) = int($btm_e[0]/4+$btm_e[1]/4+$btm_e[2]/4+$btm_e[3]/4+0.5);
|
|
556 |
@{$dta{ENSEMBLE}[$btm_e]->{BTFWT_VELOCITY}} = $beamCoords # BT from WT
|
|
557 |
? velInstrumentToEarth(\%dta,$e,
|
|
558 |
velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}}))
|
|
559 |
: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BT_VELOCITY}});
|
|
560 |
}
|
|
561 |
|
|
562 |
dumpEns($e) if defined($opt_E); # output BT profiles
|
|
563 |
if ($opt_B) { dumpBT($e); }
|
|
564 |
else { binEns($e); }
|
|
565 |
}
|
|
566 |
|
|
567 |
filterEnsembleStats() if defined($opt_F);
|
|
568 |
exit(0) if ($opt_R);
|
|
569 |
|
|
570 |
printf(STDERR "%5d BT records removed due to bad w\n",$bad_w_ref)
|
|
571 |
if defined($bad_w_ref);
|
|
572 |
printf(STDERR "%5d BT records removed due to bad err vel\n",$bad_e_vel)
|
|
573 |
if defined($bad_e_vel);
|
|
574 |
printf(STDERR "%5d BT records removed due to bad echo amplitude\n",$bad_amp)
|
|
575 |
if defined($bad_amp);
|
|
576 |
printf(STDERR "%5d BT records removed due to bad correlation\n",$bad_corr)
|
|
577 |
if defined($bad_corr);
|
|
578 |
printf(STDERR "%5d BT records removed due to bad depth\n",$bad_depth)
|
|
579 |
if defined($bad_depth);
|
|
580 |
|
|
581 |
die("$0: no good BT data\n") unless ($good);
|
|
582 |
|
|
583 |
printf(STDERR "\n%5d BT records remaining\n",$good);
|
|
584 |
|
|
585 |
dumpProf($zbias,$w_depth,-$dta{HEADING_BIAS})
|
|
586 |
unless ($opt_B);
|
|
587 |
|
|
588 |
exit(0);
|
|
589 |
|