#!/usr/bin/perl
#======================================================================
# L I S T W
# doc: Wed Mar 24 06:45:09 2004
# dlm: Thu Jul 30 17:42:33 2009
# (c) 2004 A.M. Thurnherr
# uE-Info: 205 53 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# dump vertical velocities
# NB: currently broken
# HISTORY:
# Mar 24, 2004: - created from [listens] and [mkprofile]
# Mar 27, 2004: - added elapsed field
# - floatized time
# Apr 3, 2004: - cosmetics
# Nov 8, 2005: - UNIXTIME => UNIX_TIME
# Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
# Jul 30, 2009: - NaN => nan
$0 =~ m{(.*)/[^/]+};
require "$1/WorkhorseBinRead.pl";
require "$1/WorkhorseCoords.pl";
require "$1/WorkhorseUtils.pl";
use Getopt::Std;
$USAGE = "$0 @ARGV";
die("Usage: $0 " .
"[-A)nts] " .
"[-F)ilter <script>] " .
"[bin -r)ange <bin|0,bin|*>] " .
"[-e)rr-vel <max|0.1>] [-c)orrelation <min|70>] " .
"[-S)alin <val|35>] [-t)emp <bias>] " .
"[output -f)ields <field[,...]> " .
"<RDI file>\n")
unless (&getopts("Ac:e:F:f:r:S:t:") && @ARGV == 1);
$opt_e = 0.1 unless defined($opt_e); # defaults
$opt_c = 70 unless defined($opt_c);
$opt_S = 35 unless defined($opt_S);
print(STDERR "WARNING: Using uncalibrated ADCP temperature!\n"),$opt_t = 0
unless defined($opt_t);
require $opt_F if defined($opt_F); # load filter
if ($opt_f) { # additional fields
@f = split(',',$opt_f);
foreach $f (@f) {
my($f) = $f; # copy it
$f =~ s{\[.*$}{}; # remove indices
$addFields .= " {$f}";
}
}
#----------------------------------------------------------------------
print(STDERR "Reading $ARGV[0]..."); # read data
readData($ARGV[0],\%dta);
print(STDERR "done\n");
if (defined($opt_r)) { # bin range
($minb,$maxb) = split(',',$opt_r);
die("$0: can't decode -r $opt_r\n") unless defined($maxb);
} else {
$minb = 0;
$maxb = $dta{N_BINS} - 1;
}
die("$ARGV[0]: not enough bins for choice of -r\n") # enough bins?
unless ($dta{N_BINS} >= $maxb);
if ($dta{BEAM_COORDINATES}) { # coords used
$beamCoords = 1;
} elsif (!$dta{EARTH_COORDINATES}) {
die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
}
#----------------------------------------------------------------------
# Reference-Layer w (from [mkprofile])
# - also sets W field when valid
#----------------------------------------------------------------------
sub ref_lr_w($)
{
my($ens) = @_;
my($i,$n,$w);
for ($i=$minb; $i<=$maxb; $i++) {
next if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $opt_c ||
$dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $opt_c ||
$dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $opt_c ||
$dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $opt_c);
if ($beamCoords) {
next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100 ||
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100 ||
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100 ||
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
@v = velInstrumentToEarth(\%dta,$ens,
velBeamToInstrument(\%dta,
@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
} else {
next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
$dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
@v = @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]};
}
next if (!defined($v[3]) || abs($v[3]) > $opt_e);
if (defined($v[2])) { # valid w
$dta{ENSEMBLE}[$ens]->{W}[$i] = $v[2];
$w += $v[2]; $n++;
}
}
return $n ? $w/$n : undef;
}
#----------------------------------------------------------------------
print(STDERR "Generating profile by integrating w...");
for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
checkEnsemble(\%dta,$e); # sanity checks
filterEnsemble(\%dta,$e) # filter ensemble
if (defined($opt_F) &&
$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][0] > 0);
$dta{ENSEMBLE}[$e]->{REFW} = ref_lr_w($e);
next unless defined($dta{ENSEMBLE}[$e]->{REFW});
unless (defined($firstgood)) { # init profile
$firstgood = $lastgood = $e;
$depth = 0;
}
my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
$dta{ENSEMBLE}[$lastgood]->{UNIX_TIME}; # ... last good ens
if ($dt > 120) {
printf(STDERR "\nWARNING: %d-s gap too long, profile restarted at ensemble $e...",$dt);
$firstgood = $lastgood = $e;
$dt = $depth = 0;
}
$depth += $dta{ENSEMBLE}[$lastgood]->{REFW} * $dt # integrate depth
if ($dt > 0);
$dta{ENSEMBLE}[$e]->{DEPTH} = $depth;
$atbottom = $e, $maxdepth = $depth if ($depth > $maxdepth);
my($ss) = soundSpeed($opt_S, # sound-speed corr
$dta{ENSEMBLE}[$e]->{TEMPERATURE}-$opt_t,
$depth);
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} =
$ss / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND} = $ss;
$dta{ENSEMBLE}[$e]->{REFW} *=
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION};
$lastgood = $e;
}
printf(STDERR "done (max depth = %.1fm, depth at end of cast = %.1fm)\n",
$maxdepth,$depth);
filterEnsembleStats() if defined($opt_F);
#----------------------------------------------------------------------
print(STDERR "Writing output...");
if ($opt_A) {
print("#ANTS# [] $USAGE\n");
print("#ANTS#FIELDS# {ens} {unix_time} {time} {bin} {depth} {dz} {w} {ref_w} {dw} $addFields\n");
printf("#ANTS#PARAMS# maxdepth{$max_depth} bottom_time{%d}\n",
$dta{ENSEMBLE}[$atbottom]->{UNIX_TIME} -
$dta{ENSEMBLE}[$firstgood]->{UNIX_TIME});
} else {
print("# ens-no time elapsed bin-no depth dz w ref-w dw $addFields\n");
print("#----------------------------------------------------------------------\n");
}
for ($e=$firstgood; $e<=$lastgood; $e++) {
for ($i=$minb; $i<=$maxb; $i++) { # dump valid
next unless defined($dta{ENSEMBLE}[$e]->{W}[$i]);
printf("%d %f %f %d %.1f %.1f %g %g %g ",
$e,$dta{ENSEMBLE}[$e]->{UNIX_TIME},
$dta{ENSEMBLE}[$e]->{UNIX_TIME} -
$dta{ENSEMBLE}[$firstgood]->{UNIX_TIME},$i,
$dta{ENSEMBLE}[$e]->{DEPTH} +
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
$dta{ENSEMBLE}[$e]->{W}[$i],
$dta{ENSEMBLE}[$e]->{REFW},
$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
$dta{ENSEMBLE}[$e]->{W}[$i] - $dta{ENSEMBLE}[$e]->{REFW},
);
sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }
if (defined(@f)) {
foreach $f (@f) {
my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
$fn = $f unless defined($fn);
p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
}
}
print("\n");
}
}
print(STDERR "done\n");
exit(0);