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