diff --git a/meanProf b/meanProf --- a/meanProf +++ b/meanProf @@ -2,9 +2,9 @@ #====================================================================== # M E A N P R O F # doc: Fri Feb 22 08:40:18 2008 -# dlm: Wed Mar 16 07:02:38 2016 +# dlm: Tue Apr 10 09:01:50 2018 # (c) 2008 A.M. Thurnherr -# uE-Info: 14 49 NIL 0 0 72 0 2 4 NIL ofnI +# uE-Info: 19 35 NIL 0 0 72 0 2 4 NIL ofnI #====================================================================== # extract time-averaged mean profile from ADCP data @@ -12,8 +12,14 @@ # HISTORY: # Feb 22, 2008: - created from [listBins] # Mar 16, 2016: - adapted to new Getopt library +# Apr 2, 2018: - BUG: velBeamToInstrument() used old usage +# Apr 9, 2018: - adapted to "new" readData() ensemble limits +# - added -l to set final bin +# - BUG: division by zero in empty bins +# Apr 10, 2018: - activate output # Soundspeed Correction: +# - based on first ensemble only # - applied as described in the RDI coord-trans manual # - sound-speed variation over range is ignored (valid for small gradients) # => - same simple correction for all velocity components @@ -21,12 +27,11 @@ use Getopt::Std; -$0 =~ m{(.*/)[^/]+}; -require "$1RDI_BB_Read.pl"; -require "$1RDI_Coords.pl"; -require "$1RDI_Utils.pl"; +$ADCP_tools_minVersion = 2.2; +($ADCP_TOOLS) = ($0 =~ m{(.*/)[^/]+}); +require "$ADCP_TOOLS/ADCP_tools_lib.pl"; -die("Usage: $0 [-r)ange ] " . +die("Usage: $0 [-r)ange ] [-l)ast ] " . "[-Q)uiet (stats only)] " . "[-S)oundspeed correction " . "[require -4)-beam solutions] [-d)iscard ] " . @@ -35,7 +40,7 @@ "[-M)agnetic ] " . "[-D)epth ] " . "\n") - unless (&getopts("4bd:D:M:p:r:QS:") && @ARGV == 1); + unless (&getopts("4bd:D:l:M:p:r:QS:") && @ARGV == 1); die("$0: -4 and -d are mutually exclusive\n") if ($opt_4 && defined($opt_d)); @@ -52,9 +57,6 @@ $ifn = $ARGV[0]; -($first_ens,$last_ens) = split(',',$opt_r) - if defined($opt_r); - ($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S) if defined($opt_S); die("$0: Cannot do variable soundspeed correction (implementation restriction)\n") @@ -68,7 +70,12 @@ $P{mag_decl} = $opt_M if defined($opt_M); print(STDERR "reading $ifn: "); -readData($ifn,\%dta); +if (defined($opt_r)) { # read selected range + my($fe,$le) = split(',',$opt_r); + readData($ifn,\%dta,$fe,$le,$opt_l); +} else { # read entire file + readData($ifn,\%dta,undef,undef,$opt_l); +} printf(STDERR "%d complete ensembles.\n",scalar(@{$dta{ENSEMBLE}})); $dta{HEADING_BIAS} = -$opt_M; # magnetic declination @@ -87,12 +94,8 @@ $lastGoodBin = 0; for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities - next if (defined($first_ens) && - $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens); $P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e unless defined($P{first_ens}); - last if (defined($last_ens) && - $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens); $P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER}; $le = $e; @@ -116,7 +119,7 @@ } @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords ? velInstrumentToEarth(\%dta,$e, - velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) + velBeamToInstrument(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) ) : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) unless ($opt_b); @@ -190,11 +193,18 @@ $mean_corr3[$b] = $sum_corr3[$b] / $nEns; $mean_corr4[$b] = $sum_corr4[$b] / $nEns; $mean_amp1[$b] = $sum_amp1[$b] / $nEns; $mean_amp2[$b] = $sum_amp2[$b] / $nEns; $mean_amp3[$b] = $sum_amp3[$b] / $nEns; $mean_amp4[$b] = $sum_amp4[$b] / $nEns; - $mean_pcg1[$b] = $sum_pcg1[$b] / $n_pcg1[$b]; $mean_pcg2[$b] = $sum_pcg2[$b] / $n_pcg2[$b]; - $mean_pcg3[$b] = $sum_pcg3[$b] / $n_pcg3[$b]; $mean_pcg4[$b] = $sum_pcg4[$b] / $n_pcg4[$b]; - $mean_u[$b] = $sum_u[$b] / $good_vels[$b]; $mean_v[$b] = $sum_v[$b] / $good_vels[$b]; + + $mean_pcg1[$b] = $sum_pcg1[$b] / $n_pcg1[$b] if ($n_pcg1[$b] > 0); + $mean_pcg2[$b] = $sum_pcg2[$b] / $n_pcg2[$b] if ($n_pcg2[$b] > 0); + $mean_pcg3[$b] = $sum_pcg3[$b] / $n_pcg3[$b] if ($n_pcg3[$b] > 0); + $mean_pcg4[$b] = $sum_pcg4[$b] / $n_pcg4[$b] if ($n_pcg4[$b] > 0); + + next unless ($good_vels[$b] > 0); + + $mean_u[$b] = $sum_u[$b] / $good_vels[$b]; + $mean_v[$b] = $sum_v[$b] / $good_vels[$b]; $mean_w[$b] = $sum_w[$b] / $good_vels[$b]; - $mean_e[$b] = $sum_e[$b] / ($good_vels[$b] - $three_beam[$b]); + $mean_e[$b] = ($good_vels[$b] - $three_beam[$b] > 0) ? $sum_e[$b] / ($good_vels[$b] - $three_beam[$b]) : undef; } for ($e=$fe; $e<=$le; $e++) { @@ -210,13 +220,13 @@ $sumsq_amp4[$b] += ($mean_amp4[$b] - $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b][3])**2; $sumsq_pcg1[$b] += ($mean_pcg1[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0])**2 - if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]); + if defined($mean_pcg1[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]); $sumsq_pcg2[$b] += ($mean_pcg2[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1])**2 - if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]); + if defined($mean_pcg2[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][1]); $sumsq_pcg3[$b] += ($mean_pcg3[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2])**2 - if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]); + if defined($mean_pcg3[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][2]); $sumsq_pcg4[$b] += ($mean_pcg4[$b] - $dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3])**2 - if defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]); + if defined($mean_pcg4[$b]) && defined($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]); next unless ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]); @@ -226,7 +236,8 @@ next if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]); - $sumsq_e[$b] += ($mean_e[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])**2; + $sumsq_e[$b] += ($mean_e[$b] - $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3])**2 + if defined($mean_e[$b]); } } @@ -251,6 +262,8 @@ # Calculate Beam Statistics #---------------------------------------------------------------------- +# not implemented yet + #---------------------------------------------------------------------- # Produce Output #---------------------------------------------------------------------- @@ -260,6 +273,8 @@ ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) : 1; + print("#!/usr/bin/perl -S list\n"); + chmod(0777&~umask,*STDOUT); print("#ANTS#PARAMS# "); foreach my $k (keys(%P)) { print("$k\{$P{$k}\} ");