#!/usr/bin/perl
#======================================================================
# M K P R O F I L E
# doc: Sun Jan 19 18:55:26 2003
# dlm: Fri Oct 13 13:39:55 2017
# (c) 2003 A.M. Thurnherr
# uE-Info: 230 37 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# Make an LADCP Profile by Integrating W (similar to Firing's scan*).
# HISTORY:
# Jan 19, 2003: - written in order to test the RDI libs
# Jan 20, 2003: - added ensemble number
# Jan 21, 2003: - added horizontal integration
# Jan 22, 2003: - corrected magnetic declination
# Jan 23, 2003: - added -F)ilter
# Jan 24, 2003: - added more %PARAMs; started integration from 1st bin
# - added -g, -f, battery status
# Jan 25, 2003: - added more %PARAMs
# Feb 1, 2003: - BUG: bottom-track quality checking was bad
# Feb 8, 2003: - allowed for array-indices on -f
# Feb 9, 2003: - added 50% goodvelbin
# - removed unknown-field err on -f to allow -f W
# Feb 10, 2003: - changed initialization depth to 0m
# - changed %bottom_depth to %max_depth
# Feb 11, 2003: - changed sign of magnetic declination
# Feb 12, 2003: - corrected BT-range scaling
# Feb 14, 2003: - added %pinging_hours, %min_range
# - removed magnetic declination from default
# Feb 26, 2004: - added earth coordinates
# Mar 3, 2004: - removed requirement for -M on !-Q
# - corrected range-stats on earth coordinates
# Mar 4, 2004: - added number of ensebles to output
# Mar 11, 2004: - BUG: rename ACD -> ADC
# Mar 12, 2004: - added %bottom_xmit_{current|voltage}
# Mar 16, 2004: - BUG: on -M u/v/x/y were wrong
# Mar 17, 2004: - added error estimates on u/v/x/y
# - removed battery stuff (has to be done btw casts)
# Mar 18, 2004: - totally re-did u/v integration
# Mar 19, 2004: - re-designed u/v uncertainty estimation
# Mar 28, 2004: - added MEAN_CORRELATION, MEAN_ECHO_AMPLITUDE
# Sep 15, 2005: - changed BinRead library name
# - made max gap length variable
# Sep 16, 2005: - re-did u,v,w uncertainties
# Nov 8, 2005: - UNIXTIME => UNIX_TIME
# - added unix_time, secno, z_BT to default output
# Dec 1, 2005: - moved profile-building code to [RDI_utils.pl]
# - changed -f syntax to allow name=FIELD
# - added %bin1_dist, %bin_length
# Dec 8, 2005: - remove spaces from -f argument to allow multiline
# definitions in Makefiles
# Nov 13, 2006: - BUG: end-of-cast depth had not been reported correctly
# - cosmetics
# Nov 30, 2007: - adapted to 3-beam solutions
# Dec 11, 2007: - adapted to earlier modifications (Sep 2007) of
# [RDI_BB_Read.pl]
# Dec 14, 2007: - replaced z by depth
# Dec 17, 2007: - BUG: downcast flag was set incorrectly
# Jan 24, 2008: - rotation had been output as degrees/s; to make it more
# consistent with pitch/roll, I changed it to simple degrees
# - added net rotations [deployment]/down/up/[recovery]
# Apr 9, 2008: - added profile -B)ottom depth
# - BUG: depth of first bin was reported as beginning of cast
# Oct 24, 2008: - added RANGE and RANGE_BINS fields
# Mar 18, 2009: - BUG: pitch/roll calculation had typo
# - calc pitch/roll separately for down-/upcasts
# - removed approximations in pitch/roll calcs
# Jul 30, 2009: - typo '<' removed from output
# - NaN => nan
# Dec 8, 2010: - added zmax/zend labels to output
# Dec 10, 2010: - made mkProfile exit with status 0 if no good ens found but -Q is set
# Dec 19, 2010: - finally made -A default and activated output file
# Jan 5, 2011: - made no-good-ensembles found test much more robust
# Jun 22, 2011: - added bandwith/power warnings
# - added ping-interval calculation
# - BUG: post-recovery rotations were always zero
# Sep 9, 2011: - BUG: range calculation for Earth coordinate data included bins without
# valid velocities
# Sep 21, 2010: - added %rms_heave_acceleration
# Apr 12, 2013: - added -p
# May 10, 2013: - BUG: mkProfile bombed when ADCP file is truncated at deepest location
# May 14, 2013: - added heading to output
# - added err_vel to output
# - finally removed -d/-g
# Nov 25, 2013: - expunged checkEnsemble
# Feb 13, 2014: - added support set_range_lim()
# Mar 4, 2014: - added support to allow missing PITCH/ROLL/HEADING values
# May 24, 2014: - finally added (gimbal-)pitch & roll to default output
# - renamed heading to hdg and pitchroll to tilt
# Mar 22, 2015: - made it work for moored time series as well
# Mar 17, 2015: - adapted to new Getopt library
# - removed warning
# Sep 12, 2016: - added %PD0_file
# Oct 13, 2017: - added instrument orientation
# NOTES:
# - the battery values are based on transmission voltages (different
# from battery voltages) and reported without units (raw 8-bit a2d
# values)
# - -B with the CTD max depth can be used to linearly scale the depths;
# even so, the profile can have negative depths, in particular when
# the CTD is sent to a shallow depth first and then returned to the surface
# before beginning the cast
# - in one case that I looked at (Anslope ][, cast 82), there are large
# depth errors, even when -B is used
# - this utility works only approximately for uplookers (profile is
# roughly ok, but apparently contaminated by surface reflection,
# but stats are not ok; e.g. NBP0402 037U.prof)
$0 =~ m{(.*)/[^/]+};
require "$1/RDI_BB_Read.pl";
require "$1/RDI_Coords.pl";
require "$1/RDI_Utils.pl";
use Getopt::Std;
$USAGE = "$0 @ARGV";
die("Usage: $0 " .
"[-Q)uiet] [-F)ilter <script>] " .
"[require -4)-beam solutions] [-d)iscard <beam#>] [apply beamvel-m)ask <file>] " .
"[-r)ef-layer <bin|1,bin|6>] [-n) vels <min|2>] " .
"[-e)rr-vel <max[0.1]] [-c)orrelation <min>] [-p)ct-good <min[100]>] " .
"[max -g)ap <len>] " .
"[output -f)ields <field[,...]> " .
"[-M)agnetic <declination>] [profile -B)ottom <depth>] " .
"<RDI file>\n")
unless (&getopts("4AB:F:M:Qd:g:r:n:e:c:f:m:p:") && @ARGV == 1);
$RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions
require $opt_F if defined($opt_F); # load filter
$opt_r = "1,6" unless defined($opt_r); # defaults
$opt_n = 2 unless defined($opt_n);
$opt_e = 0.1 unless defined($opt_e);
$opt_c = 70 unless defined($opt_c);
$opt_g = 120 unless defined($opt_g);
$opt_p = 100 unless defined($opt_p);
($minb,$maxb) = split(',',$opt_r); # reference layer
die("$0: can't decode -r $opt_r\n") unless defined($maxb);
print(STDERR "Reading $ARGV[0]..."); # read data
readData($ARGV[0],\%dta);
print(STDERR "done\n");
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");
}
if (defined($opt_m) && -r $opt_m) {
die("$ARGV[0]: -m only implemented for data collected in beam coordinates\n")
unless ($beamCoords);
print(STDERR "Masking beam velocities as prescribed in $opt_m...");
open(BVM,$opt_m) || die("$opt_m: $!\n");
while (<BVM>) {
s/#.*//;
s/^\s*$//;
next if ($_ eq '');
my($fe,$te,$db) = split;
die("$opt_m: cannot decode $_\n")
unless (numberp($fe) && numberp($te) && $te>=$fe && $db>=1 && $db<=4);
die("$0: assertion failed")
unless ($dta{ENSEMBLE}[$fe-1]->{NUMBER} == $fe &&
$dta{ENSEMBLE}[$te-1]->{NUMBER} == $te);
for (my($ens)=$fe-1; $ens<=$te-1; $ens++) {
$nens++;
for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
undef($dta{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$db-1]);
}
}
}
close(BVM);
print(STDERR " $nens ensembles edited\n");
}
if (defined($opt_d)) { # discard entire beam
die("$ARGV[0]: -d only implemented for data collected in beam coordinates\n")
unless ($beamCoords);
print(STDERR "Discarding beam-$opt_d velocities...");
for (my($ens)=0; $ens<=$#{$dta{ENSEMBLE}}; $ens++) {
for (my($bin)=0; $bin<$dta{N_BINS}; $bin++) {
undef($dta{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$opt_d-1]);
}
}
print(STDERR "done\n");
}
if (defined($opt_M)) { # magnetic declination
$dta{HEADING_BIAS} = -1*$opt_M;
} else {
$dta{HEADING_BIAS} = 0;
}
ensure_BT_RANGE(\%dta); # set BT_RANGE field if it is missing (old firmware bug)
if ($opt_f) { # additional fields
@f = split(',',$opt_f);
foreach $f (@f) {
$f =~ s/\s//g; # remove spaces
@def = split('=',$f);
if (@def == 2) { # name=field
$addFields .= " {$def[0]}";
$f = $def[1];
} else { # field
$addFields .= " {$f}";
}
}
# print(STDERR "addFields = $addFields\n");
# print(STDERR "\@f = @f\n");
}
#======================================================================
# Step 0: Check data & Calculate Ping Rates
#======================================================================
unless ($dta{NARROW_BANDWIDTH}) {
print(STDERR "WARNING: $0 WIDE BANDWIDTH!\n");
}
unless ($dta{TRANSMIT_POWER_HIGH}) {
print(STDERR "WARNING: $0 LOW TRANSMIT POWER!\n");
}
printf(STDERR "ADCP : %s (s/n %d) %s\n",
$dta{INSTRUMENT_TYPE},$dta{SERIAL_NUMBER},
$dta{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'UL' : 'DL');
printf(STDERR "# of ensembles : %d\n",scalar(@{$dta{ENSEMBLE}}));
my($sdt1,$sdt2,$ndt);
my($mindt1) = my($mindt2) = 9e99;
my($maxdt1) = my($maxdt2) = 0;
for (my($e)=2; $e<=$#{$dta{ENSEMBLE}}; $e+=2,$ndt++) {
my($dt1) = $dta{ENSEMBLE}[$e-1]->{UNIX_TIME} - $dta{ENSEMBLE}[$e-2]->{UNIX_TIME};
my($dt2) = $dta{ENSEMBLE}[$e-0]->{UNIX_TIME} - $dta{ENSEMBLE}[$e-1]->{UNIX_TIME};
$mindt1 = $dt1 if ($dt1 < $mindt1);
$mindt2 = $dt2 if ($dt2 < $mindt2);
$maxdt1 = $dt1 if ($dt1 > $maxdt1);
$maxdt2 = $dt2 if ($dt2 > $maxdt2);
$sdt1 += $dt1; $sdt2 += $dt2;
}
printf(STDERR "Ping intervals : %.1fs/%.1fs",$sdt1/$ndt,$sdt2/$ndt);
if ($maxdt1-$mindt1>=0.1 || $maxdt2-$mindt2>=0.1) {
printf(STDERR " (%.1fs-%.1fs/%.1fs-%.1fs)\n",$mindt1,$maxdt1,$mindt2,$maxdt2);
} else {
print(STDERR "\n");
}
#======================================================================
# Step 1: Integrate w & determine water depth
#======================================================================
($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz,$rms_heave_accel) =
mk_prof(\%dta,0,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_g,$opt_p);
if ($lastgood == $atbottom) {
print(STDERR "$ARGV[0]: truncated file (ends at max depth)\n")
} elsif (($atbottom > $firstgood) && ($lastgood > $atbottom)) {
# all good
} elsif ($lastgood > $firstgood) {
print(STDERR "$ARGV[0]: no bottom depth found\n")
unless ($atbottom > 0);
} else {
if ($opt_Q) {
print(STDERR "$ARGV[0]: no valid cast data found\n");
exit(0);
} else {
die(sprintf("$ARGV[0]: no valid cast data found (firstgood=%d atbottom=%d lastgood=%d)\n",$firstgood,$atbottom,$lastgood));
}
}
if (defined($opt_B)) { # scale Z
my($zscale) = $opt_B / ($dta{ENSEMBLE}[$atbottom]->{DEPTH} -# downcast
$dta{ENSEMBLE}[$firstgood]->{DEPTH});
# printf(STDERR "scaling downcast depths by %.2f\n",$zscale);
for (my($e)=$firstgood; $e<$atbottom; $e++) {
next unless defined($dta{ENSEMBLE}[$e]->{DEPTH});
$dta{ENSEMBLE}[$e]->{DEPTH} =
$dta{ENSEMBLE}[$firstgood]->{DEPTH} + $zscale *
($dta{ENSEMBLE}[$e]->{DEPTH}-$dta{ENSEMBLE}[$firstgood]->{DEPTH});
}
$zscale = $opt_B / ($dta{ENSEMBLE}[$atbottom]->{DEPTH} - # upcast
$dta{ENSEMBLE}[$lastgood]->{DEPTH});
# printf(STDERR "scaling upcast depths by %.2f\n",$zscale);
for (my($e)=$atbottom; $e<=$lastgood; $e++) {
next unless defined($dta{ENSEMBLE}[$e]->{DEPTH});
$dta{ENSEMBLE}[$e]->{DEPTH} =
$dta{ENSEMBLE}[$firstgood]->{DEPTH} + $zscale *
($dta{ENSEMBLE}[$e]->{DEPTH}-$dta{ENSEMBLE}[$lastgood]->{DEPTH});
}
}
set_range_lim(\%dta); # set {range_lim} field
($water_depth,$sig_wd) = # sea bed
find_seabed(\%dta,$atbottom,$beamCoords);
#======================================================================
# Step 2: Integrate u & v
#======================================================================
sub ref_lr_uv($$$) # calc ref-level u/v
{
my($ens,$z,$water_depth) = @_;
my($i,$n,@v,@goodU,@goodV);
$water_depth = 99999 unless defined($water_depth);
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 = velApplyHdgBias(\%dta,$ens,
@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]});
}
next if (!defined($v[3]) || abs($v[3]) > $opt_e);
# Martin's BT routines show strong shear just above sea bed
# => skip lowest 20m.
if (defined($v[0])) { # valid u,v
if ($dta{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) {
if ($z - $dta{DISTANCE_TO_BIN1_CENTER}
- $i*$dta{BIN_LENGTH} > 0) {
push(@goodU,$v[0]); push(@goodV,$v[1]);
$dta{ENSEMBLE}[$ens]->{U} += $v[0];
$dta{ENSEMBLE}[$ens]->{V} += $v[1];
$n++;
}
} else {
if ($z + $dta{DISTANCE_TO_BIN1_CENTER}
+ $i*$dta{BIN_LENGTH} < $water_depth-20) {
push(@goodU,$v[0]); push(@goodV,$v[1]);
$dta{ENSEMBLE}[$ens]->{U} += $v[0];
$dta{ENSEMBLE}[$ens]->{V} += $v[1];
$n++;
}
}
}
}
if ($n >= 2) {
my(@sumsq) = (0,0);
$dta{ENSEMBLE}[$ens]->{U} /= $n;
$dta{ENSEMBLE}[$ens]->{V} /= $n;
for ($i=0; $i<$n; $i++) {
$sumsq[0] += ($dta{ENSEMBLE}[$ens]->{U}-$goodU[$i])**2;
$sumsq[1] += ($dta{ENSEMBLE}[$ens]->{V}-$goodV[$i])**2;
}
$dta{ENSEMBLE}[$ens]->{U_ERR} = sqrt($sumsq[0])/($n-1);
$dta{ENSEMBLE}[$ens]->{V_ERR} = sqrt($sumsq[1])/($n-1);
} else {
$dta{ENSEMBLE}[$ens]->{U} = undef;
$dta{ENSEMBLE}[$ens]->{V} = undef;
}
}
#----------------------------------------------------------------------
($x,$y) = (0,0); # init
$dta{ENSEMBLE}[$firstgood]->{X} = $dta{ENSEMBLE}[$firstgood]->{X_ERR} = 0;
$dta{ENSEMBLE}[$firstgood]->{Y} = $dta{ENSEMBLE}[$firstgood]->{Y_ERR} = 0;
$prevgood = $firstgood;
for ($e=$firstgood+1; defined($opt_M)&&$e<=$lastgood; $e++) {
#--------------------------------------------------
# within profile: both $firstgood and $prevgood set
#--------------------------------------------------
ref_lr_uv($e,$dta{ENSEMBLE}[$e]->{DEPTH},$water_depth) # instrument vel
if (defined($dta{ENSEMBLE}[$e]->{W}));
if (!defined($dta{ENSEMBLE}[$e]->{U})) { # gap
$uv_gap_time += $dta{ENSEMBLE}[$e]->{UNIX_TIME} -
$dta{ENSEMBLE}[$e-1]->{UNIX_TIME};
next;
}
my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} - # time step since
$dta{ENSEMBLE}[$prevgood]->{UNIX_TIME}; # ...last good ens
#-----------------------------------
# The current ensemble has valid u/v
#-----------------------------------
$x -= $dta{ENSEMBLE}[$prevgood]->{U} * $dt; # integrate
$xErr += ($dta{ENSEMBLE}[$prevgood]->{U_ERR} * $dt)**2;
$dta{ENSEMBLE}[$e]->{X} = $x;
$dta{ENSEMBLE}[$e]->{X_ERR} = sqrt($xErr);
$y -= $dta{ENSEMBLE}[$prevgood]->{V} * $dt;
$yErr += ($dta{ENSEMBLE}[$prevgood]->{V_ERR} * $dt)**2;
$dta{ENSEMBLE}[$e]->{Y} = $y;
$dta{ENSEMBLE}[$e]->{Y_ERR} = sqrt($yErr);
$prevgood = $e;
}
unless (defined($dta{ENSEMBLE}[$lastgood]->{X})) { # last is bad in u/v
my($dt) = $dta{ENSEMBLE}[$lastgood]->{UNIX_TIME} - # time step since
$dta{ENSEMBLE}[$prevgood]->{UNIX_TIME}; # ...last good ens
$x -= $dta{ENSEMBLE}[$prevgood]->{U} * $dt; # integrate
$xErr += ($dta{ENSEMBLE}[$prevgood]->{U_ERR} * $dt)**2;
$dta{ENSEMBLE}[$lastgood]->{X} = $x;
$dta{ENSEMBLE}[$lastgood]->{X_ERR} = sqrt($xErr);
$y -= $dta{ENSEMBLE}[$prevgood]->{V} * $dt;
$yErr += ($dta{ENSEMBLE}[$prevgood]->{V_ERR} * $dt)**2;
$dta{ENSEMBLE}[$lastgood]->{Y} = $y;
$dta{ENSEMBLE}[$lastgood]->{Y_ERR} = sqrt($yErr);
}
$firstgood++ if ($firstgood == 0); # centered diff
$lastgood-- if ($lastgood == $#{$dta{ENSEMBLE}}); # in step 6
#======================================================================
# Step 3: Calculate Uncertainties
#======================================================================
# Time series of W_ERR indicate that errors are very large near the
# surface and near the sea bed, perhaps because of reflections.
# A reasonable estimate for typical uncertainty is therefore the mode
# of the std errors.
my(@histUErr,@histVErr,@histWErr);
my($histRez) = 1e-4;
for ($e=$firstgood; $e<=$lastgood; $e++) {
$histWErr[int($dta{ENSEMBLE}[$e]->{W_ERR}/$histRez+0.5)]++
if defined($dta{ENSEMBLE}[$e]->{W_ERR});
$histUErr[int($dta{ENSEMBLE}[$e]->{U_ERR}/$histRez+0.5)]++
if defined($dta{ENSEMBLE}[$e]->{U_ERR});
$histVErr[int($dta{ENSEMBLE}[$e]->{V_ERR}/$histRez+0.5)]++
if defined($dta{ENSEMBLE}[$e]->{V_ERR});
}
my($max) = 0; my($mode);
for (my($i)=0; $i<=$#histWErr; $i++) {
next if ($histWErr[$i] < $max);
$max = $histWErr[$i]; $mode = $i;
}
$wErr = $mode * $histRez if defined($mode);
$max = 0; $mode = undef;
for (my($i)=0; $i<=$#histUErr; $i++) {
next if ($histUErr[$i] < $max);
$max = $histUErr[$i]; $mode = $i;
}
$uErr = $mode * $histRez if defined($mode);
$max = 0; $mode = undef;
for (my($i)=0; $i<=$#histVErr; $i++) {
next if ($histVErr[$i] < $max);
$max = $histVErr[$i]; $mode = $i;
}
$vErr = $mode * $histRez if defined($mode);
#print(STDERR "u: mu = $meanU / sigma = $uErr\n");
#print(STDERR "v: mu = $meanV / sigma = $vErr\n");
#print(STDERR "w: mu = $meanW / sigma = $wErr\n");
if (defined($opt_M)) { # displacement errors
$x_err = $uErr * $uv_gap_time + $dta{ENSEMBLE}[$lastgood]->{X_ERR};
$y_err = $vErr * $uv_gap_time + $dta{ENSEMBLE}[$lastgood]->{Y_ERR};
}
$z_err = $wErr * $w_gap_time + $dta{ENSEMBLE}[$lastgood]->{DEPTH_ERR};
#printf(STDERR "x_err = $dta{ENSEMBLE}[$lastgood]->{X_ERR} + %g\n",
# $uErr * $uv_gap_time);
#printf(STDERR "y_err = $dta{ENSEMBLE}[$lastgood]->{Y_ERR} + %g\n",
# $vErr * $uv_gap_time);
#printf(STDERR "z_err = $dta{ENSEMBLE}[$lastgood]->{DEPTH_ERR} + %g\n",
# $wErr * $w_gap_time);
#======================================================================
# Step 4: Calculate Beam Range Stats
#======================================================================
my($min_good_bins) = 999;
my($worst_beam);
sub count_good_vels($) # count good vels
{
my($ens) = @_;
my($good) = -1; my($this_worst_beam);
if ($beamCoords) {
for (my($i)=0; $i<$dta{N_BINS}; $i++) {
for (my($b)=0; $b<4; $b++) {
$good=$i,$this_worst_beam=$b,$nVels[$i][$b]++
if defined($dta{ENSEMBLE}[$ens]->{VELOCITY}[$i][$b]);
}
}
} else {
for (my($i)=0; $i<$dta{N_BINS}; $i++) {
next unless defined($dta{ENSEMBLE}[$ens]->{VELOCITY}[$i][0]);
for (my($b)=0; $b<4; $b++) {
$good=$i,$this_worst_beam=$b,$nVels[$i][$b]++
if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][$b] >=
$dta{MIN_CORRELATION});
}
}
}
$min_good_ens=$ens, $min_good_bins=$good, $worst_beam=$this_worst_beam
if ((!defined($water_depth) ||
$dta{ENSEMBLE}[$ens]->{DEPTH} < $water_depth-200)
&& $good >= 0 && $good < $min_good_bins);
}
#----------------------------------------------------------------------
for ($e=$firstgood; $e<=$lastgood; $e++) { # range
my($i);
for ($i=0; $i<$dta{N_BINS}; $i++) {
last if (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][0]) +
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][1]) +
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][2]) +
defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][3]) < 3);
}
$dta{ENSEMBLE}[$e]->{RANGE_BINS} = $i;
$dta{ENSEMBLE}[$e]->{RANGE} =
$dta{DISTANCE_TO_BIN1_CENTER} + $i * $dta{BIN_LENGTH};
}
for ($e=$firstgood; $e<=$lastgood; $e++) { # mean corr/amp
$sumcor = $sumamp = $ndata = 0;
for (my($i)=0; $i<$dta{N_BINS}; $i++) {
for (my($b)=0; $b<4; $b++) {
next unless ($dta{ENSEMBLE}[$e]->{CORRELATION}[$i][$b]);
$sumcor += $dta{ENSEMBLE}[$e]->{CORRELATION}[$i][$b];
$sumamp += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b];
$ndata++;
}
}
$dta{ENSEMBLE}[$e]->{MEAN_CORRELATION} = $sumcor/$ndata;
$dta{ENSEMBLE}[$e]->{MEAN_ECHO_AMPLITUDE} = $sumamp/$ndata;
}
for ($e=$firstgood+50; $e<=$lastgood-50; $e++) { # range stats
count_good_vels($e);
}
for ($i=0; $i<$dta{N_BINS}; $i++) {
for ($b=0; $b<4; $b++) {
$maxVels = $nVels[$i][$b] unless ($maxVels > $nVels[$i][$b]);
}
}
for ($i=0; $i<$dta{N_BINS}; $i++) {
for ($b=0; $b<4; $b++) {
$gb[$b] = $i if ($nVels[$i][$b] >= 0.8*$maxVels);
}
}
$gb = ($gb[0]+$gb[1]+$gb[2]+$gb[3]) / 4;
#======================================================================
# Step 5: Remove Ship Drift (probably not useful => removed)
#======================================================================
#======================================================================
# Step 6: Pitch, Roll, Rotation
#======================================================================
# in case of PITCH/ROLL/HEADING data gaps (IMP data), the calculations
# are not entirely correct, as
# i) the rotation implied by the pre-/post-gap headings is not counted
# ii) the gappy ensembles are counted for calculating the rms vals
my($prrms,$dnprrms,$upprrms) = (0,0,0);
my($rotrms,$prerot,$dnrot,$uprot,$postrot) = (0,0,0,0,0);
sub rot($)
{
my($e) = @_;
return 0
unless defined($dta{ENSEMBLE}[$e]->{HEADING}) && defined($dta{ENSEMBLE}[$e-1]->{HEADING});
my($rot) = $dta{ENSEMBLE}[$e]->{HEADING} -
$dta{ENSEMBLE}[$e-1]->{HEADING};
$rot -= 360 if ($rot > 180);
$rot += 360 if ($rot < -180);
return $rot;
}
for ($e=1; $e<$firstgood; $e++) { # pre-deployment
$prerot += rot($e);
}
for (; $e<= $atbottom; $e++) { # downcast
$dta{ENSEMBLE}[$e]->{TILT} =
&angle_from_vertical($dta{ENSEMBLE}[$e]->{PITCH},
$dta{ENSEMBLE}[$e]->{ROLL});
$prrms += $dta{ENSEMBLE}[$e]->{TILT}**2
if numberp($dta{ENSEMBLE}[$e]->{TILT});
$dta{ENSEMBLE}[$e]->{ROTATION} = rot($e);
$dnrot += $dta{ENSEMBLE}[$e]->{ROTATION};
$rotrms += $dta{ENSEMBLE}[$e]->{ROTATION}**2;
}
$dnprrms = $prrms;
for (; $e<=$lastgood; $e++) { # upcast
$dta{ENSEMBLE}[$e]->{TILT} =
&angle_from_vertical($dta{ENSEMBLE}[$e]->{PITCH},
$dta{ENSEMBLE}[$e]->{ROLL});
$prrms += $dta{ENSEMBLE}[$e]->{TILT}**2
if numberp($dta{ENSEMBLE}[$e]->{TILT});
$dta{ENSEMBLE}[$e]->{ROTATION} = rot($e);
$uprot += $dta{ENSEMBLE}[$e]->{ROTATION};
$rotrms += $dta{ENSEMBLE}[$e]->{ROTATION}**2;
}
$upprrms = $prrms - $dnprrms;
for (; $e<=$#{$dta{ENSEMBLE}}; $e++) { # post-recovery
$postrot += rot($e);
}
$prerot /= 360; # rotations, not degrees
$dnrot /= 360;
$uprot /= 360;
$postrot /= 360;
$prrms = sqrt($prrms/($lastgood-$firstgood));
$dnprrms = sqrt($dnprrms/($atbottom-$firstgood));
if ($lastgood == $atbottom) {
print(STDERR "WARNING: $0 NO UPCAST DATA\n");
$upprrms = nan;
} else {
$upprrms = sqrt($upprrms/($lastgood-$atbottom));
}
$rotrms = sqrt($rotrms/($lastgood-$firstgood));
#======================================================================
# PRODUCE OUTPUT
#======================================================================
printf(STDERR "Start of cast : %s (#%5d) at %6.1fm\n",
$dta{ENSEMBLE}[$firstgood]->{TIME},
$dta{ENSEMBLE}[$firstgood]->{NUMBER},
$dta{ENSEMBLE}[$firstgood]->{DEPTH});
printf(STDERR "Bottom of cast (zmax) : %s (#%5d) at %6.1fm\n",
$dta{ENSEMBLE}[$atbottom]->{TIME},
$dta{ENSEMBLE}[$atbottom]->{NUMBER},
$dta{ENSEMBLE}[$atbottom]->{DEPTH});
if (defined($water_depth)) {
printf(STDERR "Seabed : at %6.1fm (+-%dm)\n",$water_depth,$sig_wd);
} else {
print(STDERR "Seabed : not found\n");
}
printf(STDERR "End of cast (zend) : %s (#%5d) at %6.1fm\n",
$dta{ENSEMBLE}[$lastgood]->{TIME},
$dta{ENSEMBLE}[$lastgood]->{NUMBER},
$dta{ENSEMBLE}[$lastgood]->{DEPTH});
printf(STDERR "Rel. Displacement : x = %d(%d)m / y = %d(%d)m\n",
$dta{ENSEMBLE}[$lastgood]->{X}, $x_err,
$dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
) if defined($opt_M);
printf(STDERR "Cast Duration : %.1f hours (pinging for %.1f hours)\n",
$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME} / 3600,
($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{UNIX_TIME} -
$dta{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);
printf(STDERR "Minimum range : %dm at ensemble %d, beam %d\n",
$dta{DISTANCE_TO_BIN1_CENTER} +
$min_good_bins*$dta{BIN_LENGTH},
$dta{ENSEMBLE}[$min_good_ens]->{NUMBER},
$worst_beam);
printf(STDERR "80%%-valid bins : %.1f\n",$gb+1);
printf(STDERR "80%%-valid range : %dm\n",
$dta{DISTANCE_TO_BIN1_CENTER} + $gb*$dta{BIN_LENGTH});
printf(STDERR "3-beam solutions : $RDI_Coords::threeBeam_1 " .
"$RDI_Coords::threeBeam_2 " .
"$RDI_Coords::threeBeam_3 " .
"$RDI_Coords::threeBeam_4\n")
unless ($opt_4);
printf(STDERR "net rotations : [%d]/%d/%d/[%d]\n",$prerot,$dnrot,$uprot,$postrot);
printf(STDERR "rms pitch/roll : %.1fdeg/%.1fdeg\n",$dnprrms,$upprrms);
printf(STDERR "rms heave acceleration: %.2fm/s^2\n",$rms_heave_accel);
exit(0) if ($opt_Q);
#----------------------------------------------------------------------
# output profile in active ANTS format
#----------------------------------------------------------------------
print("#!/usr/bin/perl -S list\n");
chmod(0777&~umask,*STDOUT);
print("#ANTS# [] $USAGE\n");
$uFields = "{u} {u_err} {v} {v_err} {x} {x_err} {y} {y_err}"
if defined($opt_M);
print("#ANTS#FIELDS# {ens} {time} {elapsed} {secno} {downcast} " .
"{w} {w_err} {err_vel} {depth} {depth_err} {seabed} " .
"{pitch} {roll} {tilt} {hdg} {rotation} " .
"$uFields $addFields\n");
printf("#ANTS#PARAMS# date{$dta{ENSEMBLE}[$firstgood]->{DATE}} " .
"start_time{$dta{ENSEMBLE}[$firstgood]->{TIME}} " .
"bottom_time{$dta{ENSEMBLE}[$atbottom]->{TIME}} " .
"end_time{$dta{ENSEMBLE}[$lastgood]->{TIME}} " .
"bottom_xmit_voltage{$dta{ENSEMBLE}[$atbottom]->{ADC_XMIT_VOLTAGE}} " .
"bottom_xmit_current{$dta{ENSEMBLE}[$atbottom]->{ADC_XMIT_CURRENT}} " .
"PD0_file{$ARGV[0]} " .
"pinging_duration{%.1f} " .
"cast_duration{%.1f} " .
"0.8_valid_bins{%.1f} " .
"0.8_valid_range{%.1f} " .
"max_depth{%.1f} " .
"depth_error{%.1f} " .
"min_range{%d} " .
"n_ensembles{%d} " .
"w_gap_time{%d} " .
"stderr_w{%.4f} " .
"rms_tilt{%.1f} " .
"downcast_rms_tilt{%.1f} " .
"upcast_rms_tilt{%.1f} " .
"rms_rotation{%.2f} " .
"deployment_rotations{%d} " .
"downcast_rotations{%d} " .
"upcast_rotations{%d} " .
"recovery_rotations{%d} " .
"rms_heave_acceleration{%.2f} " .
"bin1_dist{%.1f} " .
"bin_length{%.1f} " .
"\n",
($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{UNIX_TIME} -
$dta{ENSEMBLE}[0]->{UNIX_TIME}),
$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
$gb+1,
$dta{DISTANCE_TO_BIN1_CENTER} + $gb*$dta{BIN_LENGTH},
$dta{ENSEMBLE}[$atbottom]->{DEPTH},
$dta{ENSEMBLE}[$lastgood]->{DEPTH} -
$dta{ENSEMBLE}[$firstgood]->{DEPTH},
$dta{DISTANCE_TO_BIN1_CENTER} +
$min_good_bins*$dta{BIN_LENGTH},
scalar(@{$dta{ENSEMBLE}}),
$w_gap_time,$wErr,$prrms,$dnprrms,$upprrms,$rotrms,
$prerot,$dnrot,$uprot,$postrot,$rms_heave_accel,
$dta{DISTANCE_TO_BIN1_CENTER},
$dta{BIN_LENGTH},
);
printf("#ANTS#PARAMS# magnetic_declination{$opt_M} " .
"uv_gap_time{%d} " .
"mean_u{%.4f} " .
"stderr_u{%.4f} " .
"dx{%d} " .
"dx_err{%d} " .
"mean_v{%.4f} " .
"stderr_v{%.4f} " .
"dy{%d} " .
"dy_err{%d}\n",
$uv_gap_time,
$dta{ENSEMBLE}[$lastgood]->{X} /
$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
$uErr, $dta{ENSEMBLE}[$lastgood]->{X}, $x_err,
$dta{ENSEMBLE}[$lastgood]->{Y} /
$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
$vErr, $dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
) if defined ($opt_M);
if ($dta{TIME_BETWEEN_PINGS} == 0) {
print("#ANTS#PARAMS# pinging_rate{staggered}\n");
} else {
printf("#ANTS#PARAMS# pinging_rate{%.2f}\n",
1/$dta{TIME_BETWEEN_PINGS});
}
if (defined($water_depth)) {
printf("#ANTS#PARAMS# water_depth{%d} sig-water_depth{%d}\n",
$water_depth,$sig_wd);
} else {
print("#ANTS#PARAMS# water_depth{nan} sig-water_depth{nan}\n");
}
sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }
sub pb($) { print($_[0]?"1 ":"0 "); }
for ($e=$firstgood; $e<=$lastgood; $e++) {
p($dta{ENSEMBLE}[$e]->{NUMBER});
p($dta{ENSEMBLE}[$e]->{UNIX_TIME});
p($dta{ENSEMBLE}[$e]->{ELAPSED_TIME});
p($dta{ENSEMBLE}[$e]->{SECNO});
pb($dta{ENSEMBLE}[$e]->{UNIX_TIME} < $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME});
p($dta{ENSEMBLE}[$e]->{W});
p($dta{ENSEMBLE}[$e]->{W_ERR});
p($dta{ENSEMBLE}[$e]->{ERR_VEL});
p($dta{ENSEMBLE}[$e]->{DEPTH});
p($dta{ENSEMBLE}[$e]->{DEPTH_ERR});
p($dta{ENSEMBLE}[$e]->{seabed});
p(&gimbal_pitch($dta{ENSEMBLE}[$e]->{PITCH},$dta{ENSEMBLE}[$e]->{ROLL}));
p($dta{ENSEMBLE}[$e]->{ROLL});
p($dta{ENSEMBLE}[$e]->{TILT});
p($dta{ENSEMBLE}[$e]->{HEADING});
p($dta{ENSEMBLE}[$e]->{ROTATION});
if (defined($opt_M)) {
p($dta{ENSEMBLE}[$e]->{U}); p($dta{ENSEMBLE}[$e]->{U_ERR});
p($dta{ENSEMBLE}[$e]->{V}); p($dta{ENSEMBLE}[$e]->{V_ERR});
p($dta{ENSEMBLE}[$e]->{X}); p($dta{ENSEMBLE}[$e]->{X_ERR});
p($dta{ENSEMBLE}[$e]->{Y}); p($dta{ENSEMBLE}[$e]->{Y_ERR});
}
if (@f) {
foreach $f (@f) {
my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
$fn = $f unless defined($fn);
p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
}
}
print("\n");
}
exit(0);