--- a/HISTORY
+++ b/HISTORY
@@ -1,9 +1,9 @@
======================================================================
H I S T O R Y
doc: Tue May 15 18:04:39 2012
- dlm: Tue May 15 18:31:09 2012
+ dlm: Fri Jul 12 12:46:50 2013
(c) 2012 A.M. Thurnherr
- uE-Info: 10 18 NIL 0 0 72 3 2 8 NIL ofnI
+ uE-Info: 17 30 NIL 0 0 72 3 2 8 NIL ofnI
======================================================================
May 15, 2012:
@@ -11,3 +11,8 @@
- began history
- uploaded current version to server for use with first version
of re-implemented shear method
+
+Jul 11, 2013:
+ - V1.1 [.hg/hgrc]
+ - various minor improvements
+
--- a/RDI_BB_Read.pl
+++ b/RDI_BB_Read.pl
@@ -1,9 +1,9 @@
#======================================================================
# R D I _ B B _ R E A D . P L
# doc: Sat Jan 18 14:54:43 2003
-# dlm: Mon Mar 25 21:38:37 2013
+# dlm: Mon Apr 29 12:49:40 2013
# (c) 2003 A.M. Thurnherr
-# uE-Info: 632 0 NIL 0 0 72 74 2 4 NIL ofnI
+# uE-Info: 53 46 NIL 0 0 72 74 2 4 NIL ofnI
#======================================================================
# Read RDI BroadBand Binary Data Files (*.[0-9][0-9][0-9])
@@ -49,6 +49,8 @@
# Mar 19, 2013: - added support for WH600 data file (58 fixed leader bytes)
# Mar 20, 2013: - removed DATA_FORMAT stuff
# - added support for BT data in subset of ensembles
+# Apr 29, 2013: - changed semantics to assume EOF when unexpected number of data types
+# are present in an ensemble
# FIRMWARE VERSIONS:
# It appears that different firmware versions generate different file
@@ -532,7 +534,7 @@
($hid,$did,$ens_length,$dummy,$ndt) = unpack('CCvCC',$buf);
$hid == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Header",$hid,0));
$did == 0x7f || die(sprintf($FmtErr,$WBRcfn,"Data Source",$did,0));
- printf(STDERR "\n$WBRcfn: WARNING: unexpected number of data types (%d)\n",$ndt)
+ printf(STDERR "\n$WBRcfn: WARNING: unexpected number of data types (%d, ens=$ens)\n",$ndt),last
unless ($ndt == 6 || $ndt == 7);
$BT_present = ($ndt == 7);
read(WBRF,$buf,2*$ndt) == 2*$ndt || die("$WBRcfn: $!");
--- a/RDI_Utils.pl
+++ b/RDI_Utils.pl
@@ -1,9 +1,9 @@
#======================================================================
# R D I _ U T I L S . P L
# doc: Wed Feb 12 10:21:32 2003
-# dlm: Fri Apr 12 09:22:10 2013
+# dlm: Thu Jun 20 15:26:47 2013
# (c) 2003 A.M. Thurnherr
-# uE-Info: 44 68 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 125 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# miscellaneous RDI-specific utilities
@@ -42,6 +42,9 @@
# Mar 27, 2013: - BUG: 3-beam solutions were not used in ref_lr_w
# - disabled apparently unused code
# Apr 12, 2013: - added $min_pctg as optional parameter to mk_prof
+# May 14, 2013: - added incident-velocity, w12 & w34 to mkProfile
+# Jun 5, 2013: - BUG: incident-flow warning was printed repeatedly
+# Jun 20, 2013: - BUG: warning had used &antsInfo()
use strict;
@@ -218,10 +221,15 @@
# mk_prof($dta,$check,$filter,$lr_b0,$lr_b1,$min_corr,$max_e,$max_gap);
#======================================================================
-sub ref_lr_w($$$$$$$) # calc ref-level vert vels
+# calculate reference-layer vertical and incident velocities
+
+{ my($warned); # static scope
+
+sub ref_lr_w($$$$$$$)
{
my($dta,$ens,$rl_b0,$rl_b1,$min_corr,$max_e,$min_pctg) = @_;
- my($i,@n,@bn,@v,@vel,@bv,@w);
+ my($i,@n,@bn,@v,@vi,@vel,@veli,@bv,@w);
+ my($w,$e,$nvi,$vi12,$vi43,@vbp,@velbp,@nbp,$w12,$w34,@w12,@w34);
for ($i=$rl_b0; $i<=$rl_b1; $i++) {
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][0])
@@ -241,26 +249,44 @@
if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < $min_pctg);
undef($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3])
if ($dta->{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < $min_pctg);
- @v = velInstrumentToEarth($dta,$ens,
- velBeamToInstrument($dta,
- @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
+ @vi = velBeamToInstrument($dta,@{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]});
+ @v = velInstrumentToEarth($dta,$ens,@vi);
+ @vbp = velBeamToBPEarth($dta,$ens,@{$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] < $min_pctg);
@v = @{$dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i]};
- # NB: no need to apply heading bias, as long as we only use w!
+ unless ($warned) {
+ print(STDERR "WARNING: incident-flow & beam-pair velocities not yet implemented for earth coordinates");
+ $warned = 1;
+ }
}
### next if (!defined($v[3]) || abs($v[3]) > $max_e); # disallow 3-beam solutions
next if (defined($v[3]) && abs($v[3]) > $max_e); # allow 3-beam solutions
- if (defined($v[2])) { # valid w
- $vel[2] += $v[2]; $n[2]++;
- $vel[3] += $v[3], $n[3]++ if defined($v[3]);
- push(@w,$v[2]); # for stderr test
+ if (defined($v[2])) { # valid velocity
+ $vel[2] += $v[2]; $n[2]++; # vertical velocity
+ $vel[3] += $v[3], $n[3]++ if defined($v[3]); # error velocity
+ push(@w,$v[2]); # save for stderr calculation
+ }
+
+ if (defined($vbp[1])) { # beam-pair vertical velocities
+ $velbp[0] += $vbp[1]; $nbp[0]++;
+ push(@w12,$vbp[1]);
+ }
+ if (defined($vbp[3])) {
+ $velbp[1] += $vbp[3]; $nbp[1]++;
+ push(@w34,$vbp[1]);
}
+ if (defined($vi[0])) { # incident velocity
+ $veli[0] += $vi[0];
+ $veli[1] += $vi[1];
+ $nvi++;
+ }
+
# The following code calculates beam-averaged ref-lr velocities.
# I do not recall what this was implemented for. Disabled Mar 27, 2013.
#
@@ -274,22 +300,63 @@
# $bv[3] += $dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3], $bn[3]++
# if defined($dta->{ENSEMBLE}[$ens]->{VELOCITY}[$i][3]);
# }
+ } # loop over ref-lr bins
+
+ $w = ($n[2] > 0) ? $vel[2]/$n[2] : undef; # calc means
+ $e = ($n[3] > 0) ? $vel[3]/$n[3] : undef;
+ if ($nvi > 0) {
+ $vi12 = $veli[0] / $nvi;
+ $vi43 = $veli[1] / $nvi;
+ } else {
+ $vi12 = $vi43 = undef;
+ }
+ $w12 = ($nbp[0] > 0) ? $velbp[0]/$nbp[0] : undef;
+ $w34 = ($nbp[1] > 0) ? $velbp[1]/$nbp[1] : undef;
+
+ if (@w12) { # w uncertainty
+ my($sumsq) = 0;
+ for ($i=0; $i<=$#w12; $i++) {
+ $sumsq += ($w12-$w12[$i])**2;
+ }
+ $dta->{ENSEMBLE}[$ens]->{W12} = $w12;
+ $dta->{ENSEMBLE}[$ens]->{W12_ERR} = sqrt($sumsq)/($nbp[0]-1)
+ if ($nbp[0]>=2);
}
- my($w) = $n[2] ? $vel[2]/$n[2] : undef; # w uncertainty
- my($sumsq) = 0;
+ if (@w34) { # w uncertainty
+ my($sumsq) = 0;
+ for ($i=0; $i<=$#w34; $i++) {
+ $sumsq += ($w34-$w34[$i])**2;
+ }
+ $dta->{ENSEMBLE}[$ens]->{W34} = $w34;
+ $dta->{ENSEMBLE}[$ens]->{W34_ERR} = sqrt($sumsq)/($nbp[1]-1)
+ if ($nbp[1]>=2);
+ }
+
+ my($sumsq) = 0; # w uncertainty
for ($i=0; $i<=$#w; $i++) {
$sumsq += ($w-$w[$i])**2;
}
my($stderr) = $n[2]>=2 ? sqrt($sumsq)/($n[2]-1) : undef;
+
# The following stderr test introduces a huge gap near the bottom of
# the profiles. Without it, they seem more reasonable.
# next if ($stderr > 0.05);
- if (defined($w)) { # valid w
+ if (defined($w)) { # valid velocity
$dta->{ENSEMBLE}[$ens]->{W} = $w;
$dta->{ENSEMBLE}[$ens]->{W_ERR} = $stderr;
}
+ $dta->{ENSEMBLE}[$ens]->{ERR_VEL} = $e if (defined($e));
+
+ $dta->{ENSEMBLE}[$ens]->{W12} = $w12 if (defined($w12));
+ $dta->{ENSEMBLE}[$ens]->{W34} = $w34 if (defined($w34));
+
+ if (defined($vi12)) {
+ $dta->{ENSEMBLE}[$ens]->{INCIDENT_VEL_T12} = $vi12;
+ $dta->{ENSEMBLE}[$ens]->{INCIDENT_VEL_T43} = $vi43;
+ }
+
# The following code calculates beam-averaged ref-lr velocities.
# I do not recall what this was implemented for. Disabled Mar 27, 2013.
#
@@ -299,8 +366,11 @@
# $dta->{ENSEMBLE}[$ens]->{V3} = $bn[2]>=2 ? $bv[2]/$bn[2] : undef;
# $dta->{ENSEMBLE}[$ens]->{V4} = $bn[3]>=2 ? $bv[3]/$bn[3] : undef;
# }
+
}
+} # static scope
+
sub mk_prof(...) # make profile
{
--- a/checkIX
+++ b/checkIX
@@ -2,37 +2,47 @@
#======================================================================
# C H E C K I X
# doc: Wed Dec 12 15:58:56 2012
-# dlm: Wed Dec 12 16:27:11 2012
+# dlm: Mon Apr 22 15:17:17 2013
# (c) 2012 A.M. Thurnherr
-# uE-Info: 31 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 43 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Dec 12, 2012: - created
-die("Usage: $0 <stn>\n")
- unless (@ARGV == 1);
+die("Usage: $0 <count-spec[ ...]>\n")
+ unless (@ARGV > 0);
+
+foreach my $id (`count @ARGV`) {
-$id = $ARGV[0]; # determine station id
-$id = sprintf('%03d',$id)
- unless (-f "$id.lad");
+ $id = sprintf('%03d',$id)
+ unless (-f "$id.log");
+
+ unless (-f "$id.log") {
+ print(STDERR "$id: missing station\n");
+ next;
+ }
-die("file <$id.lad> missing\n") # ensure required output is here
- unless (-f "$id.lad");
-die("file <$id.log> missing\n")
- unless (-f "$id.log");
-die("file <$id.mat> missing\n")
- unless (-f "$id.mat");
-die("file <$id.txt> missing\n")
- unless (-f "$id.txt");
+ die("$id: file <$id.lad> missing\n") # ensure required output is here
+ unless (-f "$id.lad");
+ die("$id: file <$id.log> missing\n")
+ unless (-f "$id.log");
+ die("$id: file <$id.mat> missing\n")
+ unless (-f "$id.mat");
+ die("$id: file <$id.txt> missing\n")
+ unless (-f "$id.txt");
+
+ if (-f "${id}_11.ps" && # handle warnings figure
+ length(`grep 'LADCP profile OK' ${id}_11.ps`) == 0) {
+ print("$id: warnings produced\n");
+ system("gv ${id}_11.ps &");
+ }
+
+ print("$id: no valid BT data\n") # check validity of ancillary data
+ unless (-f "$id.bot");
+ print("$id: no valid SADCP data\n")
+ unless (length(`grep 'all SADCP values removed' $id.log`) == 0);
-if (-f "${id}_11.ps" && # handle warnings figure
- length(`grep 'LADCP profile OK' ${id}_11.ps`) == 0) {
- print("warnings produced\n");
- system("gv ${id}_11.ps &");
}
-print("no valid BT data\n") # check validity of ancillary data
- unless (-f "$id.bot");
-print("no valid SADCP data\n")
- unless (length(`grep 'all SADCP values removed' $id.log`) == 0);
+exit(0);
--- a/listBins
+++ b/listBins
@@ -2,9 +2,9 @@
#======================================================================
# L I S T B I N S
# doc: Fri Aug 25 15:57:05 2006
-# dlm: Sun Aug 22 22:40:32 2010
+# dlm: Mon Apr 29 11:50:06 2013
# (c) 2006 A.M. Thurnherr
-# uE-Info: 45 28 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 108 20 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# Split data file into per-bin time series.
@@ -43,6 +43,8 @@
# May 22, 2009: - added -B) <bias/bias/bias/bias>
# May 23, 2009: - adapted to changed beampair-velocity fun name
# Aug 22, 2010: - added -R
+# Apr 29, 2013: - cosmetics
+# - added warning on missing -S
# General Notes:
# - everything (e.g. beams) is numbered from 1
@@ -76,7 +78,7 @@
"[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
"[-P)itch/Roll <bias/bias>] [-B)eamvel <bias/bias/bias/bias>] " .
"[require -4)-beam solutions] [-d)iscard <beam#>] " .
- "[-%)good <min>] " .
+ "[-p)ct-good <min>] " .
"[output -b)eam coordinates] [output two separate -w) estimates] " .
"<RDI file>\n")
unless (&Getopts("4abB:d:f:M:p:r:P:RS:w") && @ARGV == 1);
@@ -103,9 +105,12 @@
($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 '*');
+if (defined($opt_S)) {
+ ($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S);
+ $variable_ssCorr = ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*');
+} else {
+ print(STDERR "WARNING: no soundspeed correction applied!\n");
+}
#----------------------------------------------------------------------
--- a/mkProfile
+++ b/mkProfile
@@ -2,9 +2,9 @@
#======================================================================
# M K P R O F I L E
# doc: Sun Jan 19 18:55:26 2003
-# dlm: Fri Apr 12 09:24:14 2013
+# dlm: Tue May 14 11:37:04 2013
# (c) 2003 A.M. Thurnherr
-# uE-Info: 269 72 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 169 15 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# Make an LADCP Profile by Integrating W (similar to Firing's scan*).
@@ -79,6 +79,10 @@
# 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
# NOTES:
# - the battery values are based on transmission voltages (different
@@ -104,15 +108,14 @@
die("Usage: $0 " .
"[-Q)uiet] [-F)ilter <script>] " .
"[-s)uppress checkensemble()] " .
- "[require -4)-beam solutions] " .
+ "[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]>] " .
- "[-m)ax <gap>] " .
- "[-d)rift <dx,dy>] [-g)ps <start lat,lon/end lat,lon>] " .
+ "[max -g)ap <len>] " .
"[output -f)ields <field[,...]> " .
"[-M)agnetic <declination>] [profile -B)ottom <depth>] " .
"<RDI file>\n")
- unless (&Getopts("4AB:F:M:Qd:r:n:e:c:g:f:m:sp:") && @ARGV == 1);
+ unless (&Getopts("4AB:F:M:Qd:g:r:n:e:c:f:m:sp:") && @ARGV == 1);
$RDI_Coords::minValidVels = 4 if ($opt_4); # no 3-beam solutions
@@ -122,25 +125,12 @@
$opt_n = 2 unless defined($opt_n);
$opt_e = 0.1 unless defined($opt_e);
$opt_c = 70 unless defined($opt_c);
-$opt_m = 120 unless defined($opt_m);
+$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);
-if ($opt_g) { # GPS info
- ($s_lat,$s_lon,$e_lat,$e_lon) = gps_to_deg($opt_g);
- $lat = $s_lat/2 + $e_lat/2;
- $lon = $s_lon/2 + $e_lon/2;
- $ddx = dist($lat,$s_lon,$lat,$e_lon);
- $ddy = dist($s_lat,$lon,$e_lat,$lon);
-}
-
-if ($opt_d) { # ship drift
- ($ddx,$ddy) = split(',',$opt_d);
- die("$0: can't decode -d $opt_d\n") unless defined($ddy);
-}
-
print(STDERR "Reading $ARGV[0]..."); # read data
readData($ARGV[0],\%dta);
print(STDERR "done\n");
@@ -153,6 +143,45 @@
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 {
@@ -178,56 +207,6 @@
}
#======================================================================
-# Misc funs used to decode options
-#======================================================================
-
-sub dist($$$$) # distance
-{
- my($lat1,$lon1,$lat2,$lon2) = @_;
- my($a) = 6378139; # Earth's radius
-
- $lat1 = rad($lat1); $lon1 = rad($lon1);
- $lat2 = rad($lat2); $lon2 = rad($lon2);
- $ct1 = cos($lat1); $st1 = sin($lat1); $cp1 = cos($lon1);
- $sp1 = sin($lon1); $ct2 = cos($lat2); $st2 = sin($lat2);
- $cp2 = cos($lon2); $sp2 = sin($lon2);
- $cos = $ct1*$cp1*$ct2*$cp2 + $ct1*$sp1*$ct2*$sp2 + $st1*$st2;
- $cos = 1 if ($cos > 1);
- $cos = -1 if ($cos < -1);
- return $a * acos($cos);
-}
-
-sub deg_to_dec($) # parse degrees
-{
- my($deg,$min) = split(':',$_[0]);
- return $deg + $min/60;
-}
-
-sub gps_to_deg($) # decode lat/lon
-{
- my($start,$end) = split('/',$_[0]);
- my($sa,$so,$ea,$eo);
-
- my($lat,$lon) = split(',',$start);
- if ($lat =~ m{N$}) { $sa = deg_to_dec($`); }
- elsif ($lat =~ m{S$}) { $sa = -deg_to_dec($`); }
- else { $sa = $lat; }
- if ($lon =~ m{E$}) { $so = deg_to_dec($`); }
- elsif ($lon =~ m{W$}) { $so = -deg_to_dec($`); }
- else { $so = $lon; }
-
- my($lat,$lon) = split(',',$end);
- if ($lat =~ m{N$}) { $ea = deg_to_dec($`); }
- elsif ($lat =~ m{S$}) { $ea = -deg_to_dec($`); }
- else { $ea = $lat; }
- if ($lon =~ m{E$}) { $eo = deg_to_dec($`); }
- elsif ($lon =~ m{W$}) { $eo = -deg_to_dec($`); }
- else { $eo = $lon; }
-
- return ($sa,$so,$ea,$eo);
-}
-
-#======================================================================
# Step 0: Check data & Calculate Ping Rates
#======================================================================
@@ -266,7 +245,7 @@
#======================================================================
($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz,$rms_heave_accel) =
- mk_prof(\%dta,!$opt_s,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_m,$opt_p);
+ mk_prof(\%dta,!$opt_s,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_g,$opt_p);
unless (($atbottom > $firstgood) && ($lastgood > $atbottom)) {
if ($opt_Q) {
@@ -575,27 +554,9 @@
$gb = ($gb[0]+$gb[1]+$gb[2]+$gb[3]) / 4;
#======================================================================
-# Step 5: Remove Ship Drift (probably not useful)
+# Step 5: Remove Ship Drift (probably not useful => removed)
#======================================================================
-if (defined($opt_M) && defined($ddx)) { # remove barotropic
- $du = $ddx / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};# mean drift vel
- $dv = $ddy / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
- $iu = $dta{ENSEMBLE}[$lastgood]->{X} / # mean obs vel
- $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
- $iv = $dta{ENSEMBLE}[$lastgood]->{Y} /
- $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
-
- for ($e=$firstgood; $e<=$lastgood; $e++) {
- next unless (defined($dta{ENSEMBLE}[$e]->{X}) &&
- defined($dta{ENSEMBLE}[$e]->{Y}));
- $dta{ENSEMBLE}[$e]->{U} -= $du;
- $dta{ENSEMBLE}[$e]->{V} -= $dv;
- $dta{ENSEMBLE}[$e]->{X} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($du-$iu);
- $dta{ENSEMBLE}[$e]->{Y} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($dv-$iv);
- }
-}
-
#======================================================================
# Step 6: Pitch, Roll, Rotation
#======================================================================
@@ -652,7 +613,13 @@
$prrms = sqrt($prrms/($lastgood-$firstgood));
$dnprrms = sqrt($dnprrms/($atbottom-$firstgood));
-$upprrms = sqrt($upprrms/($lastgood-$atbottom));
+
+if ($lastgood == $atbottom) {
+ print(STDERR "WARNING: $0 NO UPCAST DATA\n");
+ $upprrms = nan;
+} else {
+ $upprrms = sqrt($upprrms/($lastgood-$atbottom));
+}
$rotrms = sqrt($rotrms/($lastgood-$firstgood));
@@ -718,8 +685,8 @@
$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} {depth} {depth_err} {depth_BT} " .
- "{pitchroll} {rotation} " .
+ "{w} {w_err} {err_vel} {depth} {depth_err} {depth_BT} " .
+ "{pitchroll} {heading} {rotation} " .
"$uFields $addFields\n");
printf("#ANTS#PARAMS# date{$dta{ENSEMBLE}[$firstgood]->{DATE}} " .
@@ -784,19 +751,12 @@
$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
$vErr, $dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
) if defined ($opt_M);
-print("#ANTS#PARAMS# start_lat{$s_lat} start_lon{$s_lon} " .
- "end_lat{$e_lat} end_lon{$e_lon} " .
- "lat{$lat} lon{$lon}\n")
- if defined($lat);
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});
}
-printf("#ANTS#PARAMS# drift_x{%d} drift_y{%d} " .
- "drift_u{%.3f} drift_v{%.3f} " .
- "\n",$ddx,$ddy,$du,$dv) if defined($ddx);
if (defined($water_depth)) {
printf("#ANTS#PARAMS# water_depth{%d} sig-water_depth{%d}\n",
$water_depth,$sig_wd);
@@ -815,10 +775,12 @@
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]->{DEPTH_BT});
p($dta{ENSEMBLE}[$e]->{PITCHROLL});
+ 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});