|
1 #!/usr/bin/perl |
|
2 #====================================================================== |
|
3 # M E A N P R O F |
|
4 # doc: Fri Feb 22 08:40:18 2008 |
|
5 # dlm: Fri Feb 22 18:02:40 2008 |
|
6 # (c) 2008 A.M. Thurnherr |
|
7 # uE-Info: 232 0 NIL 0 0 72 74 2 4 NIL ofnI |
|
8 #====================================================================== |
|
9 |
|
10 # extract time-averaged mean profile from ADCP data |
|
11 |
|
12 # HISTORY: |
|
13 # Feb 22, 2008: - created from [listBins] |
|
14 |
|
15 # Soundspeed Correction: |
|
16 # - applied as described in the RDI coord-trans manual |
|
17 # - sound-speed variation over range is ignored (valid for small gradients) |
|
18 # => - same simple correction for all velocity components |
|
19 # - simple correction for cell depths |
|
20 |
|
21 require "getopts.pl"; |
|
22 $0 =~ m{(.*/)[^/]+}; |
|
23 require "$1RDI_BB_Read.pl"; |
|
24 require "$1RDI_Coords.pl"; |
|
25 require "$1RDI_Utils.pl"; |
|
26 |
|
27 die("Usage: $0 [-r)ange <first_ens,last_ens>] " . |
|
28 "[-Q)uiet (stats only)] " . |
|
29 "[-S)oundspeed correction <salin|*,temp|*,depth|*> " . |
|
30 "[require -4)-beam solutions] [-d)iscard <beam#>] " . |
|
31 "[-%)good <min>] " . |
|
32 "[output -b)eam coordinates] " . |
|
33 "[-M)agnetic <declination>] " . |
|
34 "[-D)epth <depth>] " . |
|
35 "<RDI file>\n") |
|
36 unless (&Getopts("4bd:D:M:p:r:QS:") && @ARGV == 1); |
|
37 |
|
38 die("$0: -4 and -d are mutually exclusive\n") |
|
39 if ($opt_4 && defined($opt_d)); |
|
40 |
|
41 die("$0: -p and -b are mutually exclusive\n") |
|
42 if ($opt_b && defined($opt_p)); |
|
43 |
|
44 $opt_p = 0 unless defined($opt_p); |
|
45 |
|
46 $RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions |
|
47 |
|
48 print(STDERR "WARNING: magnetic declination not set!\n") |
|
49 unless defined($opt_M) || defined($opt_b); |
|
50 |
|
51 $ifn = $ARGV[0]; |
|
52 |
|
53 ($first_ens,$last_ens) = split(',',$opt_r) |
|
54 if defined($opt_r); |
|
55 |
|
56 ($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S) |
|
57 if defined($opt_S); |
|
58 die("$0: Cannot do variable soundspeed correction (implementation restriction)\n") |
|
59 if ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*'); |
|
60 |
|
61 #---------------------------------------------------------------------- |
|
62 # Read & Check Data, Transform Velocities |
|
63 #---------------------------------------------------------------------- |
|
64 |
|
65 $P{RDI_file} = $ifn; |
|
66 $P{mag_decl} = $opt_M if defined($opt_M); |
|
67 |
|
68 print(STDERR "reading $ifn: "); |
|
69 readData($ifn,\%dta); |
|
70 printf(STDERR "%d complete ensembles.\n",scalar(@{$dta{ENSEMBLE}})); |
|
71 $dta{HEADING_BIAS} = -$opt_M; # magnetic declination |
|
72 |
|
73 if ($dta{BEAM_COORDINATES}) { # coords |
|
74 $beamCoords = 1; |
|
75 } else { |
|
76 die("$0: -b requires input in beam coordinates\n") |
|
77 if ($opt_b); |
|
78 die("$ifn: only beam and earth coordinates implemented so far\n") |
|
79 if (!$dta{EARTH_COORDINATES}); |
|
80 } |
|
81 |
|
82 for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz |
|
83 $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH}; |
|
84 } |
|
85 |
|
86 $lastGoodBin = 0; |
|
87 for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities |
|
88 next if (defined($first_ens) && |
|
89 $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens); |
|
90 $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e |
|
91 unless defined($P{first_ens}); |
|
92 last if (defined($last_ens) && |
|
93 $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens); |
|
94 $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER}; |
|
95 $le = $e; |
|
96 |
|
97 die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") |
|
98 if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4); |
|
99 die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n") |
|
100 if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR}); |
|
101 die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") |
|
102 if ($dta{ENSEMBLE}[$e]->{LOW_GAIN}); |
|
103 |
|
104 for (my($b)=0; $b<$dta{N_BINS}; $b++) { |
|
105 if (defined($opt_d)) { |
|
106 undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]); |
|
107 undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]); |
|
108 } |
|
109 for (my($i)=0; $i<4; $i++) { |
|
110 if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) { |
|
111 undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]); |
|
112 undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]); |
|
113 } |
|
114 } |
|
115 @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords |
|
116 ? velInstrumentToEarth(\%dta,$e, |
|
117 velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
|
118 ) |
|
119 : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) |
|
120 unless ($opt_b); |
|
121 |
|
122 $sum_corr1[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0]; |
|
123 $sum_corr2[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][1]; |
|
124 $sum_corr3[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][2]; |
|
125 $sum_corr4[$b] += $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][3]; |
|
126 |
|
127 $sum_amp1[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][0]; |
|
128 $sum_amp2[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][1]; |
|
129 $sum_amp3[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][2]; |
|
130 $sum_amp4[$b] += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3]; |
|
131 |
|
132 unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) { |
|
133 undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}); |
|
134 next; |
|
135 } |
|
136 |
|
137 $dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag; |
|
138 $three_beam[$b] += $RDI_Coords::threeBeamFlag; |
|
139 $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1; |
|
140 $good_vels[$b]++; |
|
141 $lastGoodBin = $b if ($b > $lastGoodBin); |
|
142 $firstGoodEns = $e unless defined($firstGoodEns); |
|
143 $lastGoodEns = $e; |
|
144 |
|
145 $sum_u[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]; |
|
146 $sum_v[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]; |
|
147 $sum_w[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]; |
|
148 $sum_e[$b] += $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] |
|
149 unless ($RDI_Coords::threeBeamFlag); |
|
150 |
|
151 $sum_pcg1[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]; |
|
152 $n_pcg1[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]); |
|
153 $sum_pcg2[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]; |
|
154 $n_pcg2[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]); |
|
155 $sum_pcg3[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]; |
|
156 $n_pcg3[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]); |
|
157 $sum_pcg4[$b] += $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]; |
|
158 $n_pcg4[$b]++ if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]); |
|
159 } |
|
160 } |
|
161 |
|
162 unless (defined($opt_r)) { |
|
163 $fe = $firstGoodEns; |
|
164 $le = $lastGoodEns; |
|
165 } |
|
166 $nEns = $le - $fe + 1; |
|
167 die("$0: insufficient data\n") if ($nEns < 2); |
|
168 $P{N_ensembles} = $nEns; |
|
169 |
|
170 $firstBin = 0; |
|
171 $lastBin = $lastGoodBin; |
|
172 |
|
173 print( STDERR "Start : $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n"); |
|
174 print( STDERR "End : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n"); |
|
175 printf(STDERR "Bins : %d-%d\n",$firstBin+1,$lastBin+1); |
|
176 printf(STDERR "3-Beam : %d %d %d %d\n",$RDI_Coords::threeBeam_1, |
|
177 $RDI_Coords::threeBeam_2, |
|
178 $RDI_Coords::threeBeam_3, |
|
179 $RDI_Coords::threeBeam_4) |
|
180 unless ($opt_b); |
|
181 |
|
182 #---------------------------------------------------------------------- |
|
183 # Calculate Stddevs |
|
184 #---------------------------------------------------------------------- |
|
185 |
|
186 for ($b=0; $b<=$lastGoodBin; $b++) { |
|
187 $mean_corr1[$b] = $sum_corr1[$b] / $nEns; $mean_corr2[$b] = $sum_corr2[$b] / $nEns; |
|
188 $mean_corr3[$b] = $sum_corr3[$b] / $nEns; $mean_corr4[$b] = $sum_corr4[$b] / $nEns; |
|
189 $mean_amp1[$b] = $sum_amp1[$b] / $nEns; $mean_amp2[$b] = $sum_amp2[$b] / $nEns; |
|
190 $mean_amp3[$b] = $sum_amp3[$b] / $nEns; $mean_amp4[$b] = $sum_amp4[$b] / $nEns; |
|
191 $mean_pcg1[$b] = $sum_pcg1[$b] / $n_pcg1[$b]; $mean_pcg2[$b] = $sum_pcg2[$b] / $n_pcg2[$b]; |
|
192 $mean_pcg3[$b] = $sum_pcg3[$b] / $n_pcg3[$b]; $mean_pcg4[$b] = $sum_pcg4[$b] / $n_pcg4[$b]; |
|
193 $mean_u[$b] = $sum_u[$b] / $good_vels[$b]; $mean_v[$b] = $sum_v[$b] / $good_vels[$b]; |
|
194 $mean_w[$b] = $sum_w[$b] / $good_vels[$b]; |
|
195 $mean_e[$b] = $sum_e[$b] / ($good_vels[$b] - $three_beam[$b]); |
|
196 } |
|
197 |
|
198 for ($e=$fe; $e<=$le; $e++) { |
|
199 for ($b=0; $b<=$lastGoodBin; $b++) { |
|
200 $sumsq_corr1[$b] += ($mean_corr1[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][0])**2; |
|
201 $sumsq_corr2[$b] += ($mean_corr2[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][1])**2; |
|
202 $sumsq_corr3[$b] += ($mean_corr3[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][2])**2; |
|
203 $sumsq_corr4[$b] += ($mean_corr4[$b] - $dta{ENSEMBLE}[$e]->{CORRELATION}[$b][3])**2; |
|
204 |
|
205 $sumsq_amp1[$b] += ($mean_amp1[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][0])**2; |
|
206 $sumsq_amp2[$b] += ($mean_amp2[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][1])**2; |
|
207 $sumsq_amp3[$b] += ($mean_amp3[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][2])**2; |
|
208 $sumsq_amp4[$b] += ($mean_amp4[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3])**2; |
|
209 |
|
210 $sumsq_pcg1[$b] += ($mean_pcg1[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0])**2 |
|
211 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]); |
|
212 $sumsq_pcg2[$b] += ($mean_pcg2[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1])**2 |
|
213 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]); |
|
214 $sumsq_pcg3[$b] += ($mean_pcg3[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2])**2 |
|
215 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]); |
|
216 $sumsq_pcg4[$b] += ($mean_pcg4[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3])**2 |
|
217 if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]); |
|
218 |
|
219 next unless ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]); |
|
220 |
|
221 $sumsq_u[$b] += ($mean_u[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])**2; |
|
222 $sumsq_v[$b] += ($mean_v[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1])**2; |
|
223 $sumsq_w[$b] += ($mean_w[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2])**2; |
|
224 |
|
225 next if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]); |
|
226 |
|
227 $sumsq_e[$b] += ($mean_e[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])**2; |
|
228 } |
|
229 } |
|
230 |
|
231 for ($b=0; $b<=$lastGoodBin; $b++) { |
|
232 $var_corr1[$b] = $sumsq_corr1[$b] / ($nEns-1); $var_corr2[$b] = $sumsq_corr2[$b] / ($nEns-1); |
|
233 $var_corr3[$b] = $sumsq_corr3[$b] / ($nEns-1); $var_corr4[$b] = $sumsq_corr4[$b] / ($nEns-1); |
|
234 $var_amp1[$b] = $sumsq_amp1[$b] / ($nEns-1); $var_amp2[$b] = $sumsq_amp2[$b] / ($nEns-1); |
|
235 $var_amp3[$b] = $sumsq_amp3[$b] / ($nEns-1); $var_amp4[$b] = $sumsq_amp4[$b] / ($nEns-1); |
|
236 $var_pcg1[$b] = $sumsq_pcg1[$b] / ($n_pcg1[$b]-1) if ($n_pcg1[$b] > 1); |
|
237 $var_pcg2[$b] = $sumsq_pcg2[$b] / ($n_pcg2[$b]-1) if ($n_pcg2[$b] > 1); |
|
238 $var_pcg3[$b] = $sumsq_pcg3[$b] / ($n_pcg3[$b]-1) if ($n_pcg3[$b] > 1); |
|
239 $var_pcg4[$b] = $sumsq_pcg4[$b] / ($n_pcg4[$b]-1) if ($n_pcg4[$b] > 1); |
|
240 next unless ($good_vels[$b] > 1); |
|
241 $var_u[$b] = $sumsq_u[$b] / ($good_vels[$b]-1); |
|
242 $var_v[$b] = $sumsq_v[$b] / ($good_vels[$b]-1); |
|
243 $var_w[$b] = $sumsq_w[$b] / ($good_vels[$b]-1); |
|
244 next unless ($good_vels[$b] - $three_beam[$b] > 1); |
|
245 $var_e[$b] = $sumsq_e[$b] / ($good_vels[$b] - $three_beam[$b] - 1); |
|
246 } |
|
247 |
|
248 #---------------------------------------------------------------------- |
|
249 # Calculate Beam Statistics |
|
250 #---------------------------------------------------------------------- |
|
251 |
|
252 #---------------------------------------------------------------------- |
|
253 # Produce Output |
|
254 #---------------------------------------------------------------------- |
|
255 |
|
256 unless ($opt_Q) { |
|
257 my($ssCorr) = defined($opt_S) |
|
258 ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) |
|
259 : 1; |
|
260 |
|
261 print("#ANTS#PARAMS# "); |
|
262 foreach my $k (keys(%P)) { |
|
263 print("$k\{$P{$k}\} "); |
|
264 } |
|
265 printf("soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!'); |
|
266 print("\n"); |
|
267 |
|
268 print("#ANTS#FIELDS# " . |
|
269 "{bin} {dz} " . |
|
270 (defined($opt_D) ? "{depth} " : "") . |
|
271 ($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") . |
|
272 ($opt_b ? "{sig_v1} {sig_v2} {sig_v3} {sig_v4} " : "{sig_u} {sig_v} {sig_w} {sig_err_vel} ") . |
|
273 "{corr1} {corr2} {corr3} {corr4} " . |
|
274 "{sig_corr1} {sig_corr2} {sig_corr3} {sig_corr4} " . |
|
275 "{amp1} {amp2} {amp3} {amp4} " . |
|
276 "{sig_amp1} {sig_amp2} {sig_amp3} {sig_amp4} " . |
|
277 "{pcg1} {pcg2} {pcg3} {pcg4} " . |
|
278 "{sig_pcg1} {sig_pcg2} {sig_pcg3} {sig_pcg4}" . |
|
279 "\n" |
|
280 ); |
|
281 |
|
282 for ($b=$firstBin; $b<=$lastBin; $b++) { |
|
283 printf("%d %.1f ",$b+1,$dz[$b]*$ssCorr); |
|
284 printf("%.1f ",$opt_D - $dz[$b]*$ssCorr) |
|
285 if defined($opt_D); |
|
286 |
|
287 printf("%s ",defined($mean_u[$b]) ? $mean_u[$b] : nan); |
|
288 printf("%s ",defined($mean_v[$b]) ? $mean_v[$b] : nan); |
|
289 printf("%s ",defined($mean_w[$b]) ? $mean_w[$b] : nan); |
|
290 printf("%s ",defined($mean_e[$b]) ? $mean_e[$b] : nan); |
|
291 |
|
292 printf("%s ",defined($var_u[$b]) ? sqrt($var_u[$b]) : nan); |
|
293 printf("%s ",defined($var_v[$b]) ? sqrt($var_v[$b]) : nan); |
|
294 printf("%s ",defined($var_w[$b]) ? sqrt($var_w[$b]) : nan); |
|
295 printf("%s ",defined($var_e[$b]) ? sqrt($var_e[$b]) : nan); |
|
296 |
|
297 printf("%g %g %g %g ",$mean_corr1[$b],$mean_corr2[$b], |
|
298 $mean_corr3[$b],$mean_corr4[$b]); |
|
299 printf("%g %g %g %g ",sqrt($var_corr1[$b]),sqrt($var_corr2[$b]), |
|
300 sqrt($var_corr3[$b]),sqrt($var_corr4[$b])); |
|
301 |
|
302 printf("%g %g %g %g ",$mean_amp1[$b],$mean_amp2[$b], |
|
303 $mean_amp3[$b],$mean_amp4[$b]); |
|
304 printf("%g %g %g %g ",sqrt($var_amp1[$b]),sqrt($var_amp2[$b]), |
|
305 sqrt($var_amp3[$b]),sqrt($var_amp4[$b])); |
|
306 |
|
307 if ($good_vels[$b] > 0) { |
|
308 printf("%g %g %g %g ",$mean_pcg1[$b],$mean_pcg2[$b], |
|
309 $mean_pcg3[$b],$mean_pcg4[$b]); |
|
310 printf("%g %g %g %g\n",sqrt($var_pcg1[$b]),sqrt($var_pcg2[$b]), |
|
311 sqrt($var_pcg3[$b]),sqrt($var_pcg4[$b])); |
|
312 } else { |
|
313 print("nan nan nan nan "); |
|
314 print("nan nan nan nan\n"); |
|
315 } |
|
316 } |
|
317 } |
|
318 |
|
319 exit(0); |