0
|
1 |
#!/usr/bin/perl
|
|
2 |
#======================================================================
|
18
|
3 |
# B E A M S T A T S
|
0
|
4 |
# doc: Fri Aug 25 15:57:05 2006
|
18
|
5 |
# dlm: Tue Mar 4 13:09:14 2014
|
0
|
6 |
# (c) 2006 A.M. Thurnherr
|
18
|
7 |
# uE-Info: 33 64 NIL 0 0 72 0 2 4 NIL ofnI
|
0
|
8 |
#======================================================================
|
|
9 |
|
|
10 |
# Split data file into per-bin time series.
|
|
11 |
|
|
12 |
# HISTORY:
|
|
13 |
# Aug 25, 2006: - created from [listEns]
|
|
14 |
# Aug 26, 2006: - added -M)agdecl
|
|
15 |
# - changed -b to -f
|
|
16 |
# Aug 27, 2006: - added %bin
|
|
17 |
# Aug 28, 2006: - BUG: options were confused
|
|
18 |
# Jan 4, 2007: - improved usage message for -a
|
|
19 |
# - added %mag_decl
|
|
20 |
# - BUG: roundoff error in %pct_good_vels
|
|
21 |
# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl]
|
|
22 |
# Feb 7, 2008: - added sound-speed correction
|
|
23 |
# - enabled 3-beam solutions
|
|
24 |
# Feb 8, 2008: - added -d)iscard <beam>
|
|
25 |
# - added -b)eam coordinate output
|
|
26 |
# Feb 12, 2008: - modified 3-beam output
|
|
27 |
# - added -p)ct_good <min>
|
|
28 |
# Feb 13, 2008: - various improvements
|
|
29 |
# Feb 19, 2008: - BUG: division by zero
|
|
30 |
# BUG: min() did not work with 1st elt undef
|
|
31 |
# Feb 21, 2008: - BUG: had forgotten to undo debugging changes
|
|
32 |
# - removed missing magdecl warning on -b
|
18
|
33 |
# Mar 4, 2014: - added support for missing PITCH/ROLL/HEADING
|
0
|
34 |
|
|
35 |
# General Notes:
|
|
36 |
# - everything (e.g. beams) is numbered from 1
|
|
37 |
# - no support for BT data
|
|
38 |
|
|
39 |
# Soundspeed Correction:
|
|
40 |
# - applied as described in the RDI coord-trans manual
|
|
41 |
# - sound-speed variation over range is ignored (valid for small gradients)
|
|
42 |
# => - same simple correction for all velocity components
|
|
43 |
# - simple correction for cell depths
|
|
44 |
|
|
45 |
# Min %-good (min_pcg):
|
|
46 |
# - nan for records w/o valid velocities
|
|
47 |
# - min(%-good) of the beams used for the velocity solution
|
|
48 |
# - min_pcg does not have to decrease monotonically with distance,
|
|
49 |
# at least when 3-beam solutions are allowed and when -p is used to
|
|
50 |
# edit the data
|
|
51 |
# - non-monotonic min_pcg is particularly obvious with the DYNAMUCK BM_ADCP
|
|
52 |
# data, where one of the beams performed much worse than the others
|
|
53 |
|
|
54 |
require "getopts.pl";
|
|
55 |
$0 =~ m{(.*/)[^/]+};
|
|
56 |
require "$1RDI_BB_Read.pl";
|
|
57 |
require "$1RDI_Coords.pl";
|
|
58 |
require "$1RDI_Utils.pl";
|
|
59 |
|
|
60 |
die("Usage: $0 [-r)ange <first_ens,last_ens>] " .
|
|
61 |
"[output -f)ile <pat[bin%d.raw]>] " .
|
|
62 |
"[-a)ll ens (not just those with good vels)] " .
|
|
63 |
"[-M)agnetic <declination>] " .
|
|
64 |
"[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
|
|
65 |
"[require -4)-beam solutions] [-d)iscard <beam#>] " .
|
|
66 |
"[-%)good <min>] " .
|
|
67 |
"[output -b)eam coordinates] " .
|
|
68 |
"<RDI file>\n")
|
|
69 |
unless (&Getopts("4abd:f:M:p:r:S:") && @ARGV == 1);
|
|
70 |
|
|
71 |
die("$0: -4 and -d are mutually exclusive\n")
|
|
72 |
if ($opt_4 && defined($opt_d));
|
|
73 |
|
|
74 |
die("$0: -p and -b are mutually exclusive\n")
|
|
75 |
if ($opt_b && defined($opt_p));
|
|
76 |
|
|
77 |
$opt_p = 0 unless defined($opt_p);
|
|
78 |
|
|
79 |
$RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions
|
|
80 |
|
|
81 |
print(STDERR "WARNING: magnetic declination not set!\n")
|
|
82 |
unless defined($opt_M) || defined($opt_b);
|
|
83 |
|
|
84 |
$opt_f = 'bin%d.raw' unless defined($opt_f);
|
|
85 |
$ifn = $ARGV[0];
|
|
86 |
|
|
87 |
($first_ens,$last_ens) = split(',',$opt_r)
|
|
88 |
if defined($opt_r);
|
|
89 |
|
|
90 |
($salin,$temp,$depth) = split(',',$opt_S)
|
|
91 |
if defined($opt_S);
|
|
92 |
$variable_ssCorr = ($salin eq '*' || $temp eq '*' || $depth eq '*');
|
|
93 |
|
|
94 |
#----------------------------------------------------------------------
|
|
95 |
|
|
96 |
sub ssCorr($$$$) # sound velocity correction factor
|
|
97 |
{
|
|
98 |
my($e,$S,$T,$D) = @_;
|
|
99 |
$S = $dta{ENSEMBLE}[$e]->{SALINITY} if ($S eq '*');
|
|
100 |
$T = $dta{ENSEMBLE}[$e]->{TEMPERATURE} if ($T eq '*');
|
|
101 |
if ($D eq '*') {
|
|
102 |
print(STDERR "WARNING: soundspeed correction using instrument pressure instead of depth!\n")
|
|
103 |
unless ($warned);
|
|
104 |
$warned = 1;
|
|
105 |
$D = $dta{ENSEMBLE}[$e]->{PRESSURE};
|
|
106 |
}
|
|
107 |
return soundSpeed($S,$T,$D) / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
|
|
108 |
}
|
|
109 |
|
|
110 |
sub min(@) # return minimum
|
|
111 |
{
|
|
112 |
my($min) = 99e99;
|
|
113 |
for (my($i)=0; $i<=$#_; $i++) {
|
|
114 |
$min = $_[$i] if defined($_[$i]) && ($_[$i] < $min);
|
|
115 |
}
|
|
116 |
return ($min == 99e99) ? nan : $min;
|
|
117 |
}
|
|
118 |
|
|
119 |
sub dumpBin($$$) # write time series of single bin
|
|
120 |
{
|
|
121 |
my($b,$fe,$le) = @_;
|
|
122 |
my($file) = sprintf($opt_f,$b+1);
|
|
123 |
|
|
124 |
open(P,">$file") || die("$file: $!\n");
|
|
125 |
print(P "#ANTS#PARAMS# ");
|
|
126 |
foreach my $k (keys(%P)) {
|
|
127 |
print(P "$k\{$P{$k}\} ");
|
|
128 |
}
|
|
129 |
my($pct3b) = ($good_vels[$b] > 0) ? 100*$three_beam[$b]/$good_vels[$b] : nan;
|
|
130 |
printf(STDERR "%.0f%%/%.0f%% ",100*$good_vels[$b]/($le-$fe+1),$pct3b);
|
|
131 |
printf(P "pct_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1));
|
|
132 |
printf(P "pct_3_beam{%.0f} ",$pct3b);
|
|
133 |
printf(P "bin{%d}",$b+1);
|
|
134 |
printf(P " soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!');
|
|
135 |
printf(P " dz{%g}",$dz[$b] *
|
|
136 |
(defined($opt_S) ? ssCorr($fe,$salin,$temp,$depth) : 1)
|
|
137 |
) unless ($variable_ssCorr);
|
|
138 |
print( P "\n");
|
|
139 |
|
|
140 |
print(P "#ANTS#FIELDS# " .
|
|
141 |
"{ensemble} {date} {time} {elapsed} " .
|
|
142 |
"{heading} {pitch} {roll} " .
|
|
143 |
"{sig_heading} {sig_pitch} {sig_roll} " .
|
|
144 |
"{xmit_current} {xmit_voltage} " .
|
|
145 |
"{temp} " .
|
|
146 |
($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") .
|
|
147 |
"{corr1} {corr2} {corr3} {corr4} " .
|
|
148 |
"{amp1} {amp2} {amp3} {amp4} " .
|
|
149 |
($opt_b ? "{pcg1} {pcg2} {pcg3} {pcg4}" : "{3_beam} {min_pcg}")
|
|
150 |
);
|
|
151 |
print(P " {dz}") if ($variable_ssCorr);
|
|
152 |
print(P "\n");
|
|
153 |
|
|
154 |
my($t0) = $dta{ENSEMBLE}[$fe]->{UNIX_TIME};
|
|
155 |
for (my($e)=$fe; $e<=$le; $e++) {
|
|
156 |
next unless ($opt_a || $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]);
|
|
157 |
|
|
158 |
my($ssCorr) = defined($opt_S) ? ssCorr($e,$salin,$temp,$depth) : 1;
|
|
159 |
|
|
160 |
print(P "$dta{ENSEMBLE}[$e]->{NUMBER} ");
|
|
161 |
print(P "$dta{ENSEMBLE}[$e]->{DATE} ");
|
|
162 |
print(P "$dta{ENSEMBLE}[$e]->{TIME} ");
|
|
163 |
printf(P "%d ",$dta{ENSEMBLE}[$e]->{UNIX_TIME}-$t0);
|
18
|
164 |
print(P defined($dta{ENSEMBLE}[$e]->{HEADING}) ? "$dta{ENSEMBLE}[$e]->{HEADING} " : 'nan ');
|
|
165 |
print(P defined($dta{ENSEMBLE}[$e]->{PITCH}) ? "$dta{ENSEMBLE}[$e]->{PITCH} " : 'nan ');
|
|
166 |
print(P defined($dta{ENSEMBLE}[$e]->{ROLL}) ? "$dta{ENSEMBLE}[$e]->{ROLL} " : 'nan ');
|
0
|
167 |
print(P "$dta{ENSEMBLE}[$e]->{HEADING_STDDEV} ");
|
|
168 |
print(P "$dta{ENSEMBLE}[$e]->{PITCH_STDDEV} ");
|
|
169 |
print(P "$dta{ENSEMBLE}[$e]->{ROLL_STDDEV} ");
|
|
170 |
print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT} ");
|
|
171 |
print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} ");
|
|
172 |
print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} ");
|
|
173 |
if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {
|
|
174 |
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] * $ssCorr);
|
|
175 |
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] * $ssCorr);
|
|
176 |
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] * $ssCorr);
|
|
177 |
printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr);
|
|
178 |
} else {
|
|
179 |
print(P "nan nan nan nan ");
|
|
180 |
}
|
|
181 |
print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");
|
|
182 |
print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");
|
|
183 |
if ($opt_b) {
|
|
184 |
print(P "@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]} ");
|
|
185 |
} else {
|
|
186 |
printf(P "%d ",$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]);
|
|
187 |
printf(P "%s ",min(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}));
|
|
188 |
}
|
|
189 |
printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr);
|
|
190 |
print(P "\n");
|
|
191 |
}
|
|
192 |
close(P);
|
|
193 |
}
|
|
194 |
|
|
195 |
#----------------------------------------------------------------------
|
|
196 |
# MAIN
|
|
197 |
#----------------------------------------------------------------------
|
|
198 |
|
|
199 |
$P{RDI_file} = $ifn;
|
|
200 |
$P{mag_decl} = $opt_M if defined($opt_M);
|
|
201 |
|
|
202 |
readData($ifn,\%dta);
|
|
203 |
printf("%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}));
|
|
204 |
$dta{HEADING_BIAS} = -$opt_M; # magnetic declination
|
|
205 |
|
|
206 |
if ($dta{BEAM_COORDINATES}) { # coords
|
|
207 |
$beamCoords = 1;
|
|
208 |
} else {
|
|
209 |
die("$0: -b requires input in beam coordinates\n")
|
|
210 |
if ($opt_b);
|
|
211 |
die("$ifn: only beam and earth coordinates implemented so far\n")
|
|
212 |
if (!$dta{EARTH_COORDINATES});
|
|
213 |
}
|
|
214 |
|
|
215 |
$P{N_ensembles} = @{$dta{ENSEMBLE}};
|
|
216 |
|
|
217 |
for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz
|
|
218 |
$dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};
|
|
219 |
}
|
|
220 |
|
|
221 |
$lastGoodBin = 0;
|
|
222 |
for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities
|
|
223 |
next if (defined($first_ens) &&
|
|
224 |
$dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
|
|
225 |
$P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e
|
|
226 |
unless defined($P{first_ens});
|
|
227 |
last if (defined($last_ens) &&
|
|
228 |
$dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);
|
|
229 |
$P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};
|
|
230 |
$le = $e;
|
|
231 |
|
|
232 |
die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
|
|
233 |
if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);
|
|
234 |
die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")
|
|
235 |
if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});
|
|
236 |
die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")
|
|
237 |
if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});
|
|
238 |
|
|
239 |
for (my($b)=0; $b<$dta{N_BINS}; $b++) {
|
|
240 |
if (defined($opt_d)) {
|
|
241 |
undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
|
|
242 |
undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
|
|
243 |
}
|
|
244 |
for (my($i)=0; $i<4; $i++) {
|
|
245 |
if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) {
|
|
246 |
undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]);
|
|
247 |
undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]);
|
|
248 |
}
|
|
249 |
}
|
|
250 |
@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords
|
|
251 |
? velInstrumentToEarth(\%dta,$e,
|
|
252 |
velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
|
|
253 |
)
|
|
254 |
: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
|
|
255 |
unless ($opt_b);
|
|
256 |
unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {
|
|
257 |
undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
|
|
258 |
next;
|
|
259 |
}
|
|
260 |
$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag;
|
|
261 |
$three_beam[$b] += $RDI_Coords::threeBeamFlag;
|
|
262 |
$dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1;
|
|
263 |
$good_vels[$b]++;
|
|
264 |
$lastGoodBin = $b if ($b > $lastGoodBin);
|
|
265 |
$firstGoodEns = $e unless defined($firstGoodEns);
|
|
266 |
$lastGoodEns = $e;
|
|
267 |
}
|
|
268 |
}
|
|
269 |
|
|
270 |
unless (defined($opt_r)) {
|
|
271 |
$fe = $firstGoodEns;
|
|
272 |
$le = $lastGoodEns;
|
|
273 |
}
|
|
274 |
|
|
275 |
$firstBin = 0;
|
|
276 |
$lastBin = $lastGoodBin;
|
|
277 |
|
|
278 |
print( STDERR "Start : $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n");
|
|
279 |
print( STDERR "End : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n");
|
|
280 |
printf(STDERR "Bins : %d-%d\n",$firstBin+1,$lastBin+1);
|
|
281 |
printf(STDERR "3-Beam : %d %d %d %d\n",$RDI_Coords::threeBeam_1,
|
|
282 |
$RDI_Coords::threeBeam_2,
|
|
283 |
$RDI_Coords::threeBeam_3,
|
|
284 |
$RDI_Coords::threeBeam_4);
|
|
285 |
|
|
286 |
print(STDERR "Good/3-beam: ");
|
|
287 |
for ($b=$firstBin; $b<=$lastBin; $b++) { # generate output
|
|
288 |
dumpBin($b,$fe,$le);
|
|
289 |
}
|
|
290 |
print(STDERR "\n");
|
|
291 |
|
|
292 |
exit(0);
|