diff --git a/listBins b/listBins new file mode 100755 --- /dev/null +++ b/listBins @@ -0,0 +1,330 @@ +#!/usr/bin/perl +#====================================================================== +# L I S T B I N S +# doc: Fri Aug 25 15:57:05 2006 +# dlm: Sat May 23 16:34:30 2009 +# (c) 2006 A.M. Thurnherr +# uE-Info: 269 0 NIL 0 0 72 2 2 4 NIL ofnI +#====================================================================== + +# Split data file into per-bin time series. + +# HISTORY: +# Aug 25, 2006: - created from [listEns] +# Aug 26, 2006: - added -M)agdecl +# - changed -b to -f +# Aug 27, 2006: - added %bin +# Aug 28, 2006: - BUG: options were confused +# Jan 4, 2007: - improved usage message for -a +# - added %mag_decl +# - BUG: roundoff error in %pct_good_vels +# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] +# Feb 7, 2008: - added sound-speed correction +# - enabled 3-beam solutions +# Feb 8, 2008: - added -d)iscard +# - added -b)eam coordinate output +# Feb 12, 2008: - modified 3-beam output +# - added -p)ct_good +# Feb 13, 2008: - various improvements +# Feb 19, 2008: - BUG: division by zero +# BUG: min() did not work with 1st elt undef +# Feb 21, 2008: - BUG: had forgotten to undo debugging changes +# - removed missing magdecl warning on -b +# Feb 22, 2008: - moved ssCorr() to [RDI_Utils.pl] +# - BUG: %d complete ensembles was written to STDOUT +# - BUG: %N_ensembles was total number of ensembles, not +# only the ones used +# - BUG: 0 errvel was output for 3-beam solutions => +# wrong time-average statistics +# May 19, 2009: - added -w to calculate vertical velocities from +# two beam pairs separately +# May 21, 2009: - added horizontal beampair velocities on -w +# - -P)itchRoll +# May 22, 2009: - added -B) +# May 23, 2009: - adapted to changed beampair-velocity fun name + +# General Notes: +# - everything (e.g. beams) is numbered from 1 +# - no support for BT data + +# Soundspeed Correction: +# - 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 +# - simple correction for cell depths + +# Min %-good (min_pcg): +# - nan for records w/o valid velocities +# - min(%-good) of the beams used for the velocity solution +# - min_pcg does not have to decrease monotonically with distance, +# at least when 3-beam solutions are allowed and when -p is used to +# edit the data +# - non-monotonic min_pcg is particularly obvious with the DYNAMUCK BM_ADCP +# data, where one of the beams performed much worse than the others + +require "getopts.pl"; +$0 =~ m{(.*/)[^/]+}; +require "$1RDI_BB_Read.pl"; +require "$1RDI_Coords.pl"; +require "$1RDI_Utils.pl"; + +die("Usage: $0 [-r)ange ] " . + "[output -f)ile ] " . + "[-a)ll ens (not just those with good vels)] " . + "[-M)agnetic ] " . + "[-S)oundspeed correction " . + "[-P)itch/Roll ] [-B)eamvel ] " . + "[require -4)-beam solutions] [-d)iscard ] " . + "[-%)good ] " . + "[output -b)eam coordinates] [output two separate -w) estimates] " . + "\n") + unless (&Getopts("4abB:d:f:M:p:r:P:S:w") && @ARGV == 1); + +($P{pitch_bias},$P{roll_bias}) = split('[,/]',$opt_P); +($P{velbias_b1},$P{velbias_b2},$P{velbias_b3},$P{velbias_b4}) = split('[,/]',$opt_B); + +die("$0: -4 and -d are mutually exclusive\n") + if ($opt_4 && defined($opt_d)); + +die("$0: -p and -b are mutually exclusive\n") + if ($opt_b && defined($opt_p)); + +$opt_p = 0 unless defined($opt_p); + +$RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions + +print(STDERR "WARNING: magnetic declination not set!\n") + unless defined($opt_M) || defined($opt_b); + +$opt_f = 'bin%d.raw' unless defined($opt_f); +$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); +$variable_ssCorr = ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*'); + +#---------------------------------------------------------------------- + +sub min(@) # return minimum +{ + my($min) = 99e99; + for (my($i)=0; $i<=$#_; $i++) { + $min = $_[$i] if defined($_[$i]) && ($_[$i] < $min); + } + return ($min == 99e99) ? nan : $min; +} + +sub dumpBin($$$) # write time series of single bin +{ + my($b,$fe,$le) = @_; + my($file) = sprintf($opt_f,$b+1); + + open(P,">$file") || die("$file: $!\n"); + print(P "#ANTS#PARAMS# "); + foreach my $k (keys(%P)) { + print(P "$k\{$P{$k}\} "); + } + if ($opt_b) { + printf(STDERR "%.0f%% ",100*$good_vels[$b]/($le-$fe+1)); + } else { + my($pct3b) = ($good_vels[$b] > 0) ? 100*$three_beam[$b]/$good_vels[$b] : nan; + printf(STDERR "%.0f%%/%.0f%% ",100*$good_vels[$b]/($le-$fe+1),$pct3b); + printf(P "pct_3_beam{%.0f} ",$pct3b); + } + printf(P "pct_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1)); + printf(P "bin{%d}",$b+1); + printf(P " soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!'); + printf(P " dz{%g}",$dz[$b] * + (defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) : 1) + ) unless ($variable_ssCorr); + print( P "\n"); + + print(P "#ANTS#FIELDS# " . + "{ensemble} {date} {time} {elapsed} " . + "{heading} {pitch} {roll} " . + "{sig_heading} {sig_pitch} {sig_roll} " . + "{xmit_current} {xmit_voltage} " . + "{temp} " . + ($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") . + ($opt_w ? "{v12} {w12} {v34} {w34} " : "" ) . + "{corr1} {corr2} {corr3} {corr4} " . + "{amp1} {amp2} {amp3} {amp4} " . + ($opt_b ? "{pcg1} {pcg2} {pcg3} {pcg4}" : "{3_beam} {min_pcg}") + ); + print(P " {dz}") if ($variable_ssCorr); + print(P "\n"); + + my($t0) = $dta{ENSEMBLE}[$fe]->{UNIX_TIME}; + for (my($e)=$fe; $e<=$le; $e++) { + next unless ($opt_a || $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]); + + my($ssCorr) = defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$e],$SS_salin,$SS_temp,$SS_depth) : 1; + + print(P "$dta{ENSEMBLE}[$e]->{NUMBER} "); + print(P "$dta{ENSEMBLE}[$e]->{DATE} "); + print(P "$dta{ENSEMBLE}[$e]->{TIME} "); + printf(P "%d ",$dta{ENSEMBLE}[$e]->{UNIX_TIME}-$t0); + print(P "$dta{ENSEMBLE}[$e]->{HEADING} "); + print(P "$dta{ENSEMBLE}[$e]->{PITCH} "); + print(P "$dta{ENSEMBLE}[$e]->{ROLL} "); + print(P "$dta{ENSEMBLE}[$e]->{HEADING_STDDEV} "); + print(P "$dta{ENSEMBLE}[$e]->{PITCH_STDDEV} "); + print(P "$dta{ENSEMBLE}[$e]->{ROLL_STDDEV} "); + print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT} "); + print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} "); + print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} "); + if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) { + printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] * $ssCorr); + printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] * $ssCorr); + printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] * $ssCorr); + if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]) { + print(P "nan "); + } else { + printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr); + } + if ($opt_w) { + printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]) ? "%g " : "nan ", + $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]); + printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]) ? "%g " : "nan ", + $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]); + printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]) ? "%g " : "nan ", + $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]); + printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]) ? "%g " : "nan ", + $dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]); + } + } else { + print(P "nan nan nan nan "); + print(P "nan nan nan nan ") if ($opt_w); + } + print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} "); + print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} "); + if ($opt_b) { + print(P "@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]} "); + } else { + printf(P "%d ",$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]); + printf(P "%s ",min(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]})); + } + printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr); + print(P "\n"); + } + close(P); +} + +#---------------------------------------------------------------------- +# MAIN +#---------------------------------------------------------------------- + +$P{RDI_file} = $ifn; +$P{mag_decl} = $opt_M if defined($opt_M); + +readData($ifn,\%dta); +printf(STDERR "%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}})); +$dta{HEADING_BIAS} = -$opt_M; # magnetic declination + +if ($dta{BEAM_COORDINATES}) { # coords + $beamCoords = 1; +} else { + die("$0: -b requires input in beam coordinates\n") + if ($opt_b); + die("$ifn: only beam and earth coordinates implemented so far\n") + if (!$dta{EARTH_COORDINATES}); +} + +for (my($b)=0; $b<$dta{N_BINS}; $b++) { # calc dz + $dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH}; +} + +$lastGoodBin = 0; +for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) { # check/transform velocities + next if (defined($first_ens) && + $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens); + + $dta{ENSEMBLE}[$e]->{PITCH} -= $P{pitch_bias}; + $dta{ENSEMBLE}[$e]->{ROLL} -= $P{roll_bias}; + + $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; + + die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") + if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4); + die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n") + if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR}); + die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n") + if ($dta{ENSEMBLE}[$e]->{LOW_GAIN}); + + for (my($b)=0; $b<$dta{N_BINS}; $b++) { + $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] -= $P{velbias_b1}; + $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] -= $P{velbias_b2}; + $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] -= $P{velbias_b3}; + $dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] -= $P{velbias_b4}; + if ($opt_w) { + @{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} = + velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}); + } + + if (defined($opt_d)) { + undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]); + undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]); + } + for (my($i)=0; $i<4; $i++) { + if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) { + undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]); + undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]); + } + } + @{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = $beamCoords + ? velInstrumentToEarth(\%dta,$e, + velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) + ) + : velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]}) + unless ($opt_b); + unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) { + undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}); + next; + } + $dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag; + $three_beam[$b] += $RDI_Coords::threeBeamFlag; + $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1; + $good_vels[$b]++; + $lastGoodBin = $b if ($b > $lastGoodBin); + $firstGoodEns = $e unless defined($firstGoodEns); + $lastGoodEns = $e; + } +} + +unless (defined($opt_r)) { + $fe = $firstGoodEns; + $le = $lastGoodEns; +} + +$P{N_ensembles} = $le - $fe + 1; + +$firstBin = 0; +$lastBin = $lastGoodBin; + +print( STDERR "Start : $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n"); +print( STDERR "End : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n"); +printf(STDERR "Bins : %d-%d\n",$firstBin+1,$lastBin+1); +if ($opt_b) { + print(STDERR "Good : "); +} else { + printf(STDERR "3-Beam : %d %d %d %d\n",$RDI_Coords::threeBeam_1, + $RDI_Coords::threeBeam_2, + $RDI_Coords::threeBeam_3, + $RDI_Coords::threeBeam_4); + + print(STDERR "Good/3-beam: "); +} +for ($b=$firstBin; $b<=$lastBin; $b++) { # generate output + dumpBin($b,$fe,$le); +} +print(STDERR "\n"); + +exit(0);