|
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 |