1 #!/usr/bin/perl |
1 #!/usr/bin/perl |
2 #====================================================================== |
2 #====================================================================== |
3 # L A D C P _ W S P E C |
3 # L A D C P _ W S P E C |
4 # doc: Thu Jun 11 12:02:49 2015 |
4 # doc: Thu Jun 11 12:02:49 2015 |
5 # dlm: Sun Mar 12 12:01:59 2017 |
5 # dlm: Thu Nov 1 10:09:48 2018 |
6 # (c) 2012 A.M. Thurnherr |
6 # (c) 2012 A.M. Thurnherr |
7 # uE-Info: 48 17 NIL 0 0 72 10 2 4 NIL ofnI |
7 # uE-Info: 39 29 NIL 0 0 72 10 2 4 NIL ofnI |
8 #====================================================================== |
8 #====================================================================== |
9 |
9 |
10 $antsSummary = 'calculate VKE window spectra from LADCP profiles'; |
10 $antsSummary = 'calculate VKE window spectra from LADCP profiles'; |
11 |
11 |
12 # HISTORY: |
12 # HISTORY: |
30 # - added w.nsamp.avg |
30 # - added w.nsamp.avg |
31 # Mar 31, 2016: - changed version %PARAM |
31 # Mar 31, 2016: - changed version %PARAM |
32 # - BUG: nspec was nan insted of 0 |
32 # - BUG: nspec was nan insted of 0 |
33 # - replaced wspec:: %PARAM-prefix with LADCP_wspec:: |
33 # - replaced wspec:: %PARAM-prefix with LADCP_wspec:: |
34 # Mar 12, 2017: - removed ANTSBIN (which is not public) |
34 # Mar 12, 2017: - removed ANTSBIN (which is not public) |
|
35 # Dec 9, 2017: - added $antsSuppressCommonOptions = 1; |
|
36 # Dec 14, 2017: - added w.rms to output |
|
37 # May 13, 2018: - BUG: Removal of higher order polynomials (-o > 0) did not work |
|
38 # May 16, 2018: - modified depth.{min,max} to respect input resolution |
|
39 # Nov 1, 2018: - cosmetics |
35 |
40 |
36 ($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
41 ($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
37 ($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
42 ($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
38 $WCALC = '.' if ($WCALC eq ''); |
43 $WCALC = '.' if ($WCALC eq ''); |
39 |
44 |
50 |
55 |
51 #---------------------------------------------------------------------- |
56 #---------------------------------------------------------------------- |
52 # Usage |
57 # Usage |
53 #---------------------------------------------------------------------- |
58 #---------------------------------------------------------------------- |
54 |
59 |
|
60 $antsSuppressCommonOptions = 1; |
55 &antsUsage('bc:dg:o:s:tuw:',0, |
61 &antsUsage('bc:dg:o:s:tuw:',0, |
56 '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]', |
62 '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]', |
57 '[-d)own/-u)pcast-only] [exclude -b)ottom window]', |
63 '[-d)own/-u)pcast-only] [exclude -b)ottom window]', |
58 '[shortwave -c)utoff <kz or lambda>]', |
64 '[shortwave -c)utoff <kz or lambda>]', |
59 '[-s)urface <layer depth to exclude[150m]>', |
65 '[-s)urface <layer depth to exclude[150m]>', |
60 '[-g)ap <max depth layer to fill with interpolation[40m]>]', |
66 '[-g)ap <max depth layer to fill with interpolation[40m]>]', |
61 '[-w)indow <power-of-two input-records[32]>]', |
67 '[-w)indow <power-of-two input-records>]', |
62 '[LADCP-profile(s)]'); |
68 '[LADCP-profile(s)]'); |
63 |
69 |
64 &antsIntOpt(\$opt_o,0); # polynomial order to remove |
70 &antsIntOpt(\$opt_o,0); # polynomial order to remove |
65 if ($opt_o >= 0) { # init model |
71 if ($opt_o >= 0) { # init model |
66 &modelUsage(); |
72 &modelUsage(); |
130 $dz = &antsXCheck($zfnr,0,$#ants_,1.01); # calc dT; 1% jitter |
136 $dz = &antsXCheck($zfnr,0,$#ants_,1.01); # calc dT; 1% jitter |
131 &antsAddParams("LADCP_wspec::input_depth_resolution",$dz); |
137 &antsAddParams("LADCP_wspec::input_depth_resolution",$dz); |
132 |
138 |
133 $opt_g = int(($opt_g - 1) / $dz); # [m] -> [records] |
139 $opt_g = int(($opt_g - 1) / $dz); # [m] -> [records] |
134 |
140 |
135 unless (defined($opt_w)) { # window size |
141 unless (defined($opt_w)) { # default window size: largest pwr-of-two <= 600m |
136 for ($opt_w=32; $opt_w*$dz>600; $opt_w/=2) {} |
142 for ($opt_w=32; $opt_w*$dz>600; $opt_w/=2) {} |
137 &antsInfo("%d-m windows ($opt_w records)",$opt_w*$dz); |
143 &antsInfo("%d-m windows ($opt_w samples)",$opt_w*$dz); |
138 } |
144 } |
139 &antsAddParams('LADCP_wspec::window_size',$opt_w,'LADCP_wspec::output_depth_resolution',$dz*$opt_w); |
145 &antsAddParams('LADCP_wspec::window_size',$opt_w,'LADCP_wspec::output_depth_resolution',$dz*$opt_w); |
140 |
146 |
141 croak(sprintf("$0: insufficient data (%d records found, %d required)\n",scalar(@ants_),$opt_w)) |
147 croak(sprintf("$0: insufficient data (%d records found, %d required)\n",scalar(@ants_),$opt_w)) |
142 unless (@ants_ >= $opt_w); |
148 unless (@ants_ >= $opt_w); |
144 $zrange = $opt_w * $dz; # NB: not equal to max-min!!! |
150 $zrange = $opt_w * $dz; # NB: not equal to max-min!!! |
145 $resolution_bandwidth = 1 / $zrange; |
151 $resolution_bandwidth = 1 / $zrange; |
146 $resolution_bandwidth *= 2*$PI; |
152 $resolution_bandwidth *= 2*$PI; |
147 &antsAddParams('LADCP_wspec::resolution_bandwidth',$resolution_bandwidth); |
153 &antsAddParams('LADCP_wspec::resolution_bandwidth',$resolution_bandwidth); |
148 |
154 |
149 push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','nspec','w.nsamp.avg'); |
155 push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','w.rms','nspec','w.nsamp.avg'); |
150 for (my($i)=0; $i<$opt_w/2+1; $i++) { |
156 for (my($i)=0; $i<$opt_w/2+1; $i++) { |
151 my($kz) = 2*$PI*$i/$zrange; |
157 my($kz) = 2*$PI*$i/$zrange; |
152 last if ($kz > $kzlim); |
158 last if ($kz > $kzlim); |
153 &antsAddParams(sprintf('LADCP_wspec::k.%d',$i),$kz); |
159 &antsAddParams(sprintf('LADCP_wspec::k.%d',$i),$kz); |
154 &antsAddParams(sprintf('LADCP_wspec::lambda.%d',$i),$i ? $zrange/$i : inf); |
160 &antsAddParams(sprintf('LADCP_wspec::lambda.%d',$i),$i ? $zrange/$i : inf); |
231 $out[0] = -1; |
237 $out[0] = -1; |
232 $fromR = @ants_ - $opt_w; |
238 $fromR = @ants_ - $opt_w; |
233 $opt_b = 1; |
239 $opt_b = 1; |
234 } |
240 } |
235 |
241 |
|
242 #-------------------------------------------------- |
|
243 # calculate rms w in window |
|
244 # - also determines if there are missing y values |
|
245 #-------------------------------------------------- |
|
246 |
|
247 my($dc_gap) = $opt_u; # exclude dc with -d, uc with -u |
|
248 my($uc_gap) = $opt_d; |
|
249 my($sumsq,$n) = (0,0); |
|
250 for (my($r)=$fromR; $r<$fromR+$opt_w; $r++) { |
|
251 if (numberp($ants_[$r][$dcwfnr])) { |
|
252 $sumsq += $ants_[$r][$dcwfnr]**2; |
|
253 $n++; |
|
254 } else { |
|
255 $dc_gap = 1; |
|
256 } |
|
257 if (numberp($ants_[$r][$ucwfnr])) { |
|
258 $sumsq += $ants_[$r][$ucwfnr]**2; |
|
259 $n++; |
|
260 } else { |
|
261 $uc_gap = 1; |
|
262 } |
|
263 } |
|
264 my($wrms) = ($n > 0) ? sqrt($sumsq/$n) : nan; |
|
265 |
236 #----------------------------------- |
266 #----------------------------------- |
237 # output nan on non-numeric y values |
267 # output nan on non-numeric y values |
238 #----------------------------------- |
268 #----------------------------------- |
239 |
|
240 my($dc_gap) = $opt_u; # exclude dc with -d, uc with -u |
|
241 my($uc_gap) = $opt_d; |
|
242 for (my($r)=$fromR; $r<$fromR+$opt_w; $r++) { |
|
243 $dc_gap |= !numberp($ants_[$r][$dcwfnr]); |
|
244 $uc_gap |= !numberp($ants_[$r][$ucwfnr]); |
|
245 } |
|
246 |
269 |
247 if ($dc_gap && $uc_gap) { |
270 if ($dc_gap && $uc_gap) { |
248 push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); |
271 push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); |
249 if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) { |
272 if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) { |
250 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); |
273 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); |
253 push(@out,$ants_[$fromR][$zfnr]); |
276 push(@out,$ants_[$fromR][$zfnr]); |
254 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); |
277 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); |
255 } |
278 } |
256 push(@out,defined($habfnr) ? # hab |
279 push(@out,defined($habfnr) ? # hab |
257 avgF($habfnr,$fromR) : nan); |
280 avgF($habfnr,$fromR) : nan); |
|
281 push(@out,$wrms); # rms w |
258 push(@out,0); # nspec |
282 push(@out,0); # nspec |
259 push(@out,nan); # w.nsamp.avg |
283 push(@out,nan); # w.nsamp.avg |
260 for ($i=0; $i<=$opt_w/2; $i++) { # power |
284 for ($i=0; $i<=$opt_w/2; $i++) { # power |
261 push(@out,nan); |
285 push(@out,nan); |
262 } |
286 } |
283 &modelInit(); # dc |
307 &modelInit(); # dc |
284 for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; } |
308 for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; } |
285 &lfit($zfnr,$dcwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1); |
309 &lfit($zfnr,$dcwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1); |
286 &modelCleanup(); |
310 &modelCleanup(); |
287 for (my($i)=0; $i<$opt_w; $i++) { |
311 for (my($i)=0; $i<$opt_w; $i++) { |
288 modelEvaluate($i,$zfnr,\@afunc); |
312 modelEvaluate($fromR+$i,$zfnr,\@afunc); |
289 for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) { |
313 for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) { |
290 $calc += $A[$p] * $afunc[$p]; |
314 $calc += $A[$p] * $afunc[$p]; |
291 } |
315 } |
292 $ants_[$fromR+$i][$dcwfnr] -= $calc; |
316 $ants_[$fromR+$i][$dcwfnr] -= $calc; |
293 } |
317 } |
296 &modelInit(); # uc |
320 &modelInit(); # uc |
297 for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; } |
321 for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; } |
298 &lfit($zfnr,$ucwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1); |
322 &lfit($zfnr,$ucwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1); |
299 &modelCleanup(); |
323 &modelCleanup(); |
300 for (my($i)=0; $i<$opt_w; $i++) { |
324 for (my($i)=0; $i<$opt_w; $i++) { |
301 modelEvaluate($i,$zfnr,\@afunc); |
325 modelEvaluate($fromR+$i,$zfnr,\@afunc); |
302 for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) { |
326 for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) { |
303 $calc += $A[$p] * $afunc[$p]; |
327 $calc += $A[$p] * $afunc[$p]; |
304 } |
328 } |
305 $ants_[$fromR+$i][$ucwfnr] -= $calc; |
329 $ants_[$fromR+$i][$ucwfnr] -= $calc; |
306 } |
330 } |
339 @dc_pwr = &pgram_onesided($opt_w,@dc_coeff) # total power |
363 @dc_pwr = &pgram_onesided($opt_w,@dc_coeff) # total power |
340 unless ($dc_gap); |
364 unless ($dc_gap); |
341 @uc_pwr = &pgram_onesided($opt_w,@uc_coeff) |
365 @uc_pwr = &pgram_onesided($opt_w,@uc_coeff) |
342 unless ($uc_gap); |
366 unless ($uc_gap); |
343 |
367 |
344 push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); # middle z |
368 # push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); # middle z |
|
369 push(@out,avgF($zfnr,$fromR)); # average z |
345 if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) { # input descending |
370 if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) { # input descending |
346 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); # z.min |
371 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]-$dz/2); # z.min |
347 push(@out,$ants_[$fromR][$zfnr]); # z.max |
372 push(@out,$ants_[$fromR][$zfnr]+$dz/2); # z.max |
348 } else { # input ascending |
373 } else { # input ascending |
349 push(@out,$ants_[$fromR][$zfnr]); # z.min |
374 push(@out,$ants_[$fromR][$zfnr]-$dz/2); # z.min |
350 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); # z.max |
375 push(@out,$ants_[$fromR+$opt_w-1][$zfnr]+$dz/2); # z.max |
351 } |
376 } |
352 push(@out,defined($habfnr) ? # hab |
377 push(@out,defined($habfnr) ? # hab |
353 avgF($habfnr,$fromR) : nan); |
378 avgF($habfnr,$fromR) : nan); |
|
379 push(@out,$wrms); # w.rms |
354 my($nspec) = !$dc_gap + !$uc_gap; # nspec |
380 my($nspec) = !$dc_gap + !$uc_gap; # nspec |
355 push(@out,$nspec); |
381 push(@out,$nspec); |
356 my($nsamp_sum) = my($nsn) = 0; # w.nsamp.avg |
382 my($nsamp_sum) = my($nsn) = 0; # w.nsamp.avg |
357 $nsamp_sum+=medianF($dcsfnr,$fromR),$nsn++ unless ($dc_gap); # median to avoid biasing by short bottle stops |
383 $nsamp_sum+=medianF($dcsfnr,$fromR),$nsn++ unless ($dc_gap); # median to avoid biasing by short bottle stops |
358 $nsamp_sum+=medianF($ucsfnr,$fromR),$nsn++ unless ($uc_gap); |
384 $nsamp_sum+=medianF($ucsfnr,$fromR),$nsn++ unless ($uc_gap); |