.
new file mode 100755
--- /dev/null
+++ b/2011_07_10/LADCP2CTDmatch
@@ -0,0 +1,333 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P 2 C T D M A T C H
+# doc: Thu Sep 23 22:01:19 2010
+# dlm: Tue Sep 28 21:43:29 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 312 33 NIL 0 0 70 10 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'match LADCP to CTD time series';
+
+# HISTORY:
+# Sep 23, 2010: - incepted
+
+($ANTS) = ($0 =~ m{^(.*)/[^/]*$});
+$PERL_TOOLS = '/Data/LADCP/Software/Thurnherr/perl-tools';
+
+require "$ANTS/ants.pl";
+require "$PERL_TOOLS/RDI_BB_Read.pl";
+require "$PERL_TOOLS/RDI_Coords.pl";
+require "$PERL_TOOLS/RDI_Utils.pl";
+
+$antsParseHeader = 0;
+&antsUsage('dw:n:',2,
+ '[-d)ebug] ',
+ '[-w)indow <size[30s]>] [-n) <windows[20]]',
+ '<RDI file> <SeaBird file>');
+
+&antsCardOpt(\$opt_w,30);
+&antsCardOpt(\$opt_n,20);
+
+$LADCP_file = &antsFileArg();
+$CTD_file = &antsFileArg();
+
+&antsAddParams('LADCP_file',$LADCP_file,'CTD_file',$CTD_file);
+@antsNewLayout = ('elapsed','CTD_w','LADCP_w');
+&antsActivateOut();
+
+#----------------------------------------------------------------------
+# Step 1: Read Data
+#----------------------------------------------------------------------
+
+print(STDERR "Reading LADCP data ($LADCP_file)...");
+readData($LADCP_file,\%LADCP);
+printf(STDERR "\n\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));
+
+print(STDERR "Reading CTD data ($CTD_file)...");
+open(F,$CTD_file) || croak("$CTD_file: $!\n");
+while (1) { # parse header
+ chomp($hdr = <F>);
+ $hdr =~ s/\r*$//;
+ croak("$0: unexpected EOF (format error)\n") unless defined($hdr);
+ last if ($hdr eq '*END*');
+
+ $CTD_nfields = $',next if ($hdr =~ /nquan = /); # Layout
+ $CTD_nrecs = $',next if ($hdr =~ /nvalues = /);
+ $pressF = $1,next if ($hdr =~ /name (\d+) = prDM:/);
+
+ &antsAddParams('start_time',$1),next # selected metadata
+ if ($hdr =~ /start_time = (.*)/);
+
+ &antsAddParams('station',$1),next
+ if ($hdr =~ /Station\s*:\s*(.*)/);
+ &antsAddParams('ship',$1),next
+ if ($hdr =~ /Ship\s*:\s*(.*)\s*$/);
+ &antsAddParams('cruise',$1),next
+ if ($hdr =~ /Cruise\s*:\s*(.*)\s*$/);
+ &antsAddParams('time',$1),next
+ if ($hdr =~ /Time\s*:\s*(.*)/);
+ &antsAddParams('date',$1),next
+ if ($hdr =~ /Date\s*:\s*(.*)/);
+
+ if ($hdr =~ /Latitude\s*[=:]\s*/) {
+ ($deg,$min,$NS) = split(/ /,$');
+ $lat = $deg + $min/60;
+ $lat *= -1 if ($NS eq 'S');
+ &antsAddParams('lat',$lat);
+ next;
+ }
+ if ($hdr =~ /Longitude\s*[=:]\s*/) {
+ ($deg,$min,$EW) = split(/ /,$');
+ $lon = $deg + $min/60;
+ $lon *= -1 if ($EW eq 'W');
+ &antsAddParams('lon',$lon);
+ next;
+ }
+
+ if ($hdr =~ /interval = seconds: /) {
+ $CTD_sampint = 1*$';
+ &antsAddParams('CTD_interval',1/$CTD_sampint);
+ next;
+ }
+
+ $CTD_badval = $',next
+ if ($hdr =~ /bad_flag = /);
+ $CTD_file_type = $',next
+ if ($hdr =~ /file_type = /);
+}
+
+croak("$CTD_file: cannot determine CTD file layout\n")
+ unless ($CTD_nfields && $CTD_nrecs);
+croak("$CTD_file: cannot determine missing value\n")
+ unless defined($CTD_badval);
+croak("$CTD_file: not a CTD time series file\n")
+ unless ($CTD_sampint);
+croak("$CTD_file: no pressure field\n")
+ unless defined($pressF);
+
+croak("$0: unknown latitude\n") unless defined($lat);
+
+if ($CTD_file_type eq 'ascii') {
+ while (1) {
+ last unless (@rec = &antsFileIn(F));
+ push(@CTD_press,($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
+ }
+} elsif ($CTD_file_type eq 'binary') {
+
+ my($fbits) = 8 * length(pack('f',0));
+ croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
+ unless ($fbits == 32);
+
+ croak("$CTD_file: can't read binary data\n")
+ unless (read(F,$dta,4*$CTD_nfields*$CTD_nrecs) == 4*$CTD_nfields*$CTD_nrecs);
+ print(STDERR "$CTD_file: WARNING: extraneous data at EOF\n") unless eof(F);
+
+ $dta = pack('V*',unpack('N*',$dta)) # big-endian CPU
+ if (unpack('h*', pack('s', 1)) =~ /01/); # c.f. perlport(1)
+
+ @dta = unpack("f*",$dta);
+
+ for ($r=0; $r<$CTD_nrecs; $r++) {
+ push(@CTD_press,($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
+ }
+} else {
+ croak("$CTD_file: unknown CTD file type $CTD_file_type\n");
+}
+
+printf(STDERR "\n\t%d scans\n",scalar(@CTD_press));
+
+#----------------------------------------------------------------------
+# Step 2a: Pre-Process CTD Data
+#----------------------------------------------------------------------
+
+printf(STDERR "Pre-processing data...");
+printf(STDERR "\n\tCTD...");
+
+#------------------------------------
+# calculate w and find deepest record
+#------------------------------------
+$CTD_maxpress = -9e99;
+for (my($r)=1; $r<$CTD_nrecs-1; $r++) {
+ print(STDERR '.') if ($r%500==0);
+ $CTD_w[$r] = 0.99*($CTD_press[$r+1] - $CTD_press[$r-1]) / (2*$CTD_sampint);
+ if ($CTD_press[$r] > $CTD_maxpress) {
+ $CTD_maxpress = $CTD_press[$r];
+ $CTD_bottom = $r;
+ }
+}
+
+#----------------------------------------------------------------------
+# Step 2b: Pre-Process LADCP Data
+#----------------------------------------------------------------------
+
+print(STDERR "\n\tLADCP...");
+
+#------------------------------------------------------
+# construct a depth-vs-time "profile" from the raw data
+#------------------------------------------------------
+
+my($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
+ mk_prof(\%LADCP,0,undef,1,6,70,0.1,120);
+croak("\n$LADCP_file: no good ensembles found\n")
+ unless defined($LADCP_start);
+
+if ($opt_d) {
+ printf(STDERR "\n\t\tStart of cast : %s (#%5d) at %6.1fm\n",
+ $LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_start]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH});
+ printf(STDERR "\t\tBottom of cast : %s (#%5d) at %6.1fm\n",
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+ if (defined($water_depth)) {
+ printf(STDERR "\t\tSeabed : at %6.1fm (+-%dm)\n",$water_depth,$sig_wd);
+ } else {
+ print(STDERR "\t\tSeabed : not found\n");
+ }
+ printf(STDERR "\t\tEnd of cast : %s (#%5d) at %6.1fm\n",
+ $LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_end]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
+}
+
+#------------------------------------------------------------------------
+# find 1st rec & ensemble >=80% down to max depth & make 1st guess at lag
+#------------------------------------------------------------------------
+
+my($CTD_80pct_down) = 0;
+$CTD_80pct_down++
+ while ($CTD_press[$CTD_80pct_down] < 0.8*$CTD_maxpress);
+
+my($LADCP_80pct_down) = 0;
+$LADCP_80pct_down++
+ while ($LADCP{ENSEMBLE}[$LADCP_80pct_down]->{DEPTH} < 0.8*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+
+
+my($first_guess_lag) = $LADCP{ENSEMBLE}[$LADCP_80pct_down]->{ELAPSED_TIME} -
+ $CTD_80pct_down*$CTD_sampint;
+
+printf(STDERR "\t\t1st guess offset = %ds\n",$first_guess_lag*$CTD_sampint)
+ if ($opt_d);
+
+#------------------------------------------------------------------------------------
+# Linearly interpolate LADCP time series onto a new grid with $CTD_sampint resolution
+# ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data
+# chunks large
+#------------------------------------------------------------------------------------
+
+my($nGaps) = 0;
+
+sub interp_LADCP_w($$)
+{
+ my($elapsed,$ens) = @_;
+ my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) /
+ ($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} -
+ $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME});
+ unless (numberp($LADCP{ENSEMBLE}[$ens-1]->{W})) {
+ $nGaps++;
+ croak("$0: cannot (yet) handle double gaps --- request feature upgrade\n")
+ unless numberp($LADCP{ENSEMBLE}[$ens]->{W});
+ return $LADCP{ENSEMBLE}[$ens]->{W};
+ }
+ return $LADCP{ENSEMBLE}[$ens-1]->{W} +
+ $sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W});
+}
+
+for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {
+ while ($r*$CTD_sampint < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
+ $LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD_sampint,$ens)
+ unless ($first_guess_lag > $r);
+ $r++;
+ }
+}
+
+print(STDERR "\t\t$nGaps gaps in w timeseries")
+ if ($opt_d);
+
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 3: Calculate Lags
+#----------------------------------------------------------------------
+
+printf(STDERR "Calculating $opt_n lags from %ds-long windows [s]: ",$opt_w);
+$opt_w = int($opt_w / $CTD_sampint);
+
+sub bestLag($)
+{
+ my($ws) = @_; # window start index
+
+ my($best);
+ my($bestmad) = 9e99; # mean absolute deviation
+ for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) {
+ my($sad) = my($nad) = 0;
+ for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
+ my($Li) = $Ci + $Llag;
+ next if ($Li<0 || $Li>=$opt_w);
+ next unless numberp($CTD_w[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
+ $sad += abs($CTD_w[$ws+$Ci] - $LADCP_w[$ws+$Li]);
+ $nad++;
+ }
+ next unless ($nad > 0);
+ if ($sad/$nad < $bestmad) {
+ $best = $Llag;
+ $bestmad = $sad/$nad;
+ }
+ }
+ return $best;
+}
+
+#---------------------------------------------------------------
+# carry out opt_n lag correlations and keep tally of the results
+#---------------------------------------------------------------
+my(@nBest);
+for (my($window)=0; $window<$opt_n; $window++) {
+ my($ws) = $window * int(@LADCP_w/$opt_n); # window start
+ $ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);
+
+ $bestLag = bestLag($ws);
+ if (defined($bestLag)) {
+ printf(STDERR "%d ",$bestLag*$CTD_sampint);
+ $nBest{$bestLag}++;
+ } else {
+ printf(STDERR "nan ");
+ }
+}
+
+#----------------------
+# find most popular lag
+#----------------------
+my($best_lag);
+foreach my $i (keys(%nBest)) {
+ $best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
+}
+croak("\n$0: cannot determine a valid lag\n")
+ unless ($nBest{$best_lag} > 1);
+print(STDERR "\n\n\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
+ if ($nBest{$best_lag} < $opt_n/2);
+
+if ($nBest{$best_lag} == $opt_n) {
+ printf(STDERR "\n\tunanimous lag = %ds\n",$best_lag*$CTD_sampint);
+} else {
+ printf(STDERR "\n\tmost popular lag = %ds\n",$best_lag*$CTD_sampint);
+}
+
+#----------------------------------------------------------------------
+# Step 4: associate CTD data with each LADCP ensemble & write time series
+#----------------------------------------------------------------------
+
+print(STDERR "Associating CTD data with LADCP ensembles...\n");
+
+my($overall_dt) = $CTD_sampint * ($best_lag + $first_guess_lag);
+for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - $overall_dt) / $CTD_sampint);
+ $LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD_w[$r];
+ &antsOut($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{CTD_W},
+ $LADCP{ENSEMBLE}[$ens]->{W});
+}
+
+&antsExit();
+
new file mode 100755
--- /dev/null
+++ b/2011_07_10/LADCPintsh
@@ -0,0 +1,446 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P I N T S H
+# doc: Thu Oct 14 21:22:50 2010
+# dlm: Thu Jul 7 06:30:50 2011
+# (c) 2010 A.M. Thurnherr & E. Firing
+# uE-Info: 32 28 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'integrate LADCP shear';
+
+# NOTES:
+# - THE CORE OF THIS CODE IS A SIMPLE TRANSCRIPTION OF AVG_SH.M AND INT_SH.M WRITTEN BY ERIC FIRING
+# - comments beginning with ## are taken from Eric's code
+# - cubic velocity interpolation across PPI gap from Eric's code has not been implemented (yet?)
+# - Eric's ndata criteria had to be scaled by 5/8 due to the different gridding used; I decided,
+# conservatively, to scale only the -s criterion (n1 in Eric's code)
+
+# WEIRDNESSES:
+# - in Eric's [avg_sh.m] the calculation of output shear stddev incorrectly assumes that the 4th column
+# in the shear profile is stddev, rather than variance. However, as far as I can tell, this output
+# is not used anywhere in Eric's code
+
+# HISTORY:
+# Oct 14, 2010: - created
+# Oct 20, 2010: - first working version
+# Oct 23, 2010: - added support for -b)
+# Oct 24, 2010: - fix spuriously small variances that can occur for BT velocities based on very small
+# samples (i.e. primarily when chosing a small -r)
+# Dec 9, 2010: - allowed for empty BT file
+# Feb 16, 2011: - BUG: gaps in shear data were not handled correctly in baroclinic solution
+# Jul 7, 2011: - added -m
+
+($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
+require "$ANTS/ants.pl";
+require "$ANTS/libstats.pl";
+
+&antsUsage('b:dm:r:s:w:u:',0,
+ '[-d)ebug]',
+ '[reference with -b)ottom-track <file> [-m)in BT samp[10]]]',
+ '[-u)plooker <file>]',
+ '[min -r)aw <shear samp[10]>] [min -s)moothed <shear samp[20]>]',
+ '[-w)eighted avg <below[1200]>]',
+ '[LADCP shear file]');
+
+croak("$0: -m requires -b\n")
+ if defined($opt_m) && !defined($opt_b);
+
+&antsCardOpt(\$opt_r,10); ## minimum number of samples for shear
+&antsCardOpt(\$opt_m,10); # minimum number of samples for BT data
+&antsCardOpt(\$opt_s,20); ## minimum number of samples for smoothed shear
+&antsCardOpt(\$opt_w,1200); ## probably 1200 for midlatitude-equatorial, 2400 for high lat
+
+&antsFileOpt($opt_b); # BT file
+
+&antsFileOpt($opt_u); # UL shear file
+open(ULF,$opt_u) || croak("$opt_u: $!\n")
+ if defined($opt_u);
+
+#======================================================================
+# Step 1: Read and Average Shear Data
+# - depth bins with less than $opt_r values are blanked out
+#======================================================================
+
+$depthF = fnr('depth'); # layout of [LADCPproc] output
+$dc_nshF = fnr('dc_nshear');
+$dc_uzF = fnr('dc_u_z');
+$dc_uzsF = fnrNoErr('dc_u_z.sig');
+$dc_uzsF = fnr('dc_u_z_sig')
+ unless defined($dc_uzsF);
+$dc_vzF = fnr('dc_v_z');
+$dc_vzsF = fnrNoErr('dc_v_z.sig');
+$dc_vzsF = fnr('dc_v_z_sig')
+ unless defined($dc_vzsF);
+$dc_wzF = fnr('dc_w_z');
+$dc_wzsF = fnrNoErr('dc_w_z.sig');
+$dc_wzsF = fnr('dc_w_z_sig')
+ unless defined($dc_wzsF);
+$uc_nshF = fnr('uc_nshear');
+$uc_uzF = fnr('uc_u_z');
+$uc_uzsF = fnrNoErr('uc_u_z.sig');
+$uc_uzsF = fnr('uc_u_z_sig')
+ unless defined($uc_uzsF);
+$uc_vzF = fnr('uc_v_z');
+$uc_vzsF = fnrNoErr('uc_v_z.sig');
+$uc_vzsF = fnr('uc_v_z_sig')
+ unless defined($uc_vzsF);
+$uc_wzF = fnr('uc_w_z');
+$uc_wzsF = fnrNoErr('uc_w_z.sig');
+$uc_wzsF = fnr('uc_w_z_sig')
+ unless defined($uc_wzsF);
+
+my(@gaps); my($curGap) = 0;
+
+for (my($r)=0; &antsIn(); $r++) {
+ my(@UL_);
+ if (defined($opt_u)) {
+ @UL_ = &antsFileIn(ULF); # read UL shear data
+ undef($opt_u) unless (@UL_); # cheap trick
+ }
+
+ $depth[$r] = $ants_[0][$depthF]; ## depth grid values
+ croak("$opt_u: inconsistent depth (DL: $depth[$r]; UL: $UL_[$depthF])\n")
+ if defined($opt_u) && ($UL_[$depthF] != $depth[$r]);
+
+ $dc_ndata = $ants_[0][$dc_nshF]; # number of shear samples
+ $uc_ndata = $ants_[0][$uc_nshF];
+ if (defined($opt_u)) {
+ $dc_ndata += $UL_[$dc_nshF];
+ $uc_ndata += $UL_[$uc_nshF];
+ }
+ $ndata[$r] = $dc_ndata + $uc_ndata;
+
+ if (defined($opt_u)) {
+ if ($dc_ndata > 0) { # downcast shear
+ my($DLf) = $ants_[0][$dc_nshF] / $dc_ndata;
+ my($ULf) = $UL_[$dc_nshF] / $dc_ndata;
+ if ($DLf>0 && $Ulf>0) {
+ $dc_uz[$r] = $DLf*$ants_[0][$dc_uzF] + $ULf*$UL_[$dc_uzF];
+ $dc_vz[$r] = $DLf*$ants_[0][$dc_vzF] + $ULf*$UL_[$dc_vzF];
+ $dc_wz[$r] = $DLf*$ants_[0][$dc_wzF] + $ULf*$UL_[$dc_wzF];
+ } elsif ($DLf > 0) {
+ $dc_uz[$r] = $ants_[0][$dc_uzF];
+ $dc_vz[$r] = $ants_[0][$dc_vzF];
+ $dc_wz[$r] = $ants_[0][$dc_wzF];
+ } else {
+ $dc_uz[$r] = $UL_[$dc_uzF];
+ $dc_vz[$r] = $UL_[$dc_vzF];
+ $dc_wz[$r] = $UL_[$dc_wzF];
+ }
+ } else {
+ $dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = nan;
+ }
+ if ($uc_ndata > 0) { # upcast shear
+ my($DLf) = $ants_[0][$uc_nshF] / $uc_ndata;
+ my($ULf) = $UL_[$uc_nshF] / $uc_ndata;
+ if ($DLf>0 && $Ulf>0) {
+ $uc_uz[$r] = $DLf*$ants_[0][$uc_uzF] + $ULf*$UL_[$uc_uzF];
+ $uc_vz[$r] = $DLf*$ants_[0][$uc_vzF] + $ULf*$UL_[$uc_vzF];
+ $uc_wz[$r] = $DLf*$ants_[0][$uc_wzF] + $ULf*$UL_[$uc_wzF];
+ } elsif ($DLf > 0) {
+ $uc_uz[$r] = $ants_[0][$uc_uzF];
+ $uc_vz[$r] = $ants_[0][$uc_vzF];
+ $uc_wz[$r] = $ants_[0][$uc_wzF];
+ } else {
+ $uc_uz[$r] = $UL_[$uc_uzF];
+ $uc_vz[$r] = $UL_[$uc_vzF];
+ $uc_wz[$r] = $UL_[$uc_wzF];
+ }
+ } else {
+ $uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = nan;
+ }
+ } else {
+ if ($dc_ndata >= $opt_r) { # downcast shear
+ $dc_uz[$r] = $ants_[0][$dc_uzF];
+ $dc_vz[$r] = $ants_[0][$dc_vzF];
+ $dc_wz[$r] = $ants_[0][$dc_wzF];
+ } else {
+ $dc_uz[$r] = $dc_vz[$r] = $dc_wz[$r] = nan;
+ }
+ if ($uc_ndata >= $opt_r) { # upcast shear
+ $uc_uz[$r] = $ants_[0][$uc_uzF];
+ $uc_vz[$r] = $ants_[0][$uc_vzF];
+ $uc_wz[$r] = $ants_[0][$uc_wzF];
+ } else {
+ $uc_uz[$r] = $uc_vz[$r] = $uc_wz[$r] = nan;
+ }
+ }
+
+ if ($depth[$r] < $opt_w) { ## use average of up and down casts at top of profile
+ if ($dc_ndata>0 && $uc_ndata>0 && $ndata[$r]>=$opt_r) {
+ $uz[$r] = $dc_uz[$r]/2 + $uc_uz[$r]/2;
+ $vz[$r] = $dc_vz[$r]/2 + $uc_vz[$r]/2;
+ $wz[$r] = $dc_wz[$r]/2 + $uc_wz[$r]/2;
+ } elsif ($dc_ndata >= $opt_r) {
+ $uz[$r] = $dc_uz[$r];
+ $vz[$r] = $dc_vz[$r];
+ $wz[$r] = $dc_wz[$r];
+ } elsif ($uc_ndata >= $opt_r) {
+ $uz[$r] = $uc_uz[$r];
+ $vz[$r] = $uc_vz[$r];
+ $wz[$r] = $uc_wz[$r];
+ } else {
+ $uz[$r] = $vz[$r] = $wz[$r] = nan;
+ }
+ } else { ## use weighted average of up and down cast data at bottom of profile
+ if ($ndata[$r] >= $opt_r) {
+ my($dcf) = $dc_ndata / $ndata[$r];
+ my($ucf) = $uc_ndata / $ndata[$r];
+ if ($dcf>0 && $ucf>0) {
+ $uz[$r] = $dcf*$dc_uz[$r] + $ucf*$uc_uz[$r];
+ $vz[$r] = $dcf*$dc_vz[$r] + $ucf*$uc_vz[$r];
+ $wz[$r] = $dcf*$dc_wz[$r] + $ucf*$uc_wz[$r];
+ } elsif ($dcf > 0) {
+ $uz[$r] = $dc_uz[$r];
+ $vz[$r] = $dc_vz[$r];
+ $wz[$r] = $dc_wz[$r];
+ } else {
+ $uz[$r] = $uc_uz[$r];
+ $vz[$r] = $uc_vz[$r];
+ $wz[$r] = $uc_wz[$r];
+ }
+ } else {
+ $uz[$r] = $vz[$r] = $wz[$r] = nan;
+ }
+ }
+
+ if (numberp($uz[$r]) && $curGap>0) { # end of gap
+ push(@gaps,$curGap) unless ($r == $curGap); # do not report "gap" at beginning of profile
+# print(STDERR "$curGap-gap at $depth[$r]m\n");
+ $curGap = 0;
+ } elsif (!numberp($uz[$r])) { # currently in gap
+# print(STDERR "in gap at $depth[$r]m (ndata = $ndata[$r], $dc_ndata,$uc_ndata)\n");
+ $curGap++;
+ }
+
+ if ($ndata[$r] >= $opt_r) { ## back to sum of squares
+ if (defined($opt_u)) {
+ $uzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_uzsF]**2 +
+ $UL_[$dc_nshF] * $UL_[$dc_uzsF]**2 +
+ $ants_[0][$uc_nshF] * $ants_[0][$uc_uzsF]**2 +
+ $UL_[$uc_nshF] * $UL_[$uc_uzsF]**2) / $ndata[$r]);
+ $vzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_vzsF]**2 +
+ $UL_[$dc_nshF] * $UL_[$dc_vzsF]**2 +
+ $ants_[0][$uc_nshF] * $ants_[0][$uc_vzsF]**2 +
+ $UL_[$uc_nshF] * $UL_[$uc_vzsF]**2) / $ndata[$r]);
+ $wzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_wzsF]**2 +
+ $UL_[$dc_nshF] * $UL_[$dc_wzsF]**2 +
+ $ants_[0][$uc_nshF] * $ants_[0][$uc_wzsF]**2 +
+ $UL_[$uc_nshF] * $UL_[$uc_wzsF]**2) / $ndata[$r]);
+ } else {
+ $uzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_uzsF]**2 +
+ $ants_[0][$uc_nshF] * $ants_[0][$uc_uzsF]**2) / $ndata[$r]);
+ $vzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_vzsF]**2 +
+ $ants_[0][$uc_nshF] * $ants_[0][$uc_vzsF]**2) / $ndata[$r]);
+ $wzsig[$r] = sqrt(($ants_[0][$dc_nshF] * $ants_[0][$dc_wzsF]**2 +
+ $ants_[0][$uc_nshF] * $ants_[0][$uc_wzsF]**2) / $ndata[$r]);
+ }
+ } else {
+ $uzsig[$r] = $vzsig[$r] = $wzsig[$r] = nan;
+ }
+}
+
+if (@gaps) {
+ &antsAddParams('shear_gaps',"@gaps");
+ print(STDERR "shear gaps: @gaps\n");
+} else {
+ &antsAddParams('shear_gaps',0);
+}
+
+#===============================================================================
+# Step 2: Low-Pass filter high-quality shear data; not yet implemented
+#===============================================================================
+
+for (my($r)=0; $r<@depth; $r++) { # create high-quality shear profile
+ if ($ndata[$r] >= $opt_s) {
+ $hq_uz[$r] = $uz[$r]; $hq_vz[$r] = $vz[$r]; $hq_wz[$r] = $wz[$r];
+ } else {
+ $hq_uz[$r] = $hq_vz[$r] = $hq_wz[$r] = nan;
+ }
+}
+
+#@antsNewLayout = ('depth','u_z','hq_u_z','v_z','hq_v_z','w_z','hq_w_z');
+#for (my($r)=0; $r<@depth; $r++) {
+# &antsOut($depth[$r],$uz[$r],$hq_uz[$r],$vz[$r],$hq_vz[$r],$wz[$r],$hq_wz[$r]);
+#}
+
+#======================================================================
+# Step 3: Integrate Shear
+# - z(vel) = z(sh) + DZ/2
+#======================================================================
+
+my($DZ) = $depth[1] - $depth[0];
+
+for (my($r)=my($u)=my($v)=my($w)=my($dc_u)=my($dc_v)=my($dc_w)=my($uc_u)=my($uc_v)=my($uc_w)=0;
+ $r<@depth; $r++) {
+ $u = $u[$r] = $u + $DZ*$uz[$r] if numberp($uz[$r]);
+ $v = $v[$r] = $v + $DZ*$vz[$r] if numberp($vz[$r]);
+ $w = $w[$r] = $w + $DZ*$wz[$r] if numberp($wz[$r]);
+
+ $dc_u = $dc_u[$r] = $dc_u + $DZ*$dc_uz[$r] if numberp($dc_uz[$r]);
+ $dc_v = $dc_v[$r] = $dc_v + $DZ*$dc_vz[$r] if numberp($dc_vz[$r]);
+ $dc_w = $dc_w[$r] = $dc_w + $DZ*$dc_wz[$r] if numberp($dc_wz[$r]);
+
+ $uc_u = $uc_u[$r] = $uc_u + $DZ*$uc_uz[$r] if numberp($uc_uz[$r]);
+ $uc_v = $uc_v[$r] = $uc_v + $DZ*$uc_vz[$r] if numberp($uc_vz[$r]);
+ $uc_w = $uc_w[$r] = $uc_w + $DZ*$uc_wz[$r] if numberp($uc_wz[$r]);
+}
+
+#======================================================================
+# Step 4: Reference Velocities
+#======================================================================
+
+my($refU,$refV,$refW,$dc_refU,$dc_refV,$dc_refW,$uc_refU,$uc_refV,$uc_refW);
+
+if (defined($opt_b)) { # reference to bottom-track profile
+ print(STDERR "Loading BT data from $opt_b...\n")
+ if ($opt_d);
+ open(BTF,$opt_b) || croak("$opt_b: $!\n");
+
+ my(@BTL) = &antsFileLayout(BTF);
+ my($BTdF,$BTndF,$BTuF,$BTvF,$BTwF,$BTusF,$BTvsF,$BTwsF);
+ for (my($f)=0; $f<@BTL; $f++) {
+ $BTdF = $f if ($BTL[$f] eq 'depth');
+ $BTuF = $f if ($BTL[$f] eq 'u');
+ $BTvF = $f if ($BTL[$f] eq 'v');
+ $BTwF = $f if ($BTL[$f] eq 'w');
+ $BTusF = $f if ($BTL[$f] eq 'u.sig');
+ $BTvsF = $f if ($BTL[$f] eq 'v.sig');
+ $BTwsF = $f if ($BTL[$f] eq 'w.sig');
+ $BTndF = $f if ($BTL[$f] eq 'ndata');
+ }
+ croak("$opt_b: not a valid BT file\n")
+ unless defined($BTdF) && defined($BTuF) && defined($BTvF) && defined($BTwF) &&
+ defined($BTusF) && defined($BTvsF) && defined($BTwsF) && defined($BTndF);
+
+ while (my(@BTr) = &antsFileIn(BTF)) {
+ my($gi) = int($BTr[$BTdF] / $DZ);
+ next unless ($BTr[$BTndF] >= $opt_m);
+ $BT_nsamp[$gi] = $BTr[$BTndF];
+ $BT_u[$gi] = $BTr[$BTuF];
+ $BT_v[$gi] = $BTr[$BTvF];
+ $BT_w[$gi] = $BTr[$BTwF];
+ $BT_u_var[$gi] = $BTr[$BTusF]**2;
+ $BT_v_var[$gi] = $BTr[$BTvsF]**2;
+ $BT_w_var[$gi] = $BTr[$BTwsF]**2;
+ }
+
+ &fixLowSampStat(\@BT_u_var,@BT_nsamp); # remove spurious small variances
+ &fixLowSampStat(\@BT_v_var,@BT_nsamp);
+ &fixLowSampStat(\@BT_w_var,@BT_nsamp);
+
+ my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW, # average integrated-shear velocities
+ $uc_sumU,$uc_sumV,$uc_sumW);
+ my($nSumVel,$dc_nSumVel,$uc_nSumVel);
+ my($wSumBTu,$wSumBTv,$wSumBTw); # weighted sums of BT-ref'd velocities
+ my($dc_wSumBTu,$dc_wSumBTv,$dc_wSumBTw);
+ my($uc_wSumBTu,$uc_wSumBTv,$uc_wSumBTw);
+ my($sumVarBTu,$sumVarBTv,$sumVarBTw); # sum of variances of BT-ref'd vels
+ my($dc_sumVarBTu,$dc_sumVarBTv,$dc_sumVarBTw);
+ my($uc_sumVarBTu,$uc_sumVarBTv,$uc_sumVarBTw);
+
+ for (my($r)=0; $r<@depth; $r++) {
+ if (numberp($BT_u[$r]) && numberp($u[$r])) {
+ $nSumVel++;
+ $sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r];
+ $wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $sumVarBTu += 1/$BT_u_var[$r];
+# printf(STDERR "v = $BT_v[$r] (w = %.1f)\n",1/$BT_v_var[$r]);
+ $wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $sumVarBTv += 1/$BT_v_var[$r];
+ $wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $sumVarBTw += 1/$BT_w_var[$r];
+ }
+ if (numberp($BT_u[$r]) && numberp($dc_u[$r])) {
+ $dc_nSumVel++;
+ $dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r];
+ $dc_wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $dc_sumVarBTu += 1/$BT_u_var[$r];
+ $dc_wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $dc_sumVarBTv += 1/$BT_v_var[$r];
+ $dc_wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $dc_sumVarBTw += 1/$BT_w_var[$r];
+ }
+ if (numberp($BT_u[$r]) && numberp($uc_u[$r])) {
+ $uc_nSumVel++;
+ $uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r];
+ $uc_wSumBTu += $BT_u[$r] / $BT_u_var[$r]; $uc_sumVarBTu += 1/$BT_u_var[$r];
+ $uc_wSumBTv += $BT_v[$r] / $BT_v_var[$r]; $uc_sumVarBTv += 1/$BT_v_var[$r];
+ $uc_wSumBTw += $BT_w[$r] / $BT_w_var[$r]; $uc_sumVarBTw += 1/$BT_w_var[$r];
+ }
+ }
+
+ if ($nSumVel > 0) {
+ $refU = $sumU/$nSumVel - $wSumBTu/$sumVarBTu;
+ $refV = $sumV/$nSumVel - $wSumBTv/$sumVarBTv;
+ $refW = $sumW/$nSumVel - $wSumBTw/$sumVarBTw;
+
+ $dc_refU = $dc_sumU/$dc_nSumVel - $dc_wSumBTu/$dc_sumVarBTu;
+ $dc_refV = $dc_sUmV/$dc_nSumVel - $dc_wSumBTv/$dc_sumVarBTv;
+ $dc_refW = $dc_sumW/$dc_nSumVel - $dc_wSumBTw/$dc_sumVarBTw;
+
+ $uc_refU = $uc_sumU/$uc_nSumVel - $uc_wSumBTu/$uc_sumVarBTu;
+ $uc_refV = $uc_sUmV/$uc_nSumVel - $uc_wSumBTv/$uc_sumVarBTv;
+ $uc_refW = $uc_sumW/$uc_nSumVel - $uc_wSumBTw/$uc_sumVarBTw;
+ } else {
+ &antsInfo("$opt_b: no valid BT data --- calculating baroclinic solution only");
+ }
+}
+
+unless (defined($refU)) {
+ my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW,$uc_sumU,$uc_sumV,$uc_sumW);
+ my($nSumVel,$dc_nSumVel,$uc_nSumVel);
+
+ for (my($r)=0; $r<@depth; $r++) {
+ if (numberp($u[$r])) {
+ $nSumVel++;
+ $sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r];
+ }
+ if (numberp($dc_u[$r])) {
+ $dc_nSumVel++;
+ $dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r];
+ }
+ if (numberp($uc_u[$r])) {
+ $uc_nSumVel++;
+ $uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r];
+ }
+ }
+
+ $refU = $sumU / $nSumVel; $refV = $sumV / $nSumVel; $refW = $sumW / $nSumVel;
+ $dc_refU = $dc_sumU / $dc_nSumVel; $dc_refV = $dc_sumV / $dc_nSumVel; $dc_refW = $dc_sumW / $dc_nSumVel;
+ $uc_refU = $uc_sumU / $uc_nSumVel; $uc_refV = $uc_sumV / $uc_nSumVel; $uc_refW = $uc_sumW / $uc_nSumVel;
+}
+
+for (my($r)=0; $r<@depth; $r++) { # reference velocities
+ $u[$r] -= $refU if defined($u[$r]);
+ $v[$r] -= $refV if defined($v[$r]);
+ $w[$r] -= $refW if defined($w[$r]);
+ $dc_u[$r] -= $dc_refU if defined($dc_u[$r]);
+ $dc_v[$r] -= $dc_refV if defined($dc_v[$r]);
+ $dc_w[$r] -= $dc_refW if defined($dc_w[$r]);
+ $uc_u[$r] -= $uc_refU if defined($uc_u[$r]);
+ $uc_v[$r] -= $uc_refV if defined($uc_v[$r]);
+ $uc_w[$r] -= $uc_refW if defined($uc_w[$r]);
+}
+
+#======================================================================
+# Determine X Factor
+#======================================================================
+
+my($first_w,$last_w);
+for (my($r)=0; !defined($first_w) || !defined($last_w); $r++) {
+ $first_w = $dc_w[$r] unless defined($first_w);
+ $last_w = $uc_w[$r] unless defined($last_w);
+}
+
+my($X_Factor) = 100 * abs($last_w-$first_w) / sqrt(@depth / $DZ);
+&antsAddParams('X-Factor',$X_Factor);
+printf(STDERR "X-Factor = %.1f\n",$X_Factor);
+
+#======================================================================
+# Output Velocity Profile
+#======================================================================
+
+@antsNewLayout = ('depth','u','v','w','dc_u','dc_v','dc_w','uc_u','uc_v','uc_w');
+for (my($r)=0; $r<@depth; $r++) {
+ &antsOut($depth[$r]+$DZ/2,$u[$r],$v[$r],$w[$r],
+ $dc_u[$r],$dc_v[$r],$dc_w[$r],
+ $uc_u[$r],$uc_v[$r],$uc_w[$r]);
+}
+
+&antsExit();
new file mode 100755
--- /dev/null
+++ b/2011_07_10/LADCPproc
@@ -0,0 +1,625 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P P R O C
+# doc: Thu Sep 16 20:36:10 2010
+# dlm: Fri Jul 8 03:33:21 2011
+# (c) 2010 A.M. Thurnherr & E. Firing
+# uE-Info: 491 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'process LADCP data to get shear, time series';
+
+# NOTES:
+# - this code is based on merge.c written by Eric Firing
+# - comments starting with ## are taken from Eric's code
+# - for SeaBird files, CTD elapsed time is estimated from recno * CTD{sampint}
+# - CTD{elapsed} is undefined for records before instrument is in the water
+# - ITS-90 temp field in degC required
+# - salin field prequired
+# - pressure field in dbar required
+# - -i should be set to the number that's added to LADCP_elapsed to make the two
+# time series overplot nicely
+
+# HISTORY:
+# Sep 16, 2010: - incepted
+# Oct 13, 2010: - first working version
+# Oct 14, 2010: - renamed from LADCPshear
+# Oct 19, 2010: - added -a)coustic backscatter profiles
+# Oct 20, 2010: - added -2)dary CTD sensors
+# Oct 23, 2010: - added magnetic-declination correction
+# Oct 26, 2010: - added tilt calculation
+# Dec 9, 2010: - added support for ASCII CTD files
+# Dec 10, 2010: - change -w) default to 120s
+# - changed nshear output to 0 from nan when there are no samples
+# Dec 27, 2010: - changed sign of -l to accept lag output from [LADCP_w]
+# Jan 10, 2011: - -o => -k added new -o
+# - added code to fill CTD sound vel gaps
+# Jan 22, 2011: - added -g) lat,lon
+# - added -c)ompass corr
+# Jun 15, 2011: - added mean backscatter profile to default output
+# Jul 7, 2011: - added support for $BT_* processing parameters
+# - replaced old per-bin acoustic backscatter profile output by
+# acoustic backscatter depth-time-series
+# - disabled seabed and acoustic backscatter code when not required (e.g. UL)
+# - made non-diagnostic output terser
+
+($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
+($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
+($LADCPPROC) = ($0 =~ m{^(.*)/[^/]*$});
+
+require "$ANTS/ants.pl";
+require "$ANTS/libEOS83.pl";
+require "$ANTS/libstats.pl";
+require "$LADCPPROC/LADCPproc.loadCTD";
+require "$LADCPPROC/LADCPproc.bestLag";
+require "$LADCPPROC/LADCPproc.BT";
+require "$LADCPPROC/LADCPproc.backscatter";
+require "$LADCPPROC/LADCPproc.UHcode";
+require "$PERL_TOOLS/RDI_BB_Read.pl";
+require "$PERL_TOOLS/RDI_Coords.pl";
+require "$PERL_TOOLS/RDI_Utils.pl";
+
+$antsParseHeader = 0;
+&antsUsage('24a:b:c:df:g:i:kl:n:o:ps:t:w:',2,
+ '[use -2)dary CTD sensor pair]',
+ '[require -4)-beam LADCP solutions]',
+ '[-s)etup <file>] [-g)ps <lat,lon>]',
+ '[-c)ompass-corr <offset,cos-fac,sin-fac>]',
+ '[enable -p)PI editing]',
+ '[-o)utput grid <resolution[5m]>]',
+ '[-i)nitial LADCP time lag <guestimate>]',
+ '[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
+ '[-d)iagnostic output]',
+ 'output: [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
+ ' [-a)coustic backscatter <depth-time-series file] [bottom-trac-k) profs]',
+ '<RDI file> <SeaBird file>');
+
+$RDI_Coords::minValidVels = 4 if ($opt_4);
+
+&antsFloatOpt($opt_l);
+&antsCardOpt(\$opt_w,120);
+ # old default of -w 30 does not work if there are significant ambiguity-velocity
+ # problems, as is the case, e.g., with 2010_DIMES_US2 station 142
+ # old default of -w 60 did not work for DIMES_UK2 station 4 (DL), possibly again
+ # related to ambiguity velocity
+&antsCardOpt(\$opt_n,20);
+&antsFileOpt($opt_s);
+&antsFloatOpt($opt_i);
+&antsCardOpt($opt_o);
+
+if (defined($opt_g)) {
+ ($CTD{lat},$CTD{lon}) = split(',',$opt_g);
+ croak("$0: cannot decode -g $opt_g\n")
+ unless numberp($CTD{lat}) && numberp($CTD{lon});
+}
+
+if (defined($opt_c)) {
+ ($CC_offset,$CC_cos_fac,$CC_sin_fac) = split(',',$opt_c);
+ croak("$0: cannot decode -c $opt_c\n")
+ unless numberp($CC_offset) && numberp($CC_cos_fac) && numberp($CC_sin_fac);
+}
+
+
+$LADCP_file = &antsFileArg();
+$CTD_file = &antsFileArg();
+
+&antsAddParams('LADCP_file',$LADCP_file,'CTD_file',$CTD_file);
+&antsActivateOut();
+
+#----------------------------------------------------------------------
+# Step 1: Read LADCP Data
+#----------------------------------------------------------------------
+
+print(STDERR "Reading LADCP data ($LADCP_file)...");
+readData($LADCP_file,\%LADCP);
+printf(STDERR "\n\t%d ensembles",scalar(@{$LADCP{ENSEMBLE}})) if ($opt_d);
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 2: Set Processing Parameters
+#----------------------------------------------------------------------
+
+print(STDERR "Setting processing parameters...\n");
+
+printf(STDERR "\tloading $LADCPPROC/LADCPproc.defaults...\n");
+require "$LADCPPROC/LADCPproc.defaults";
+
+if (defined($opt_s)) {
+ print(STDERR "\tloading $opt_s...\n");
+ require $opt_s;
+}
+
+if ($LADCP{BLANKING_DISTANCE} == 0) {
+ print(STDERR "\t\tBLANKING_DISTANCE == 0 => excluding all data from bin 1\n")
+ if ($opt_d);
+ $wbin_start = 2 unless ($wbin_start > 2);
+ $ubin_start = 2 unless ($ubin_start > 2);
+ $shbin_start = 2 unless ($shbin_start > 2);
+ $BT_bin_start = 2 unless ($BT_bin_start > 2);
+}
+
+&antsAddParams('ADCP_orientation',
+ $LADCP{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'uplooker' : 'downlooker');
+
+$SHEAR_PREGRID_DZ = 20;
+$GRID_DZ = defined($opt_o) ? $opt_o : 5;
+
+my($year) = substr($LADCP{ENSEMBLE}[0]->{DATE},6,4);
+my($month) = substr($LADCP{ENSEMBLE}[0]->{DATE},0,2);
+my($dau ) = substr($LADCP{ENSEMBLE}[0]->{DATE},3,2);
+my($magdec,$maginc,$h_strength,$v_strength) = split('\s+',`magdec $CTD{lon} $CTD{lat} $year $month $day`);
+
+croak("$0: unknown magnetic declination\n")
+ unless defined($magdec);
+
+&antsAddParams('magnetic_declination',$magdec);
+
+#----------------------------------------------------------------------
+# Step 3: Read CTD data
+#----------------------------------------------------------------------
+
+print(STDERR "Reading CTD data ($CTD_file)...");
+readCTD($CTD_file,\%CTD);
+printf(STDERR "\n\t%d scans",scalar(@{$CTD{press}})) if ($opt_d);
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 4: Pre-Process CTD & LADCP Data
+#----------------------------------------------------------------------
+
+printf(STDERR "Pre-processing data...");
+printf(STDERR "\n\tCTD...") if ($opt_d);
+
+#------------------------
+# clean CTD pressure data
+#------------------------
+my($pSpikes) = 0;
+for (my($r)=1; $r<@{$CTD{press}}; $r++) {
+ $pSpikes++,$CTD{press}[$r]=nan
+ if (abs($CTD{press}[$r]-$CTD{press}[$r-1])/$CTD{sampint} > 2);
+}
+print(STDERR "\n\t\t$pSpikes pressure spikes removed")
+ if ($pSpikes>0 && $opt_d);
+
+#------------------------------------
+# calculate w and find deepest record
+#------------------------------------
+$CTD{maxpress} = -9e99;
+for (my($r)=1; $r<@{$CTD{press}}-1; $r++) {
+ $CTD{w}[$r] = 0.99*($CTD{press}[$r+1] - $CTD{press}[$r-1]) / (2*$CTD{sampint});
+ if ($CTD{press}[$r] > $CTD{maxpress}) {
+ $CTD{maxpress} = $CTD{press}[$r];
+ $CTD{atbottom} = $r;
+ }
+}
+printf(STDERR "\n\t\tmax pressure [%ddbar] at scan#%d",$CTD{maxpress},$CTD{atbottom})
+ if $opt_d;
+
+#----------------------------------------------------------------------
+# Step 4b: Pre-Process LADCP Data
+#----------------------------------------------------------------------
+
+print(STDERR "\n\tLADCP...") if ($opt_d);
+
+#-------------------------------------------
+# transform to earth coordinates if required
+#-------------------------------------------
+
+$U = 0; # velocity indices
+$V = 1;
+$W = 2;
+$E = 3;
+
+$LADCP{HEADING_BIAS} = -$magdec;
+
+if ($LADCP{BEAM_COORDINATES}) {
+ print(STDERR "\n\t\ttransforming beam to Earth coordinates...")
+ if ($opt_d);
+ for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{TILT} = &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
+ velInstrumentToEarth(\%LADCP,$ens,velBeamToInstrument(\%LADCP,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
+ @{$LADCP{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$bin]} = # fake it to fool ref_lr_w
+ (0,0,0,defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]) ? 100 : 0);
+ }
+ }
+ $LADCP{BEAM_COORDINATES} = 0;
+ $LADCP{EARTH_COORDINATES} = 1;
+ unless ($opt_4) {
+ print(STDERR "\n\t\t\t3-beam solutions: $RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4\n")
+ if ($opt_d);
+ &antsAddParams('3_beam_solutions',"$RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4");
+ }
+} elsif ($LADCP{EARTH_COORDINATES}) {
+ if ($opt_d) {
+ if ($opt_c) {
+ printf(STDERR "\n\t\tcalculating tilt and correcting for compass error and magnetic declination of %.1f deg...\n",$magdec);
+ } else {
+ printf(STDERR "\n\t\tcalculating tilt and correcting for magnetic declination of %.1f deg...\n",$magdec);
+ }
+ }
+ for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
+ $LADCP{ENSEMBLE}[$ens]->{TILT} = &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+ my($hdg) = rad($LADCP{ENSEMBLE}[$ens]->{HEADING});
+ $LADCP{HEADING_BIAS} = ($CC_offset + $CC_cos_fac*cos($hdg) + $CC_sin_fac*sin($hdg)) - $magdec
+ if ($opt_c);
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ @{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
+ velApplyHdgBias(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
+ }
+ }
+} else {
+ croak("$0: can only handle beam or earth coordinates\n")
+}
+
+#------------------------------------------------------
+# construct a depth-vs-time "profile" from the raw data
+#------------------------------------------------------
+
+print(STDERR "\t\tconstructing profile time series...")
+ if ($opt_d);
+
+($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
+ mk_prof(\%LADCP,0,undef,1,6,70,0.1,9999);
+croak("\n$LADCP_file: no good ensembles found\n")
+ unless defined($LADCP_start);
+
+if ($opt_d) {
+ printf(STDERR "\n\t\t\tStart of cast : %s (#%5d) at %6.1fm\n",
+ $LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_start]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH});
+ printf(STDERR "\t\t\tBottom of cast : %s (#%5d) at %6.1fm\n",
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+ printf(STDERR "\t\t\tEnd of cast : %s (#%5d) at %6.1fm",
+ $LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
+ $LADCP{ENSEMBLE}[$LADCP_end]->{NUMBER},
+ $LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
+}
+
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 5: Add CTD to LADCP Data & correct velocities for sound speed
+# - {DEPTH} field is overwritten!
+#----------------------------------------------------------------------
+
+print(STDERR "Matching CTD to LADCP time series...");
+
+$opt_l = defined($opt_l) ? -$opt_l : &lagLADCP2CTD();
+
+print(STDERR "Associating CTD data with LADCP ensembles...");
+
+for (my($min_depth)=9e99,my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ my($lastSvel);
+ my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - $opt_l) / $CTD{sampint});
+ if ($r < 0 && $ens == $LADCP_start) {
+ $r = int(($LADCP{ENSEMBLE}[++$ens]->{ELAPSED_TIME} - $opt_l) / $CTD{sampint})
+ while ($r < 0);
+ printf(STDERR "\n\tCTD data begin with instrument already in water => skipping %ds of LADCP data",
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME});
+ $LADCP_start = $ens;
+ }
+ if ($r > $#{$CTD{press}}) {
+ printf(STDERR "\n\tCTD data end while instrument is still in water => truncating %ds of LADCP data",
+ $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME})
+ if ($opt_d);
+ $LADCP_end = $ens - 1;
+ last;
+ }
+ my($dr);
+ for ($dr=0; !numberp($CTD{press}[$r+$dr]); $dr--) {}
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH} = depth($CTD{press}[$r+$dr],$CTD{lat});
+ if ($LADCP{ENSEMBLE}[$ens]->{DEPTH} < $min_depth) {
+ $min_depth = $LADCP{ENSEMBLE}[$ens]->{DEPTH};
+ $LADCP_top = $ens;
+ }
+ $LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD{w}[$r];
+ $LADCP{ENSEMBLE}[$ens]->{CTD_TEMP} = $CTD{temp}[$r];
+ $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} = sVel($CTD{salin}[$r],$CTD{temp}[$r],$CTD{press}[$r+$dr]);
+ if (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_SVEL})) {
+ $lastSvel = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL};
+ } else {
+ $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} = $lastSvel;
+ }
+ my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]);
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U] *= $sscorr;
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V] *= $sscorr;
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] *= $sscorr;
+ }
+}
+
+&antsAddParams('top_depth',,round($LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH}),
+ 'bottom_depth',round($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}),
+ 'start_date',$LADCP{ENSEMBLE}[$LADCP_start]->{DATE},
+ 'start_time',$LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
+ 'bottom_date',$LADCP{ENSEMBLE}[$LADCP_bottom]->{DATE},
+ 'bottom_time',$LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
+ 'end_date',$LADCP{ENSEMBLE}[$LADCP_end]->{DATE},
+ 'end_time',$LADCP{ENSEMBLE}[$LADCP_end]->{TIME});
+
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 6: Calculate Acoustic Backscatter Profile
+#----------------------------------------------------------------------
+
+print(STDERR "Calculating acoustic backscatter profiles...");
+mk_backscatter_profs($LADCP_start,$LADCP_end);
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 7: Find Seabed
+#----------------------------------------------------------------------
+
+if ($LADCP{ENSEMBLE}[$LADCP_start]->{XDUCER_FACING_DOWN}) {
+
+ print(STDERR "Finding seabed...");
+
+ print(STDERR "\n\tin acoustic backscatter profiles...") if ($opt_d);
+ ($water_depth,$sig_water_depth) =
+ find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+
+ $min_hab = $water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
+ printf(STDERR "\n\twater depth = %d(+-%.1f)m",$water_depth,$sig_water_depth);
+ printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
+
+ print(STDERR "\n\tin RDI BT data...") if ($opt_d);
+ ($BT_water_depth,$sig_BT_water_depth) =
+ find_seabed(\%LADCP,$LADCP_bottom,$LADCP{BEAM_COORDINATES});
+
+ if (defined($BT_water_depth)) {
+ $min_hab = $BT_water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
+ printf(STDERR "\n\twater depth = %d(+-%.1f)m",$BT_water_depth,$sig_BT_water_depth);
+ printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
+# $water_depth = $BT_water_depth; # assume BT data are better
+# $sig_water_depth = $sig_BT_water_depth; # (at least they are higher vertical resolution)
+ }
+
+ unless (defined($water_depth)) {
+ print(STDERR "\n\tno seabed found\n");
+ print(STDERR "\n\tunknown water depth => PPI editing disabled\n")
+ if ($opt_d);
+ $clip_margin = 0;
+ }
+
+ print(STDERR "\n");
+
+} else { # UPLOOKER
+ $water_depth = $sig_water_depth = $min_hab = nan;
+ print(STDERR "Uplooker data => PPI editing disabled\n")
+ if ($opt_d);
+}
+
+#----------------------------------------------------------------------
+# Step 8: Edit Data
+#----------------------------------------------------------------------
+
+print(STDERR "Calculating shear profiles...");
+
+$LADCP_start = 1 if ($LADCP_start == 0); # ensure that there is previous ensemble
+
+print(STDERR "\n\tdowncast...") if ($opt_d);
+edit_velocity($LADCP_start,$LADCP_bottom); # downcast
+calc_shear($LADCP_start,$LADCP_bottom,$SHEAR_PREGRID_DZ,0); # pre-grid shear @SHEAR_PREGRID_DZm resolution
+calc_shear($LADCP_start,$LADCP_bottom,$GRID_DZ,1); # calculate final gridded shear profile
+
+@dc_sh_n = @sh_n; # save downcast results
+@dc_ush_mu = @ush_mu; @dc_ush_sig = @ush_sig;
+@dc_vsh_mu = @vsh_mu; @dc_vsh_sig = @vsh_sig;
+@dc_wsh_mu = @wsh_mu; @dc_wsh_sig = @wsh_sig;
+
+print(STDERR "\n\tupcast...") if ($opt_d);
+edit_velocity($LADCP_end,$LADCP_bottom); # upcast
+calc_shear($LADCP_end,$LADCP_bottom,$SHEAR_PREGRID_DZ,0);
+calc_shear($LADCP_end,$LADCP_bottom,$GRID_DZ,1);
+
+@uc_sh_n = @sh_n; # save upcast results
+@uc_ush_mu = @ush_mu; @uc_ush_sig = @ush_sig;
+@uc_vsh_mu = @vsh_mu; @uc_vsh_sig = @vsh_sig;
+@uc_wsh_mu = @wsh_mu; @uc_wsh_sig = @wsh_sig;
+
+print(STDERR "\n\tcombined...") if ($opt_d);
+for (my($gi)=0; $gi<@dc_ush_mu; $gi++) {
+ if ($dc_sh_n[$gi]>0 && $uc_sh_n[$gi]>0) {
+ $sh_n[$gi] = $dc_sh_n[$gi] + $uc_sh_n[$gi];
+ $ush_mu[$gi] = ($dc_sh_n[$gi]*$dc_ush_mu[$gi] + $uc_sh_n[$gi]*$uc_ush_mu[$gi]) / $sh_n[$gi];
+ $vsh_mu[$gi] = ($dc_sh_n[$gi]*$dc_vsh_mu[$gi] + $uc_sh_n[$gi]*$uc_vsh_mu[$gi]) / $sh_n[$gi];
+ $wsh_mu[$gi] = ($dc_sh_n[$gi]*$dc_wsh_mu[$gi] + $uc_sh_n[$gi]*$uc_wsh_mu[$gi]) / $sh_n[$gi];
+ $ush_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_ush_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_ush_sig[$gi]**2) / $sh_n[$gi]);
+ $vsh_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_vsh_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_vsh_sig[$gi]**2) / $sh_n[$gi]);
+ $wsh_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_wsh_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_wsh_sig[$gi]**2) / $sh_n[$gi]);
+ } elsif ($dc_sh_n[$gi] > 0) {
+ $sh_n[$gi] = $dc_sh_n[$gi];
+ $ush_mu[$gi] = $dc_ush_mu[$gi]; $vsh_mu[$gi] = $dc_vsh_mu[$gi]; $wsh_mu[$gi] = $dc_wsh_mu[$gi];
+ $ush_sig[$gi] = $dc_ush_sig[$gi]; $vsh_sig[$gi] = $dc_vsh_sig[$gi]; $wsh_sig[$gi] = $dc_wsh_sig[$gi];
+ } elsif ($uc_sh_n[$gi] > 0) {
+ $sh_n[$gi] = $uc_sh_n[$gi];
+ $ush_mu[$gi] = $uc_ush_mu[$gi]; $vsh_mu[$gi] = $uc_vsh_mu[$gi]; $wsh_mu[$gi] = $uc_wsh_mu[$gi];
+ $ush_sig[$gi] = $uc_ush_sig[$gi]; $vsh_sig[$gi] = $uc_vsh_sig[$gi]; $wsh_sig[$gi] = $uc_wsh_sig[$gi];
+ } else {
+ $sh_n[$gi] = 0;
+ $ush_mu[$gi] = $vsh_mu[$gi] = $wsh_mu[$gi] = nan;
+ $ush_sig[$gi] = $vsh_sig[$gi] = $wsh_sig[$gi] = nan;
+ }
+}
+
+print(STDERR "\n");
+
+#----------------------------------------------------------------------
+# Step 9: Get bottom track profile
+#----------------------------------------------------------------------
+
+if ($LADCP{ENSEMBLE}[$LADCP_start]->{XDUCER_FACING_DOWN}) {
+ print(STDERR "Getting BT profile...");
+ getBTprof($LADCP_start,$LADCP_end);
+ print(STDERR "\n");
+}
+
+#----------------------------------------------------------------------
+# Step 10: Write Output Files
+#----------------------------------------------------------------------
+
+print(STDERR "Writing shear profiles...");
+
+@antsNewLayout = ('depth','dc_nshear','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
+ 'uc_nshear','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
+ 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n');
+
+$commonParams = $antsCurParams;
+&antsAddParams('ubin_start',$ubin_start,'ubin_end',$ubin_end, # record processing params
+ 'wbin_start',$wbin_start,'wbin_end',$wbin_end,
+ 'shbin_start',$shbin_start,'shbin_end',$shbin_end,
+ 'w_ref_bin',$w_ref_bin,'w_dif',$w_dif,
+ 'wake_hd_dif',$wake_hd_dif,'wake_ang_min',$wake_ang_min,
+ 'min_wake_w',$min_wake_w,'n_wake_bins',$n_wake_bins,
+ 'e_max',$e_max,'min_cor',$min_cor,
+ 'max_shdev',$max_shdev,'max_shdev_sum',$max_shdev_sum,
+ 'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth),
+ 'min_hab',round($min_hab),
+ 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin,
+ 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end,
+ 'BT_bin_start',$BT_bin_start,'BT_bin_search_above',$BT_bin_search_above,
+ 'BT_max_bin_spread',$BT_max_bin_spread,'BT_max_w_difference',$BT_max_w_difference,
+);
+if (defined($BT_min_depth)) {
+ &antsAddParams('BT_min_depth',$BT_min_depth,'BT_max_depth',$BT_max_depth);
+} else {
+ &antsAddParams('BT_max_depth_error',$BT_max_depth_error);
+}
+
+for (my($gi)=0; $gi<@ush_mu; $gi++) {
+ &antsOut(depthOfGI($gi), # depth in center of bin
+ numberp($dc_sh_n[$gi])?$dc_sh_n[$gi]:0, # downcast
+ $dc_ush_mu[$gi],$dc_ush_sig[$gi],
+ $dc_vsh_mu[$gi],$dc_vsh_sig[$gi],
+ $dc_wsh_mu[$gi],$dc_wsh_sig[$gi],
+ numberp($uc_sh_n[$gi])?$uc_sh_n[$gi]:0, # upcast
+ $uc_ush_mu[$gi],$uc_ush_sig[$gi],
+ $uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
+ $uc_wsh_mu[$gi],$uc_wsh_sig[$gi],
+ $sh_n[$gi], # combined
+ $ush_mu[$gi],$ush_sig[$gi],
+ $vsh_mu[$gi],$vsh_sig[$gi],
+ $wsh_mu[$gi],$wsh_sig[$gi],
+ $nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan,
+ $nSv_prof[$gi],
+ );
+}
+
+print(STDERR "\n");
+
+#---------------------------------------
+# Acoustic backscatter depth-time-series
+#---------------------------------------
+
+if (defined($opt_a)) {
+ print(STDERR "Writing acoustic backscatter depth-time-series to <$opt_a>...");
+
+
+ @antsNewLayout = ('ens','elapsed','CTD_depth','depth','bin','downcast',
+ 'Sv_beam1','Sv_beam2','Sv_beam3','Sv_beam4','Sv.median');
+ &antsOut('EOF');
+ close(STDOUT);
+ open(STDOUT,">$opt_a") || croak("$opt_a: $!\n");
+
+ $antsCurParams = $commonParams;
+ &antsAddParams('min_elapsed',$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME},
+ 'max_elapsed',$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+ 'min_depth',$LADCP{ENSEMBLE}[$LADCP_top]->{XDUCER_FACING_UP} ?
+ &depthOfBin($LADCP_top,$LADCP{N_BINS}-1) : $LADCP{ENSEMBLE}[$LADCP_top]->{DEPTH},
+ 'max_depth',$LADCP{ENSEMBLE}[$LADCP_bottom]->{XDUCER_FACING_UP} ?
+ $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH} : &depthOfBin($LADCP_bottom,$LADCP{N_BINS}-1)
+ );
+
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ &antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH},
+ &depthOfBin($ens,$bin),$bin,
+ ($ens <= $LADCP_bottom) ? 1 : 0,
+ @{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]},
+ median(@{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]})
+ );
+ }
+ }
+
+ print(STDERR "\n");
+}
+
+#----------------------------------------------------------------------
+
+if (defined($opt_t)) {
+ print(STDERR "Writing time series to <$opt_t>...");
+
+ @antsNewLayout = ('ens','elapsed','depth','CTD_w','LADCP_w');
+ &antsOut('EOF');
+ $antsCurParams = $commonParams;
+ close(STDOUT);
+ open(STDOUT,">$opt_t") || croak("$opt_t: $!\n");
+
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ &antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH},
+ $LADCP{ENSEMBLE}[$ens]->{CTD_W},
+ $LADCP{ENSEMBLE}[$ens]->{W});
+ }
+ print(STDERR "\n");
+}
+
+#----------------------------------------------------------------------
+
+if (defined($opt_b)) {
+ print(STDERR "Writing bottom-track data to <$opt_b>...");
+
+ @antsNewLayout = ('depth','u','v','w','u.sig','v.sig','w.sig','ndata');
+ &antsOut('EOF');
+ $antsCurParams = $commonParams;
+ close(STDOUT);
+ open(STDOUT,">$opt_b") || croak("$opt_b: $!\n");
+
+ my($skipped);
+ for (my($gi)=0; $gi<@BT_nsamp; $gi++) {
+ $skipped = 1 if ($BT_nsamp[$gi] > 0);
+ next unless ($skipped);
+ &antsOut(depthOfGI($gi),$BTu[$gi],$BTv[$gi],$BTw[$gi],$BTu_sig[$gi],$BTv_sig[$gi],$BTw_sig[$gi],$BT_nsamp[$gi]);
+ }
+ print(STDERR "\n");
+}
+
+#----------------------------------------------------------------------
+
+if (defined($opt_f)) {
+ print(STDERR "Writing data flags to <$opt_f>...");
+
+ @antsNewLayout = ('ens');
+ for (my($i)=1; $i<=$LADCP{N_BINS}; $i++) {
+ $antsNewLayout[$i] = "bin$i";
+ }
+ &antsOut('EOF');
+ $antsCurParams = $commonParams;
+
+ close(STDOUT);
+ open(STDOUT,">$opt_f") || croak("$opt_f: $!\n");
+
+ &antsPrintHeaders(STDOUT,@antsNewLayout);
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ printf('%4d ',$LADCP{ENSEMBLE}[$ens]->{NUMBER});
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ printf("%02x ",$edit_flags[$ens][$bin]);
+ }
+ print($opt_R);
+ }
+
+ print(STDERR "\n");
+}
+
+&antsExit();
+
new file mode 100644
--- /dev/null
+++ b/2011_07_10/LADCPproc.BT
@@ -0,0 +1,187 @@
+#======================================================================
+# L A D C P P R O C . B T
+# doc: Wed Oct 20 21:05:37 2010
+# dlm: Fri Jul 8 03:24:56 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 148 69 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Oct 20, 2010: - created
+# Jan 10, 2010: - -o => -k
+# Jul 7, 2010: - added $DEBUG
+# - added BTrangeFlag
+# - added $BT processing parameters
+# - changed from echo amplitude to Sv
+
+my($BEAM1) = 0;
+my($BEAM2) = 1;
+my($BEAM3) = 2;
+my($BEAM4) = 3;
+
+my($nBTfound,$nBTrangeFlag,$nBTdepthFlag,$nBTvalidVelFlag,$nBTwFlag) = (0,0,0,0,0);
+
+my($DEBUG) = 0;
+
+sub binBTprof($)
+{
+ my($ens) = @_;
+
+ my(@Sv_max) = (-9e99,-9e99,-9e99,-9e99); my(@Sv_max_bin) = (nan,nan,nan,nan);
+ for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
+ if (defined($BT_min_depth)) { # manually supplied bottom depth range
+ my($dob) = &depthOfBin($ens,$bin);
+ next unless ($dob >= $BT_min_depth && $dob <= $BT_max_depth);
+ }
+ $Sv_max[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1],
+ $Sv_max_bin[$BEAM1] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1] > $Sv_max[$BEAM1]);
+ $Sv_max[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2],
+ $Sv_max_bin[$BEAM2] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2] > $Sv_max[$BEAM2]);
+ $Sv_max[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3],
+ $Sv_max_bin[$BEAM3] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3] > $Sv_max[$BEAM3]);
+ $Sv_max[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4],
+ $Sv_max_bin[$BEAM4] = $bin
+ if ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4] > $Sv_max[$BEAM4]);
+ }
+
+ print(STDERR "@Sv_max | @Sv_max_bin\n") if ($DEBUG);
+ $nBTfound++;
+
+ $nBTrangeFlag++,return # inconsistent range (&, impliclity, large tilt)
+ unless (max(@Sv_max_bin)-min(@Sv_max_bin) <= $BT_max_bin_spread);
+
+ my($range_bin) = round(avg(@Sv_max_bin));
+ printf(STDERR "water_depth = $water_depth; BT peak depth = %d in bin $range_bin\n",depthOfBin($ens,$range_bin))
+ if ($DEBUG);
+
+ $nBTdepthFlag++,return # BT range inconsistent with water depth
+ unless defined($BT_min_depth) ||
+ (abs($water_depth-depthOfBin($ens,$range_bin)) < $sig_water_depth + $BT_max_depth_error);
+
+ # try bin of max plus one above and below
+ # this does not really work because, often, only one of the bins has valid velocities
+ my($w1) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];
+ my($w2) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
+ my($w3) = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
+
+ printf(STDERR "w123 = %.1f,%.1f,%.1f\n",$w1,$w2,$w3)
+ if ($DEBUG);
+
+ $w1 = 9e99 unless numberp($w1); # no valid velocities
+ $w2 = 9e99 unless numberp($w1);
+ $w3 = 9e99 unless numberp($w1);
+
+ my($CTD_u,$CTD_v,$CTD_w);
+
+ if (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2) &&
+ abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w3)) {
+ $CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$U];
+ $CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$V];
+ $CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin-1][$W];
+ } elsif (abs($LADCP{ENSEMBLE}[$ens]->{W}-$w1) < abs($LADCP{ENSEMBLE}[$ens]->{W}-$w2)) {
+ $CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$U];
+ $CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$V];
+ $CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin+1][$W];
+ } else {
+ $nBTvalidVelFlag++,return if ($w2 == 9e99); # none of 3 bins has valid velocity
+ $CTD_u = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$U];
+ $CTD_v = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$V];
+ $CTD_w = $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$range_bin][$W];
+ }
+
+ $nBTwFlag++,return if (abs($CTD_w-$LADCP{ENSEMBLE}[$ens]->{W}) > $BT_max_w_difference);
+
+ printf(STDERR "good BT [%5.2f %5.2f %5.2f] found at ens $ens\n",$CTD_u,$CTD_v,$CTD_w)
+ if ($DEBUG);
+
+ if ($opt_k) {
+ for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
+ next if ($edit_flags[$ens][$bin]);
+ printf(BTF "%d %d %d %f %f %f %f %f %f %f %f %f %f %f\n",
+ $LADCP{ENSEMBLE}[$ens]->{NUMBER},
+ depthOfBin($ens,$bin),$LADCP{ENSEMBLE}[$ens]->{DEPTH},
+ $LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL},
+ $CTD_u,$CTD_v,$CTD_w,
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U],
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V],
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W],
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u,
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$CTD_v,
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]-$CTD_w);
+ }
+ print(BTF "nan nan nan nan\n");
+ }
+
+ for (my($bin)=$BT_bin_start-1; $bin<$LADCP{N_BINS}; $bin++) {
+ next if ($edit_flags[$ens][$bin]);
+ my($gi) = depthOfBin($ens,$bin) / $GRID_DZ;
+ push(@{$BTu_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U]-$CTD_u);
+ push(@{$BTv_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V]-$CTD_v);
+ push(@{$BTw_vals[$gi]},$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]-$CTD_w);
+ }
+
+}
+
+sub getBTprof($$)
+{
+ my($LADCP_start,$LADCP_end) = @_;
+
+ if ($opt_k) {
+ open(BTF,">BT.profs");
+ print(BTF "#ANTS#FIELDS# {ens} {depth} {CTD_depth} {pitch} {roll} {CTD_u} {CTD_v} {CTD_w} {u} {v} {w} {BT_u} {BT_v} {BT_w}\n");
+ }
+
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ next unless ($water_depth-$LADCP{ENSEMBLE}[$ens]->{DEPTH} < $BT_begin_search_above);
+ binBTprof($ens);
+ }
+
+ if ($opt_d) {
+ print(STDERR "\n\t$nBTfound BT ensembles found\n");
+ print(STDERR "\t\t$nBTrangeFlag flagged bad because of inconsistent range to seabed\n");
+ print(STDERR "\t\t$nBTdepthFlag flagged bad because of wrong bottom depth\n");
+ print(STDERR "\t\t$nBTvalidVelFlag flagged bad because of lack of valid velocities\n");
+ print(STDERR "\t\t$nBTwFlag flagged bad because of incorrect vertical velocities");
+ printf(STDERR "\n\t=> %d velocities from %d BT ensembles used",
+ scalar(@BTu_vals),
+ $nBTfound-$nBTrangeFlag-$nBTdepthFlag-$nBTvalidVelFlag-$nBTwFlag);
+ }
+
+ @BTu = @BTv = @BTw = ();
+ @BTu_sig = @BTv_sig = @BTw_sig = ();
+ @BT_nsamp = ();
+
+ for (my($gi)=0; $gi<@BTu_vals; $gi++) { # calc grid means & stddev
+ my($sum_u,$sum_v,$sum_w);
+
+ $BT_nsamp[$gi] = @{$BTu_vals[$gi]};
+
+ for (my($vi)=0; $vi<$BT_nsamp[$gi]; $vi++) {
+ $sum_u += $BTu_vals[$gi][$vi];
+ $sum_v += $BTv_vals[$gi][$vi];
+ $sum_w += $BTw_vals[$gi][$vi];
+ }
+ $BTu[$gi] = $BT_nsamp[$gi] ? $sum_u/$BT_nsamp[$gi] : nan;
+ $BTv[$gi] = $BT_nsamp[$gi] ? $sum_v/$BT_nsamp[$gi] : nan;
+ $BTw[$gi] = $BT_nsamp[$gi] ? $sum_w/$BT_nsamp[$gi] : nan;
+ }
+
+ for (my($gi)=0; $gi<@BTu_vals; $gi++) { # calc & grid stddevs
+ my($sumsq_u,$sumsq_v,$sumsq_w);
+ for (my($vi)=0; $vi<$BT_nsamp[$gi]; $vi++) {
+ $sumsq_u += ($BTu_vals[$gi][$vi] - $BTu[$gi])**2;
+ $sumsq_v += ($BTv_vals[$gi][$vi] - $BTv[$gi])**2;
+ $sumsq_w += ($BTw_vals[$gi][$vi] - $BTw[$gi])**2;
+ }
+ $BTu_sig[$gi] = $BT_nsamp[$gi]>1 ? sqrt($sumsq_u/($BT_nsamp[$gi]-1)) : nan;
+ $BTv_sig[$gi] = $BT_nsamp[$gi]>1 ? sqrt($sumsq_v/($BT_nsamp[$gi]-1)) : nan;
+ $BTw_sig[$gi] = $BT_nsamp[$gi]>1 ? sqrt($sumsq_w/($BT_nsamp[$gi]-1)) : nan;
+ }
+
+ close(BTF) if ($opt_k);
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/2011_07_10/LADCPproc.UHcode
@@ -0,0 +1,439 @@
+#======================================================================
+# L A D C P P R O C . U H C O D E
+# doc: Fri Sep 17 20:27:53 2010
+# dlm: Mon Jan 10 13:23:13 2011
+# (c) 2010 A.M. Thurnherr & E. Firing
+# uE-Info: 295 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# PERLified functions from Eric's [merge.c]; with mods
+
+# NOTES:
+# - velocity integration removed
+# - percent-good flag removed (no point for single-ping ensembles)
+# - need the following per-ensemble fields from previous steps:
+# CTD_DEPTH
+# CTD_SVEL
+# - negative depths are not allowed (should not happen given CTD_DEPTH)
+
+# WEIRDNESSES:
+# - w reference layer in set_misc_flags is calculated without taking
+# the ref layer bins into account
+# - u,v ref lr vels in set_wake_flag use w ref layer bins
+# - distance to bin 1 center is not sound-speed corrected [CORRECTED]
+# - $tilt calculation is wrong
+
+# HISTORY:
+# Sep 17, 2010: - created
+# Oct 13, 2010: - first working version
+# Oct 14, 2010: - renamed from LADCPshear.UHcode
+# Oct 15, 2010: - added support for -pPI edit suppresion
+# Oct 19, 2010: - reversed semantics of -p
+# Oct 25, 2010: - added W_CTD_BIT, renamed W_BIT to W_OUTLIER_BIT
+# Oct 26, 2010: - added TILT_BIT
+# Dec 10, 2010: - modified assertion to allow processing of UL data
+
+#======================================================================
+# VELOCITY EDITING
+#======================================================================
+
+my($BADVEL_BIT) = 0x01;
+my($ERRVEL_BIT) = 0x02;
+my($CORREL_BIT) = 0x04;
+my($W_OUTLIER_BIT) = 0x08;
+my($SHEAR_BIT) = 0x10;
+my($SIDELOBE_BIT) = 0x20;
+my($WAKE_BIT) = 0x40;
+my($PPI_BIT) = 0x80;
+my($W_CTD_BIT) = 0x100;
+my($TILT_BIT) = 0x200;
+my($DELTA_TILT_BIT) = 0x400;
+
+my(%flag_count);
+
+# apparently, in Eric's code, DISTANCE_TO_BIN1_CENTER is not sound-speed corrected
+sub dzToBin($$)
+{
+ my($ens,$bin) = @_;
+ my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
+ return $sscorr * ($LADCP{DISTANCE_TO_BIN1_CENTER} + $bin*$LADCP{BIN_LENGTH});
+}
+
+sub depthOfBin($$)
+{
+ my($ens,$bin) = @_;
+ return $LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} ?
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH} - &dzToBin($ens,$bin) :
+ $LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBin($ens,$bin);
+}
+
+sub set_wake_flag($$)
+{
+ my($ens,$De) = @_;
+
+ my($n) = 0;
+ my($uref) = my($vref) = my($wref) = 0;
+
+ for (my($bin)=$wbin_start-1; $bin<$wbin_end; $bin++) { # calc crude ref lr vel from all data
+ next if ($edit_flags[$ens][$bin]);
+ $uref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U];
+ $vref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V];
+ $wref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W];
+ $n++;
+ }
+ return if ($n==0);
+ $uref /= $n;
+ $vref /= $n;
+ $wref /= $n;
+
+ ## if upward (=negative) velocity greater than minimum, calculate wake
+ ## heading and inclination
+ if ($wref < -$min_wake_w) {
+ my($wake_hd) = 180 / 3.14159265358979 * atan2($uref,$vref);
+ my($speed) = sqrt($uref*$uref + $vref*$vref);
+ my($wake_ang)= abs(180 / 3.14159265358979 * atan($speed/$wref));
+
+ $wake_hd += 360
+ if ($wake_hd < 0);
+ $LADCP{ENSEMBLE}[$ens]->{HEADING} += 360
+ if ($LADCP{ENSEMBLE}[$ens]->{HEADING} < 0);
+
+ my($wake_mod) = $wake_hd % 90; # % returns integer part of remainder, but that's sufficient
+ my($adcp_mod) = $LADCP{ENSEMBLE}[$ens]->{HEADING} % 90;
+
+ if (((abs($wake_mod - $adcp_mod) < $wake_hd_dif) ||
+ (abs($wake_mod - $adcp_mod) > (90 - $wake_hd_dif))) &&
+ ($wake_ang > $wake_ang_min)) {
+ for (my($bin)=0; $bin<$n_wake_bins; $bin++) {
+ $flag_count{$WAKE_BIT}[$bin]++;
+ $edit_flags[$ens][$bin] |= $WAKE_BIT;
+ }
+ }
+ } ## if ($wref < -min_wake_w)
+
+ ## This does not make a lot of sense, because it trims points
+ ## on only one side of the wake error, and that side depends
+ ## on whether the integration is forward or backward.
+
+ if (($edit_flags[$ens+$De][0]&$WAKE_BIT) &&
+ ($edit_flags[$ens][0]&$WAKE_BIT == 0)) {
+ for (my($bin)=0; $bin<$n_wake_bins; $bin++) {
+ $flag_count{$WAKE_BIT}[$bin]++;
+ $edit_flags[$ens][$bin] |= $WAKE_BIT;
+ }
+ }
+}
+
+sub set_misc_flags($$)
+{
+ my($ens,$De) = @_;
+ my($ww) = my($n) = 0;
+ my($SLIfac) = 1 - cos(rad($LADCP{BEAM_ANGLE}));
+
+ for (my($bin)=0; $bin<$w_ref_bin; $bin++) { # ref-lr w
+ if (!$edit_flags[$ens][$bin]) {
+ $ww += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W];
+ $n++;
+ }
+ }
+ $ww /= $n if ($n > 0);
+
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ next if ($edit_flags[$ens][$bin]);
+
+ ## We use the standard criterion for bottom interference; e.g. for
+ ## 30-degree beams, the last 15% of the profile is contaminated
+ ## by the sidelobe bounce. 1.5 bin length is added to allow for
+ ## the length of the bin and pulse, that is, contamination of part of a
+ ## bin. Profiler tilt does not require a more stringent criterion.
+ if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {
+ if (numberp($water_depth) &&
+ $water_depth - &depthOfBin($ens,$bin) <=
+ $SLIfac * ($water_depth - $LADCP{ENSEMBLE}[$ens]->{DEPTH})
+ + 1.5 * $LADCP{BIN_LENGTH}) {
+ $edit_flags[$ens][$bin] |= $SIDELOBE_BIT;
+ $flag_count{$SIDELOBE_BIT}++;
+ }
+ } else { ## upward-looking
+ if (&depthOfBin($ens,$bin) <=
+ $SLIfac * $LADCP{ENSEMBLE}[$ens]->{DEPTH}
+ + 1.5 * $LADCP{BIN_LENGTH}) {
+ $edit_flags[$ens][$bin] |= $SIDELOBE_BIT;
+ $flag_count{$SIDELOBE_BIT}++;
+ }
+ }
+
+ if ($ww != 0 &&
+ abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] - $ww) > $w_dif) {
+ $flag_count{$W_OUTLIER_BIT}++;
+ $edit_flags[$ens][$bin] |= $W_OUTLIER_BIT;
+ }
+
+ if (abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$E]) > $e_max) {
+ $flag_count{$ERRVEL_BIT}++;
+ $edit_flags[$ens][$bin] |= $ERRVEL_BIT;
+ }
+
+ for (my($beam)=0; $beam<=3; $beam++) {
+ if ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] < $min_cor) {
+ $flag_count{$CORREL_BIT}++;
+ $edit_flags[$ens][$bin] |= $CORREL_BIT;
+ }
+ }
+
+ if ($bin < $shbin_start-1 || $bin >= $shbin_end) {
+ $edit_flags[$ens][$bin] |= $BADVEL_BIT;
+ $flag_count{$BADVEL_BIT}++;
+ }
+ } # for ($bin=0...
+}
+
+## The following is for editing out the second bottom bounce.
+sub set_PPI_flags($$)
+{
+ my($ens,$De) = @_;
+
+ my($dt_ping) = $LADCP{ENSEMBLE}[$ens]->{UNIX_TIME} - $LADCP{ENSEMBLE}[$ens-1]->{UNIX_TIME};
+
+ if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {
+ if (numberp($water_depth)) {
+ $clip_z1 = $water_depth
+ - $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
+ * cos(rad($LADCP{BEAM_ANGLE} + $LADCP{ENSEMBLE}[$ens]->{TILT}))
+ + $clip_margin;
+ $clip_z0 = $water_depth
+ - $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
+ * cos(rad($LADCP{BEAM_ANGLE} - $LADCP{ENSEMBLE}[$ens]->{TILT}))
+ - $clip_margin;
+ }
+ } else { # upward-looking
+ $clip_z1 = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
+ * cos(rad($LADCP{BEAM_ANGLE} - $LADCP{ENSEMBLE}[$ens]->{TILT}))
+ + $clip_margin;
+ $clip_z0 = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping
+ * cos(rad($LADCP{BEAM_ANGLE} + $LADCP{ENSEMBLE}[$ens]->{TILT}))
+ - $clip_margin;
+ }
+
+ for (my($bin)=$first_clip_bin-1; $bin<$LADCP{N_BINS}; $bin++) {
+ next if ($edit_flags[$ens][$bin]);
+
+ my($dob) = depthOfBin($ens,$bin);
+ if ($dob >= $clip_z0 && $dob <= $clip_z1) {
+ $edit_flags[$ens][$bin] |= $PPI_BIT;
+ $flag_count{$PPI_BIT}++;
+ }
+ }
+}
+
+sub edit_velocity($$)
+{
+ my($start,$end) = @_; # ensemble indices
+ my($De) = $start<$end ? 1 : -1; # downcast: De = 1; upcast: De = -1
+
+ $flag_count{$WAKE_BIT} = $flag_count{$W_OUTLIER_BIT} = $flag_count{$ERRVEL_BIT} =
+ $flag_count{$CORREL_BIT} = $flag_count{$SHEAR_BIT} = $flag_count{$BADVEL_BIT} =
+ $flag_count{$SIDELOBE_BIT} = $flag_count{$PPI_BIT} = $flag_count{$W_CTD_BIT} =
+ $flag_count{$TILT_BIT} = $flag_count{$DELTA_TILT_BIT} = 0;
+
+ for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) { # loop over all ens from start to end
+ if (abs($LADCP{ENSEMBLE}[$ens]->{CTD_W}-$LADCP{ENSEMBLE}[$ens]->{W}) > $w_max_err) { # get rid of aliased vels (ambiguity)
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ $edit_flags[$ens][$bin] |= $W_CTD_BIT;
+ $flag_count{$W_CTD_BIT}++;
+ }
+ next;
+ }
+ if ($LADCP{ENSEMBLE}[$ens]->{TILT} > $max_tilt) { # get rid ensembles with large tilt
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ $edit_flags[$ens][$bin] |= $TILT_BIT;
+ $flag_count{$TILT_BIT}++;
+ }
+ next;
+ } # get rid ensembles after large rotation
+ if (abs($LADCP{ENSEMBLE}[$ens]->{TILT}-$LADCP{ENSEMBLE}[$ens-$De]->{TILT}) > $max_delta_tilt) {
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ $edit_flags[$ens][$bin] |= $DELTA_TILT_BIT;
+ $flag_count{$DELTA_TILT_BIT}++;
+ }
+ next;
+ }
+ for (my($bin)=$shbin_start-1; $bin<$shbin_end; $bin++) { # flag bad velocities
+ $edit_flags[$ens][$bin] |= $BADVEL_BIT,$flag_count{$BADVEL_BIT}++
+ unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]);
+ }
+ set_wake_flag($ens,$De);
+ set_misc_flags($ens,$De);
+ set_PPI_flags($ens,$De)
+ if $opt_p && ($clip_margin > 0); # PPI editing is off by default
+ }
+}
+
+#======================================================================
+# CALCULATE VELOCITY SHEAR
+# - output in @ush_mu,@vsh_mu,@wsh_mu,@ush_sig,@vsh_sig,@wsh_sig
+# - @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh are defined "local" in calc_shear
+#======================================================================
+
+sub uv_to_shear($)
+{
+ my($ens) = @_;
+ my($nvel) = 0;
+
+ @sh_i0 = @sh_i1 = ();
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ next if ($edit_flags[$ens][$bin]);
+ $nvel++;
+ push(@sh_i1,$bin) if (@sh_i0);
+ push(@sh_i0,$bin) if ($bin < $LADCP{N_BINS}-1);
+ }
+
+ @dsh = ();
+ for (my($i)=0; $i<@sh_i1; $i++) { # calc and bin shears
+ my($d0) = &depthOfBin($ens,$sh_i0[$i]);
+ my($d1) = &depthOfBin($ens,$sh_i1[$i]);
+ die("$0: assertion failed (ens=$ens i=$i depth=$LADCP{ENSEMBLE}[$ens]->{DEPTH} sh_i0[$i]=$sh_i0[$i] sh_i1[$i]=$sh_i1[$i] d0=$d0 d1=$d1)")
+ unless ($d0>=0 && $d1>=0);
+ die("$0: assertion failed (ens=$ens i=$i sh_i0[$i]=$sh_i0[$i] sh_i1[$i]=$sh_i1[$i] d0=$d0 d1=$d1)")
+ unless (($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN} && $d1-$d0>0) ||
+ ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} && $d1-$d0<0));
+ $dsh[$i] = ($d1 + $d0) / 2;
+ $ush[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$U] -
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$U]) / ($d1-$d0);
+ $vsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$V] -
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$V]) / ($d1-$d0);
+ $wsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$W] -
+ $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$W]) / ($d1-$d0);
+ }
+
+ return $nvel;
+}
+
+# here ush_mu, sh_n, etc are still set from the pre-gridding pass
+sub set_shear_flag($)
+{
+ my($ens) = @_;
+ my(@ibad,@ush_dev,@vsh_dev);
+
+ for (my($i)=0; $i<@dsh; $i++) {
+ die("$0: assertion failed") unless numberp($dsh[$i]);
+ my($bsgi) = int($dsh[$i] / $SHEAR_PREGRID_DZ);
+
+ $ush_dev[$i] = $ush[$i] - $ush_mu[$bsgi];
+ $vsh_dev[$i] = $vsh[$i] - $vsh_mu[$bsgi];
+
+ push(@ibad,$i) if ($ush_sig[$i] > 0 &&
+ (abs($ush_dev[$i]/$ush_sig[$bsgi]) > $max_shdev ||
+ abs($vsh_dev[$i]/$vsh_sig[$bsgi]) > $max_shdev));
+ } ## end of loop through shears
+
+ ## Look for internal glitches: a positive shear followed
+ ## immediately by a compensating negative shear, for
+ ## example. When one is found, flag the common velocity
+ ## sample, and untag the two shears by setting each ibad
+ ## to -1.
+ for (my($bi)=0; $bi<@ibad-1; $bi++) {
+ next unless ($ibad[$bi]+1 == $ibad[$bi+1]);
+
+ my($i) = $ibad[$bi];
+ my($bsgi) = int($dsh[$i] / $SHEAR_PREGRID_DZ);
+
+ if ($ush_sig[$bsgi] > 0 && $vsh_sig[$bsgi] > 0 &&
+ sqrt(($ush_dev[$i]+$ush_dev[$i+1])**2/$ush_sig[$bsgi]**2) < $max_shdev_sum &&
+ sqrt(($vsh_dev[$i]+$vsh_dev[$i+1])**2/$vsh_sig[$bsgi]**2) < $max_shdev_sum) {
+ $flag_count{$SHEAR_BIT}++;
+ $edit_flags[$ens][$sh_i1[$i]] |= $SHEAR_BIT;
+ $ibad[$bi] = $ibad[$bi+1] = -1;
+ }
+ } ## end of first loop through bad shears
+
+ ## Now flag all remaining velocities involved in the shears
+ ## listed by ibad.
+ for (my($bi)=0; $bi<@ibad; $bi++) {
+ next if ($ibad[$bi] < 0);
+ $flag_count{$SHEAR_BIT} += 2;
+ $edit_flags[$ens][$sh_i0[$ibad[$bi]]] |= $SHEAR_BIT;
+ $edit_flags[$ens][$sh_i1[$ibad[$bi]]] |= $SHEAR_BIT;
+ }
+}
+
+sub calc_shear($$$$)
+{
+ my($start,$end,$grid_dz,$edit_shear) = @_;
+ my($De) = $start<$end ? 1 : -1; # downcast: De = 1; upcast: De = -1
+
+ my(@ush_vals,@vsh_vals,@wsh_vals);
+
+ my($nvel,$nsh) = (0,0);
+ for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) { # loop over all ens from start to end
+
+ local(@sh_i0,@sh_i1);
+ local(@dsh,@ush,@vsh,@wsh);
+
+ uv_to_shear($ens);
+ if ($edit_shear) {
+ set_shear_flag($ens);
+ $nvel += uv_to_shear($ens);
+ $nsh += @dsh;
+ }
+
+ for (my($i)=0; $i<@dsh; $i++) { # save shears for binning calculations
+ my($gi) = int($dsh[$i] / $grid_dz);
+ push(@{$ush_vals[$gi]},$ush[$i]);
+ push(@{$vsh_vals[$gi]},$vsh[$i]);
+ push(@{$wsh_vals[$gi]},$wsh[$i]);
+ }
+ } # $ens loop
+
+ @ush_mu = @vsh_mu = @wsh_mu = ();
+ @ush_sig = @vsh_sig = @wsh_sig = ();
+
+ for (my($gi)=0; $gi<@ush_vals; $gi++) { # calc grid means & stddev
+ my($sum_ush,$sum_vsh,$sum_wsh);
+
+ $sh_n[$gi] = @{$ush_vals[$gi]};
+
+ for (my($vi)=0; $vi<$sh_n[$gi]; $vi++) {
+ $sum_ush += $ush_vals[$gi][$vi];
+ $sum_vsh += $vsh_vals[$gi][$vi];
+ $sum_wsh += $wsh_vals[$gi][$vi];
+ }
+ $ush_mu[$gi] = $sh_n[$gi] ? $sum_ush/$sh_n[$gi] : nan;
+ $vsh_mu[$gi] = $sh_n[$gi] ? $sum_vsh/$sh_n[$gi] : nan;
+ $wsh_mu[$gi] = $sh_n[$gi] ? $sum_wsh/$sh_n[$gi] : nan;
+ }
+
+ for (my($gi)=0; $gi<@ush_vals; $gi++) { # calc & grid stddevs
+ my($sumsq_ush,$sumsq_vsh,$sumsq_wsh);
+ for (my($vi)=0; $vi<$sh_n[$gi]; $vi++) {
+ $sumsq_ush += ($ush_vals[$gi][$vi] - $ush_mu[$gi])**2;
+ $sumsq_vsh += ($vsh_vals[$gi][$vi] - $vsh_mu[$gi])**2;
+ $sumsq_wsh += ($wsh_vals[$gi][$vi] - $wsh_mu[$gi])**2;
+ }
+ $ush_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_ush/($sh_n[$gi]-1)) : nan;
+ $vsh_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_vsh/($sh_n[$gi]-1)) : nan;
+ $wsh_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_wsh/($sh_n[$gi]-1)) : nan;
+ }
+
+ if ($edit_shear && $opt_d) {
+ print(STDERR "\n\t\t$nvel valid velocities");
+ print(STDERR "\n\t\t$nsh valid shears");
+ print(STDERR "\n\t\tflag counts:");
+ print(STDERR "\n\t\t\tBADVEL_BIT = $flag_count{$BADVEL_BIT}");
+ print(STDERR "\n\t\t\tERRVEL_BIT = $flag_count{$ERRVEL_BIT}");
+ print(STDERR "\n\t\t\tCORREL_BIT = $flag_count{$CORREL_BIT}");
+ print(STDERR "\n\t\t\tW_OUTLIER_BIT = $flag_count{$W_OUTLIER_BIT}");
+ print(STDERR "\n\t\t\tSHEAR_BIT = $flag_count{$SHEAR_BIT}");
+ print(STDERR "\n\t\t\tSIDELOBE_BIT = $flag_count{$SIDELOBE_BIT}");
+ print(STDERR "\n\t\t\tWAKE_BIT = $flag_count{$WAKE_BIT}");
+ print(STDERR "\n\t\t\tPPI_BIT = $flag_count{$PPI_BIT}");
+ printf(STDERR "\n\t\t\tW_CTD_BIT = $flag_count{$W_CTD_BIT} (%d ensembles)",
+ $flag_count{$W_CTD_BIT}/$LADCP{N_BINS});
+ printf(STDERR "\n\t\t\tTILT_BIT = $flag_count{$TILT_BIT} (%d ensembles)",
+ $flag_count{$TILT_BIT}/$LADCP{N_BINS});
+ printf(STDERR "\n\t\t\tDELTA_TILT_BIT = $flag_count{$DELTA_TILT_BIT} (%d ensembles)",
+ $flag_count{$DELTA_TILT_BIT}/$LADCP{N_BINS});
+ }
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/2011_07_10/LADCPproc.backscatter
@@ -0,0 +1,137 @@
+#======================================================================
+# L A D C P P R O C . B A C K S C A T T E R
+# doc: Wed Oct 20 13:02:27 2010
+# dlm: Thu Jul 7 08:55:45 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 96 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Oct 20, 2010: - created
+# Dec 10, 2010: - BUG: backscatter above sea surface made code bomb
+# when run with uplooker data
+# Jun 15, 2011: - added calculation of backscatter profiles from
+# subset of bins
+# Jul 7, 2011: - use $BEAM? vars to clarify code
+# - save SV values to use in BT code
+
+my($BEAM1) = 0;
+my($BEAM2) = 1;
+my($BEAM3) = 2;
+my($BEAM4) = 3;
+
+#----------------------------------------------------------------------
+# Volume Scattering Coefficient, following Deines (IEEE 1999)
+# NOTES:
+# - instrument specific! (300kHz Workhorse)
+# - no sound-speed correction applied
+# - R in bin center, instead of last quarter
+# - transmit power assumes 33V batteries
+#----------------------------------------------------------------------
+
+# NB:
+# - correction seems to work for a subset of bins (~bins 3-9 for
+# 2010 P403 station 46)
+# - this may imply that noise level depends on bin
+# - far bins are important for seabed detection, i.e. cannot simply
+# be discarded at this stage
+
+sub log10 {
+ my $n = shift;
+ return log($n)/log(10);
+}
+
+sub Sv($$$$$)
+{
+ my($temp,$PL,$Er,$R,$EA) = @_;
+ my($C) = -143; # RDI Workhorse monitor
+ my($Ldbm) = 10 * log10($PL);
+ my($PdbW) = 14.0;
+ my($alpha) = 0.069;
+ my($Kc) = 0.45;
+
+ return $C + 10*log10(($temp+273)*$R**2) - $Ldbm - $PdbW
+ + 2*$alpha*$R + $Kc*($EA-$Er);
+}
+
+sub mk_backscatter_profs($$)
+{
+ my($LADCP_start,$LADCP_end) = @_;
+
+ my(@Er) = (1e99,1e99,1e99,1e99); # echo intensity reference level
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ $Er[$BEAM1] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM1] < $Er[$BEAM1]);
+ $Er[$BEAM2] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM2]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM2] < $Er[$BEAM2]);
+ $Er[$BEAM3] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM3]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM3] < $Er[$BEAM3]);
+ $Er[$BEAM4] = $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM4]
+ if ($LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$LADCP{N_BINS}-1][$BEAM4] < $Er[$BEAM4]);
+ }
+ print(STDERR "\n\t\@per-beam noise levels = @Er") if ($opt_d);
+
+ for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
+ my($gi) = int(&depthOfBin($ens,$bin) / $GRID_DZ);
+ next if ($gi < 0);
+ my($range_to_bin) = &dzToBin($ens,$bin) / cos(rad($LADCP{BEAM_ANGLE}));
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM1] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM1],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM1]);
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM2] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM2],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM2]);
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM3] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM3],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM3]);
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin][$BEAM4] = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[$BEAM4],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][$BEAM4]);
+ my($Sv) = avg(@{$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]});
+
+ $sSv[$gi][$bin] += $Sv;
+ $nSv[$gi][$bin]++;
+
+ if ($bin>=$Svbin_start && $bin<=$Svbin_end) {
+ $sSv_prof[$gi] += $Sv;
+ $nSv_prof[$gi]++;
+ }
+ }
+ }
+}
+
+sub depthOfGI($) { return $_[0]*$GRID_DZ + $GRID_DZ/2; } # depth corresponding to particular grid index
+
+sub find_backscatter_seabed($)
+{
+ my($search_below) = @_;
+ my($mdgi) = int($search_below/$GRID_DZ); # grid index to begin search
+ my(@wdepth_gi); # water_depth indices
+
+ print(STDERR "\n\tseabed-max grid indices:") if ($opt_d);
+
+ for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { # find backscatter min/max below $search_below in each bin
+ my($min,$max,$gimax,$lastvalid) = (1e99,-1e99,-1,-1);
+ for (my($gi)=$mdgi; $gi<@nSv; $gi++) {
+ next unless ($nSv[$gi][$bin] > 0);
+ my($avg) = $sSv[$gi][$bin] / $nSv[$gi][$bin];
+ $lastvalid = $gi;
+ $min = $avg if ($avg < $min);
+ $max = $avg, $gimax = $gi if ($avg > $max);
+ }
+ if ($max-$min>10 && $gimax!=$lastvalid) { # ignore boundary maxima & scatter
+ printf(STDERR " %d",$gimax-$mdgi) if ($opt_d);
+ push(@wdepth_gi,$gimax);
+ }
+ }
+
+ return (depthOfGI(avg(@wdepth_gi)),stddev(@wdepth_gi)*$GRID_DZ);
+
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/2011_07_10/LADCPproc.bestLag
@@ -0,0 +1,189 @@
+#======================================================================
+# L A D C P P R O C . B E S T L A G
+# doc: Tue Sep 28 21:58:48 2010
+# dlm: Fri Jul 8 02:54:29 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 55 101 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# TODO:
+# - first lag is always(?) nan unless CTD is turned on with LADCP in water
+
+# HISTORY:
+# Sep 28, 2010: - created
+# Dec 9, 2010: - adapted to %CTD
+# Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags
+# Jan 5, 2011: - changed first guess from 80% down to 10% down
+# - added LADCP time lag to %PARAMs
+# - added support of -i
+# Jul 7, 2011: - added code to remove window-mean of w before lagging to
+# make it work in regions of crazy ocean w (IWISE 16007)
+
+sub interp_LADCP_w($$)
+{
+ my($elapsed,$ens) = @_;
+ my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) /
+ ($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} -
+ $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME});
+ unless (numberp($LADCP{ENSEMBLE}[$ens-1]->{W})) {
+ $nGaps++;
+ return $LADCP{ENSEMBLE}[$ens]->{W};
+ }
+ return $LADCP{ENSEMBLE}[$ens-1]->{W} +
+ $sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W});
+}
+
+sub bestLag($)
+{
+ my($ws) = @_; # window start index
+
+ my($best);
+ my($bestmad) = 9e99; # mean absolute deviation
+ for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) {
+ my($mCw,$mLw,$nw) = (0,0,0); # first calc means
+ for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
+ my($Li) = $Ci + $Llag;
+ next if ($Li<0 || $Li>=$opt_w);
+ next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
+ $mCw += $CTD{w}[$ws+$Ci];
+ $mLw += $LADCP_w[$ws+$Li];
+ $nw++;
+ }
+ next unless ($nw > 0);
+ $mCw /= $nw; $mLw /= $nw;
+
+ my($sad) = my($nad) = 0; # calc mad with means removed
+ for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
+ my($Li) = $Ci + $Llag;
+ next if ($Li<0 || $Li>=$opt_w);
+ next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
+ $sad += abs($CTD{w}[$ws+$Ci]-$mCw - ($LADCP_w[$ws+$Li]-$mLw));
+ $nad++;
+ }
+ if ($sad/$nad < $bestmad) {
+ $best = $Llag;
+ $bestmad = $sad/$nad;
+ }
+ }
+ return $best;
+}
+
+sub lagLADCP2CTD()
+{
+ #------------------------------------------------------------------------
+ # find 1st rec & ensemble >=10% down to max depth & make 1st guess at lag
+ #------------------------------------------------------------------------
+
+ my($first_guess_lag); # in units of CTD records
+
+ if (defined($opt_i)) {
+ $first_guess_lag = -$opt_i / $CTD{sampint};
+ } else {
+ my($CTD_10pct_down) = 0;
+ $CTD_10pct_down++ # "until" formulation allows for missing pressures
+ until ($CTD{press}[$CTD_10pct_down]-$CTD{press}[0] >= 0.1*($CTD{maxpress}-$CTD{press}[0]));
+
+ my($LADCP_10pct_down) = 0;
+ $LADCP_10pct_down++
+ while ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} < 0.1*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
+
+ $first_guess_lag = ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED_TIME} -
+ $CTD_10pct_down*$CTD{sampint}) / $CTD{sampint};
+
+ printf(STDERR "\n\t1st guess offset [CTD pressure, LADCP estimated depth] = %ds [%ddbar, %dm]\n",
+ $first_guess_lag*$CTD{sampint},$CTD{press}[$CTD_10pct_down],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH})
+ if ($opt_d);
+ croak("\n$0: cannot determine valid 1st guess offset ($LADCP_10pct_down,$CTD_10pct_down,$CTD{press}[$CTD_10pct_down])\n")
+ if ($first_guess_lag*$CTD{sampint} > 1200);
+ }
+
+ #------------------------------------------------------------------------------------
+ # Linearly interpolate LADCP time series onto a new grid with $CTD{sampint} resolution
+ # ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data
+ # chunks large
+ #------------------------------------------------------------------------------------
+
+ $nGaps = 0;
+
+ for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {
+ while ($r*$CTD{sampint} < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
+ $LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD{sampint},$ens)
+ unless ($first_guess_lag > $r);
+ $r++;
+ }
+ }
+
+ print(STDERR "\t$nGaps gaps in w timeseries")
+ if ($opt_d);
+
+ print(STDERR "\n");
+
+ #----------------------------------------------------------------------
+ # Calculate lags
+ #----------------------------------------------------------------------
+
+ printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]:",$opt_w);
+ $opt_w = int($opt_w / $CTD{sampint});
+
+ #---------------------------------------------------------------
+ # carry out opt_n lag correlations and keep tally of the results
+ #---------------------------------------------------------------
+ my(%nBest);
+ my($nLags) = 0;
+ my($lags) = '';
+ my($lastLag) = 9e99; my($nSame) = 1;
+ for (my($window)=0; $window<$opt_n; $window++) {
+ my($ws) = $window * int(@LADCP_w/$opt_n); # window start
+ $ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);
+ $bestLag = bestLag($ws);
+ if (defined($bestLag)) {
+ if (defined($lastLag) && $bestLag == $lastLag) {
+ $nSame++;
+ } else {
+ printf(STDERR "(x%d)",$nSame)
+ if ($nSame > 1);
+ printf(STDERR " %d",$bestLag*$CTD{sampint});
+ $nSame = 1;
+ $lastLag = $bestLag;
+ }
+ $lags .= sprintf(" %s",$bestLag*$CTD{sampint});
+ $nBest{$bestLag}++;
+ $nLags++;
+ } else {
+ if (!defined($lastLag)) {
+ $nSame++;
+ } else {
+ printf(STDERR "(x%d)",$nSame)
+ if ($nSame > 1);
+ printf(STDERR " nan");
+ $nSame = 1;
+ $lastLag = $bestLag;
+ }
+ }
+ }
+ printf(STDERR "(x%d)",$nSame) if ($nSame > 1);
+ &antsAddParams('time_lags',$lags);
+
+ #----------------------
+ # find most popular lag
+ #----------------------
+ my($best_lag);
+ foreach my $i (keys(%nBest)) {
+ $best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
+ }
+ croak("\n$0: cannot determine a valid lag\n")
+ unless ($nBest{$best_lag} > $opt_n/3);
+ print(STDERR "\n\n\t\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
+ if ($nBest{$best_lag} < $opt_n/2);
+
+ if ($nBest{$best_lag} == $nLags) {
+ printf(STDERR "\n\t\tunanimous lag = %ds\n",$best_lag*$CTD{sampint});
+ } else {
+ printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD{sampint});
+ }
+
+ &antsAddParams('LADCP_time_lag',($first_guess_lag + $best_lag) * $CTD{sampint});
+ return ($first_guess_lag + $best_lag) * $CTD{sampint};
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/2011_07_10/LADCPproc.defaults
@@ -0,0 +1,276 @@
+#======================================================================
+# L A D C P P R O C . D E F A U L T S
+# doc: Fri Sep 17 09:44:21 2010
+# dlm: Fri Jul 8 00:06:07 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 91 16 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+# default parameters for [mkShearProf]
+
+# NOTES:
+# - defaults are taken:
+# 1) from my current UH processing merge control files
+# 2) from the defaults set in Eric's [merge.c]
+# - the default version in the ANTS bin dir is always loaded
+# - if there is a version in the current processing directory it is loaded
+# afterwards
+# - per-station parameters can be chosed based on $LADCP_file and $CTD_file
+# - for additional notes, see [LADCPproc]
+
+# HISTORY:
+# Sep 17, 2010: - created
+# Dec 9, 2010: - added doc for ASCII CTD file support
+# Jun 15, 2011: - added Svbin_start, Svbin_end
+# Jul 7, 2011: - added BT processing parameters
+
+#----------------------------------------------------------------------
+# ASCII CTD file support
+#----------------------------------------------------------------------
+
+# By default, [LADCPproc] reads a standard seabird .cnv file. Alternatively,
+# the CTD data can be supplied in a plain ASCII file. To do so, the following
+# variables must be defined in the setup file (-s). Each variable defines
+# the field (column) number of the corresponding data field. Fields are
+# numbered beginning with 1.
+
+# $CTD_ASCII_sampfreq = 1; ## in Hz
+
+# $CTD_ASCII_press_field = 1;
+# $CTD_ASCII_temp_field = 2;
+# $CTD_ASCII_salin_field = 3;
+# $CTD_ASCII_lat_field = 4;
+# $CTD_ASCII_lon_field = 5;
+
+# The following variable is optional.
+
+# $CTD_ASCII_badval = 999;
+
+#----------------------------------------------------------------------
+# Backscatter Profile Parameters
+#----------------------------------------------------------------------
+
+# The output profile of volume-scattering coefficient is calculated
+# from a subset of the bins only. This is based on the observation,
+# from on a single P403 profile, that the volume-scattering coefficient
+# correction of Deines (IEEE 1999) yield similar results only for
+# bins 3-9.
+
+# The seabed-finding routine, on the other hand, uses acoustic
+# backscatter from all bins, as it is possible that the seabed is only
+# seen in the farthest bins.
+
+$Svbin_start = 3;
+$Svbin_end = 9;
+
+#----------------------------------------------------------------------
+# Bottom Track Parameters
+#----------------------------------------------------------------------
+
+# First bin to consider when looking for seabed to calculated post-
+# processed BT. NB: For consistency with EF's code, bin 1 is the
+# first bin.
+
+$BT_bin_start = 1;
+
+# CTD depth above seabed where BT search begins
+
+$BT_begin_search_above = 300;
+
+# Maximum allowed spread of bin number where max echo is found. Large values
+# imply sloping seabed and/or large instrument tilt. The default value of
+# 3 is probably too tight. The range spread calculation should be modified
+# to take instrument tilt into consideration.
+
+$BT_max_bin_spread = 3;
+
+# For tricky BT cases, the acoustic backscatter from the seabed can be plotted
+# and the bottom range set manually. This has been found to be very powerful
+# in case of 2011_IWISE yoyo profile 160, where the bottom was around 150m
+# from the lower turning points and where the method of Deines [1999] clearly
+# does not work.
+
+#$BT_min_depth = xxx;
+#$BT_max_depht = xxx;
+
+# Maximum difference between water depth and average distance of echo max.,
+# if $BT_min_depth and $BT_max_depth are not set. The stddev of the detected
+# water depth is added To this number.
+
+$BT_max_depth_error = 20;
+
+# Maximum allowed difference between reference-layer w and BT w. Note
+# that this must be relaxed in regions of significant near-bottom
+# vertical velocity, e.g. due to cross-slope currents.
+
+$BT_max_w_difference = 0.03;
+
+#----------------------------------------------------------------------
+# Shear Processing Parameters
+#----------------------------------------------------------------------
+
+## u_bin0, u_bin1, w_bin0, w_bin1: These set the first and
+## last bin indices used when integrating horizontal and
+## vertical velocity, respectively. The indices start
+## from 1 and are inclusive. These parameters do not
+## affect the calculation of shear. They should specify
+## a good reference layer (for the calculation of the
+## barotropic component) such as bins 1 to 4.
+
+$wbin_start = 1; ## These parameters start from 1, not zero
+$wbin_end = 5;
+$ubin_start = 1;
+$ubin_end = 5;
+
+## sh_bin0, sh_bin1: These are like the above except that
+## they affect the shear calculation only. They can
+## normally be left at the default range of 1 to 128
+## (disabled).
+
+$shbin_start = 1;
+$shbin_end = $LADCP{N_BINS};
+
+## w_ref_bin, w_dif: These control one of the editing
+## criteria recommended by Fischer and Visbeck, or rather
+## a modification of it. All velocity data are rejected
+## below the point at which the vertical velocity
+## estimate is w_dif larger or smaller than the estimate
+## in w_ref_bin. I have not found this criterion
+## helpful. 94/12/22: changed, so that only those points
+## where w actually deviates from the mean from 0 to w_ref_bin
+## are flagged.
+
+$w_ref_bin = 10;
+$w_dif = 0.05;
+
+## wake_hd_dif, wake_ang_min, min_wake_w, n_wake_bins: These control
+## editing based upon the direction and inclination of the package wake,
+## calculated from reference layer U,V, and W. wake_hd_dif
+## sets how close to the heading of the wake must be to
+## the heading of any beam for interference (in degrees).
+## wake_ang_min sets the minimum wake angle from the
+## vertical for interference (in degrees). min_wake_w sets
+## the minimum package speed required for wake interference
+## (although upward speed is measured negative, this parameter
+## is input as a positive value). All three criteria
+## must be satisfied for the wake flag to be set. n_wake_bins
+## determines how many bins from the top are removed from a
+## flagged profile, the default is 1 (top bin only).
+## Additionally, if the previous ensemble met all
+## interference criteria, the present ensemble
+## will be flagged even if the criteria are not met.
+## Default values are wake_hd_dif=0.0, wake_ang_min=90.0,
+## and min_wake_w=0.1, which results in no wake editing.
+
+$wake_hd_dif = 0.0; ## set wake editing defaults so that no wake editing
+$wake_ang_min = 90.0; ## occurs
+$min_wake_w = 0.1;
+$n_wake_bins = 1; ## wake editing default - top bin only
+
+## e_max is the maximum error velocity. An error velocity
+## greater than e_max will flag the other velocity
+## components. For the BB with 2-ping ensembles, 0.01
+## m/s looks about right for this parameter. It
+## knocks out the incorrect ambiguity glitches as well
+## as some smaller but still significant glitches.
+## Note that for the BB, the PG criteria operate on
+## percent 4-beam solutions, so PG combined with the
+## e_max provides a reliable filter against big
+## glitches.
+## For the NB with 18-ping ensembles, it is not yet
+## clear whether there is any benefit in using e_max
+## at all. To be effective, the value will have to be
+## small, something like 0.015 or 0.02. Note also
+## that for the NB, the pg array holds pg counting
+## both 4-beam and 3-beam solutions, so 3-beam
+## glitches will slip through the e_max net.
+## Default is e_max= 10.0, which effectively disables
+## it.
+
+$e_max = 0.1;
+
+## min_correlation (BB only) is the minimum correlation, in
+## counts, for each beam in a given bin.
+
+$min_cor = 70;
+
+## Shear editing: This requires 2 passes through the
+## database, one using the option
+## binned_shear_time_range:
+## to set the time range for a pass during which the
+## shear statistics will be calculated on a relatively
+## coarsely grid by binning rather than interpolation,
+## then using the usual option
+## time_range:
+## for the normal pass in which the results of the
+## previous pass can be used to flag bad velocity values
+## based on anomalous shears. The first pass can include
+## both up and down casts if desired, in which case it
+## would be followed by a second pass for the up and
+## downcasts separately. The binned shear statistics are
+## saved until they are explicitly recalculated with the
+## "binned_shear_time_range:" option.
+## The shear editing is controlled by these parameters:
+## shear_dev_max= x.x, where x.x is a floating point
+## number giving a threshold in standard
+## deviations. Any shear component deviating from
+## the binned mean by more than this times the
+## local standard deviation will raise a flag.
+## Suitable values for this parameter are probably
+## in the range 3-5. Less than 3 is likely to
+## start trimming too many valid samples, more than
+## 4 or 5 is likely to do nothing at all.
+## shear_sum_dev_max= x.x gives a threshold used in
+## detecting isolated bad velocity points. If the
+## sum of two successive bad shears (based on
+## shear_dev_max), divided by the
+## standard deviation, is LESS THAN this amount,
+## then the common velocity point is considered to
+## be an isolated glitch, and it alone is flagged.
+## Otherwise, both velocity samples contributing to
+## a flagged shear will be flagged. A reasonable
+## value for this parameter is probably around 1-2,
+## but I don't yet have enough experience to be
+## sure.
+## Warning: don't forget to reset the first_x=,
+## first_y=, and first_z= parameters after pass 1
+## (with binned_shear_time_range:). If you don't, and
+## the depths are not in the database, then the depths
+## will be completely wrong in pass 2 (with
+## time_range:). You may also want to use a different
+## output file name for pass 1, so that the binned
+## shear statistics file will not be overwritten
+## during pass 2.
+
+$max_shdev = 3.5; ## to disable, set to nan
+$max_shdev_sum = 1.5;
+
+## previous ping bottom bounce interference editing is controlled by:
+## clip_margin 0.0 turns off this editing; otherwise, it is the margin in
+## meters (on each side) by which the calculated range of the
+## interference is expanded before clipping. A reasonable value would
+## be something like 32 (2 depth bins).
+
+# NOTES:
+# - default value (90m) taken from comment in merge control files
+# - clipping is disabled if water_depth cannot be determined
+
+$clip_margin = 90; ## default: no previous ping bottom bounce editing
+$first_clip_bin = 1; ## default: apply previous ping clipping to all bins */
+
+# tilt editing as in Visbeck's code
+
+$max_tilt = 22; # max allowed angle from vertical
+$max_delta_tilt = 4; # max allowed ping-to-ping tilt difference
+
+# On DIMES US2 stations in Drake Passage it was found that the ambiguity velocity had been
+# set too low. It appears, however, that aliased velocities are easy to detect because aliasing
+# causes a large positive velocity to appear as a large negative one, and vice versa.
+# The following defines the maximum allowed discrepancy between the vertical velocity from CTD
+# pressure and the LADCP reference-layer w. This parameter should probably be set to something
+# similar as the ambiguity velocity.
+
+$w_max_err = 2.5;
+
+
+1;
new file mode 100644
--- /dev/null
+++ b/2011_07_10/LADCPproc.loadCTD
@@ -0,0 +1,183 @@
+#======================================================================
+# L A D C P P R O C . L O A D C T D
+# doc: Thu Dec 9 18:39:01 2010
+# dlm: Sat Jan 22 21:40:12 2011
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 174 0 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Dec 9, 2010: - exported from LADCPproc
+# - added support for ASCII files
+# Dec 16, 2010: - BUG cnv read did not work any more
+# Jan 10, 2011: - added code to skip ANTS header
+# Jan 22, 2011: - adapted to new -g
+
+sub readCTD_ASCII($$)
+{
+ my($fn,$dtaR) = @_;
+
+ croak("$fn: unknown pressure field\n") unless defined($CTD_ASCII_press_field);
+ croak("$fn: unknown temperature field\n") unless defined($CTD_ASCII_temp_field);
+ croak("$fn: unknown salinity field\n") unless defined($CTD_ASCII_salin_field);
+ unless (numberp($dtaR->{lat})) {
+ croak("$fn: unknown latitude field\n") unless defined($CTD_ASCII_lat_field);
+ croak("$fn: unknown longitude field\n") unless defined($CTD_ASCII_lon_field);
+ }
+
+ $CTD_ASCII_badval = 9e99 unless defined($CTD_ASCII_badval);
+
+ open(F,$fn) || croak("$fn: $!\n");
+ my($sumLat,$sumLon); my($nPos) = 0;
+ my($ds);
+ while (chomp($ds = <F>)) {
+ next if ($ds =~ /^#/);
+ my(@rec) = split('\s+',$ds);
+ push(@{$dtaR->{press}},($rec[$CTD_ASCII_press_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_press_field-1]);
+ push(@{$dtaR->{temp}}, ($rec[$CTD_ASCII_temp_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_temp_field-1]);
+ push(@{$dtaR->{salin}},($rec[$CTD_ASCII_salin_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_salin_field-1]);
+ unless (!defined($CTD_ASCII_lat_field) || $rec[$CTD_ASCII_lat_field-1] == $CTD_ASCII_badval) {
+ $nPos++;
+ $sumLat += $rec[$CTD_ASCII_lat_field-1];
+ $sumLon += $rec[$CTD_ASCII_lon_field-1];
+ }
+ }
+ close(F);
+
+ if ($nPos > 0) {
+ $dtaR->{lat} = $sumLat / $nPos;
+ $dtaR->{lon} = $sumLon / $nPos;
+ }
+
+ $dtaR->{sampint} = 1 / $CTD_ASCII_sampfreq;
+}
+
+sub readCTD_CNV($$)
+{
+ my($fn,$dtaR) = @_;
+ my($CTD_nrecs,$CTD_nfields,$pressF,$tempF,$salinF);
+ my($CTD_badval,$CTD_file_type);
+
+ open(F,$fn) || croak("$fn: $!\n");
+ while (1) { # parse header
+ my($hdr);
+ chomp($hdr = <F>);
+ $hdr =~ s/\r*$//;
+ croak("$0: unexpected EOF (format error)\n") unless defined($hdr);
+ last if ($hdr eq '*END*');
+
+ $CTD_nfields = $',next if ($hdr =~ /nquan = /); # Layout
+ $CTD_nrecs = $',next if ($hdr =~ /nvalues = /);
+ $pressF = $1,next if ($hdr =~ /name (\d+) = prDM:/);
+ if ($opt_2) {
+ $tempF = $1,next if ($hdr =~ /name (\d+) = t190C:/);
+ $salinF = $1,next if ($hdr =~ /name (\d+) = sal11:/);
+ } else {
+ $tempF = $1,next if ($hdr =~ /name (\d+) = t090C:/);
+ $salinF = $1,next if ($hdr =~ /name (\d+) = sal00:/);
+ }
+
+ &antsAddParams('start_time',$1),next # selected metadata
+ if ($hdr =~ /start_time = (.*)/);
+
+ &antsAddParams('station',$1),next
+ if ($hdr =~ /Station\s*:\s*(.*)/);
+ &antsAddParams('ship',$1),next
+ if ($hdr =~ /Ship\s*:\s*(.*)\s*$/);
+ &antsAddParams('cruise',$1),next
+ if ($hdr =~ /Cruise\s*:\s*(.*)\s*$/);
+ &antsAddParams('time',$1),next
+ if ($hdr =~ /Time\s*:\s*(.*)/);
+ &antsAddParams('date',$1),next
+ if ($hdr =~ /Date\s*:\s*(.*)/);
+
+ if ($hdr =~ /Latitude\s*[=:]\s*/) {
+ ($deg,$min,$NS) = split(/ /,$');
+ $dtaR->{lat} = $deg + $min/60;
+ $dtaR->{lat} *= -1 if ($NS eq 'S');
+ next;
+ }
+ if ($hdr =~ /Longitude\s*[=:]\s*/) {
+ ($deg,$min,$EW) = split(/ /,$');
+ $dtaR->{lon} = $deg + $min/60;
+ $dtaR->{lon} *= -1 if ($EW eq 'W');
+ next;
+ }
+
+ if ($hdr =~ /interval = seconds: /) {
+ $dtaR->{sampint} = 1*$';
+ next;
+ }
+
+ $CTD_badval = $',next
+ if ($hdr =~ /bad_flag = /);
+ $CTD_file_type = $',next
+ if ($hdr =~ /file_type = /);
+ }
+
+ croak("$CTD_file: cannot determine CTD file layout\n")
+ unless ($CTD_nfields && $CTD_nrecs);
+ croak("$CTD_file: cannot determine missing value\n")
+ unless defined($CTD_badval);
+ croak("$CTD_file: not a CTD time series file\n")
+ unless ($dtaR->{sampint});
+ croak("$CTD_file: no pressure field\n")
+ unless defined($pressF);
+ croak("$CTD_file: no suitable temperature field\n")
+ unless defined($tempF);
+ croak("$CTD_file: no suitable salinity field\n")
+ unless defined($salinF);
+
+ if ($CTD_file_type eq 'ascii') {
+ while (1) {
+ last unless (my(@rec) = &antsFileIn(F));
+ push(@{$dtaR->{press}},($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
+ push(@{$dtaR->{temp}}, ($rec[$tempF] == $CTD_badval) ? nan : $rec[$tempF]);
+ push(@{$dtaR->{salin}},($rec[$salinF] == $CTD_badval) ? nan : $rec[$salinF]);
+ }
+ } elsif ($CTD_file_type eq 'binary') {
+
+ my($fbits) = 8 * length(pack('f',0));
+ croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
+ unless ($fbits == 32);
+
+ croak("$fn: can't read binary data\n")
+ unless (read(F,my($dta),4*$CTD_nfields*$CTD_nrecs) == 4*$CTD_nfields*$CTD_nrecs);
+ print(STDERR "$fn: WARNING: extraneous data at EOF\n") unless eof(F);
+
+ $dta = pack('V*',unpack('N*',$dta)) # big-endian CPU
+ if (unpack('h*', pack('s', 1)) =~ /01/); # c.f. perlport(1)
+
+ my(@dta) = unpack("f*",$dta);
+
+ for (my($r)=0; $r<$CTD_nrecs; $r++) {
+ push(@{$dtaR->{press}},($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
+ push(@{$dtaR->{temp}}, ($dta[$r*$CTD_nfields+$tempF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$tempF]);
+ push(@{$dtaR->{salin}},($dta[$r*$CTD_nfields+$salinF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$salinF]);
+ }
+ } else {
+ croak("$fn: unknown CTD file type $CTD_file_type\n");
+ }
+ close(F);
+}
+
+sub readCTD($$)
+{
+ my($fn,$dtaR) = @_;
+
+ if (defined($CTD_ASCII_sampfreq)) {
+ readCTD_ASCII($fn,$dtaR);
+ } else {
+ readCTD_CNV($fn,$dtaR);
+ }
+
+ croak("$0: unknown latitude\n") unless defined($dtaR->{lat});
+ &antsAddParams('lat',$dtaR->{lat});
+ croak("$0: unknown longitude\n") unless defined($dtaR->{lon});
+ &antsAddParams('lon',$dtaR->{lon});
+
+ &antsAddParams('CTD_sampfreq',1/$dtaR->{sampint});
+ &antsAddParams('ITS',$P{ITS} = 90);
+}
+
+1;
new file mode 100644
--- /dev/null
+++ b/2011_07_10/README
@@ -0,0 +1,14 @@
+======================================================================
+ R E A D M E
+ doc: Wed Oct 13 12:29:50 2010
+ dlm: Sun Jul 10 16:29:37 2011
+ (c) 2010 A.M. Thurnherr
+ uE-Info: 9 68 NIL 0 0 72 0 2 4 NIL ofnI
+======================================================================
+
+I froze this version before adding output for shear Visbeck diagram.
+
+=Differences with UH Software=
+
+- much simpler shear gridding
+- no PPI editing
new file mode 100644
--- /dev/null
+++ b/2011_07_10/TODO
@@ -0,0 +1,58 @@
+======================================================================
+ T O D O
+ doc: Thu Oct 14 09:15:02 2010
+ dlm: Wed Jun 15 15:16:20 2011
+ (c) 2010 A.M. Thurnherr
+ uE-Info: 40 0 NIL 0 0 72 3 2 4 NIL ofnI
+======================================================================
+
+=LADCPintsh=
+
+- determine what constitutes a gap
+- look into cross-gap shear integration
+
+
+=LADCPproc=
+
+- allow setting of GRID_DZ (-o) in ProcessingParams
+
+- add warnings:
+ - wide band
+ - CTD_W_BIT (ambiguity velocity)
+ - bad time match
+ - seabed not found
+ - inconsistent water depth
+
+- add diagnostic plots
+
+- output shear vals to make histograms
+
+- try median instead of mean for shear binning
+
+- add pitch/roll to BT
+
+- check interaction between 3-beam solutions and correlation limit
+
+- implement TILT_BIT (others?)
+
+- calculate real acoustic backscatter profile
+ - currently, profiles are arbitarily calculated from bins 3-9
+
+- clean up LADCPproc.UHcode:
+ - remove weirdnesses
+ - avoid multiple calculation of &depthOfBin()
+ - improve "staging":
+ 1) basic velocity editing (BADVEL, ERRVEL, CORREL) NB: before
+ beam-to-Earth transformation!!!
+ 2) ref-lr w (also handle W_BIT editing?)
+ 3) integrate w to get z_approx
+ 4) match LADCP to CTD time series
+ 5) make accoustic backscatter profile
+ 6) find seabed
+ 7) sidelobe & PPI editing
+ 8) ref-lr u,v,w
+ 9) WAKE editing
+ 10) SHEAR editing
+ 11) ref-lr u,v,w
+
+- re-add time-series calculation
new file mode 100755
--- /dev/null
+++ b/2011_07_10/splitCNV.olderVersion
@@ -0,0 +1,48 @@
+#!/usr/bin/perl
+#======================================================================
+# S P L I T C N V
+# doc: Mon Oct 25 22:24:21 2010
+# dlm: Fri Oct 29 20:04:42 2010
+# (c) 2010 A.M. Thurnherr
+# uE-Info: 16 29 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# split SeaBird CNV file using scan values found with yoyo
+# e.g. splitCNV P403_015_1sec.cnv `importCNV -s P403_015_1sec.cnv | yoyo -QFscan -ut`
+
+# HISTORY:
+# Oct 25, 2010: - created
+# Oct 26, 2010: - BUG: *.00 file always only contained the 1st record
+# Oct 29, 2010: - cosmetics
+
+# scan number must be the first field of each record
+
+die("Usage: $0 <CNV-file> <val> [...]\n")
+ unless (@ARGV >= 2 && -f $ARGV[0]);
+
+chomp($basename = `basename $ARGV[0] .cnv`);
+
+$state = 0; # init finite state machine
+$next = 0; # next extension number
+$trg = 2; # NB: yoyo -ut includes first & last scans
+
+open(IN,$ARGV[0]);
+while (<IN>) {
+ if ($state == 0) { # state 0: reading header
+ push(@HDR,$_);
+ $state = 1 if /^\*END\*/;
+ next;
+ }
+ if ($state == 1) { # state 1: begin new file
+ close(OUT);
+ printf(STDERR "writing $basename.%02d\n",$next);
+ open(OUT,sprintf(">$basename.%02d",$next++));
+ foreach my $h (@HDR) { print(OUT $h); }
+ $state = 2;
+ }
+ if ($state == 2) { # state 2: copy data until target scan
+ print(OUT);
+ my(@f) = split;
+ $trg++,$state=1 if ($f[0] == $ARGV[$trg]);
+ }
+}
--- a/LADCPproc
+++ b/LADCPproc
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P P R O C
# doc: Thu Sep 16 20:36:10 2010
-# dlm: Fri Jul 8 03:33:21 2011
+# dlm: Wed Jul 13 11:43:18 2011
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 491 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 422 55 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'process LADCP data to get shear, time series';
@@ -42,6 +42,7 @@
# acoustic backscatter depth-time-series
# - disabled seabed and acoustic backscatter code when not required (e.g. UL)
# - made non-diagnostic output terser
+# Jul 11, 2011: - changed default output to .tds and added -p)rofile option
($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -60,7 +61,7 @@
require "$PERL_TOOLS/RDI_Utils.pl";
$antsParseHeader = 0;
-&antsUsage('24a:b:c:df:g:i:kl:n:o:ps:t:w:',2,
+&antsUsage('24a:b:c:df:g:i:kl:n:o:p:s:t:w:',2,
'[use -2)dary CTD sensor pair]',
'[require -4)-beam LADCP solutions]',
'[-s)etup <file>] [-g)ps <lat,lon>]',
@@ -69,9 +70,9 @@
'[-o)utput grid <resolution[5m]>]',
'[-i)nitial LADCP time lag <guestimate>]',
'[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
- '[-d)iagnostic output]',
- 'output: [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
- ' [-a)coustic backscatter <depth-time-series file] [bottom-trac-k) profs]',
+ '[-d)iagnostic screen output]',
+ 'output: [shear-p)rofile <file>] [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
+ ' [-a)coustic backscatter <dts-file] [bottom-trac-k) profs]',
'<RDI file> <SeaBird file>');
$RDI_Coords::minValidVels = 4 if ($opt_4);
@@ -85,7 +86,7 @@
&antsCardOpt(\$opt_n,20);
&antsFileOpt($opt_s);
&antsFloatOpt($opt_i);
-&antsCardOpt($opt_o);
+&antsCardOpt(\$opt_o,5);
if (defined($opt_g)) {
($CTD{lat},$CTD{lon}) = split(',',$opt_g);
@@ -142,7 +143,7 @@
$LADCP{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'uplooker' : 'downlooker');
$SHEAR_PREGRID_DZ = 20;
-$GRID_DZ = defined($opt_o) ? $opt_o : 5;
+$GRID_DZ = $opt_o;
my($year) = substr($LADCP{ENSEMBLE}[0]->{DATE},6,4);
my($month) = substr($LADCP{ENSEMBLE}[0]->{DATE},0,2);
@@ -397,10 +398,70 @@
}
#----------------------------------------------------------------------
-# Step 8: Edit Data
+# Step 8: Edit Data & also produce depth-times-series output vial callback
#----------------------------------------------------------------------
-print(STDERR "Calculating shear profiles...");
+print(STDERR "Calculating shear profiles & producing time-depth-series (.tds) output...");
+
+ #--------------------------------------------------------------------------------
+ # callback routine to output .tds data, called once each for down-/upcasts after
+ # shear has been gridded. Note that only the nominal (bin center) depth is known
+ # for each sample in this version.
+
+ sub outTDseries($)
+ {
+ my($downcast) = @_; # also use local @[uvw]sh_vals[][]
+
+ if ($downcast) {
+ my($mingi);
+ for ($mingi=0; $mingi<@ush_vals; $mingi++) {
+ last if @{$ush_vals[$mingi]};
+ }
+ &antsAddParams('min_ens',$LADCP_start,'min_elapsed',$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME},
+ 'max_ens',$LADCP_end,'max_elapsed',$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
+ 'min_depth',depthOfGI($mingi),'max_depth',depthOfGI($#ens_vals),'foo',$#ens_vals);
+ for (my($gi)=0; $gi<@ush_vals; $gi++) {
+ for (my($i)=0; $i<@{$ush_vals[$gi]}; $i++) {
+ &antsOut($ens_vals[$gi][$i],$LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME},
+ $downcast,depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ }
+ }
+ } else {
+ for (my($gi)=$#ush_vals; $gi>=0; $gi--) {
+ for (my($i)=0; $i<@{$ush_vals[$gi]}; $i++) {
+ &antsOut($ens_vals[$gi][$i],$LADCP{ENSEMBLE}[$ens_vals[$gi][$i]]->{ELAPSED_TIME},
+ $downcast,depthOfGI($gi),$ush_vals[$gi][$i],$vsh_vals[$gi][$i],$wsh_vals[$gi][$i]);
+ }
+ }
+ }
+ } # sub outDSseries
+
+ #--------------------------------------------------------------------------------
+
+@antsNewLayout = ('ens','elapsed','downcast','depth','u_z','v_z','w_z');
+$commonParams = $antsCurParams;
+
+&antsAddParams('ubin_start',$ubin_start,'ubin_end',$ubin_end, # record processing params
+ 'wbin_start',$wbin_start,'wbin_end',$wbin_end,
+ 'shbin_start',$shbin_start,'shbin_end',$shbin_end,
+ 'w_ref_bin',$w_ref_bin,'w_dif',$w_dif,
+ 'wake_hd_dif',$wake_hd_dif,'wake_ang_min',$wake_ang_min,
+ 'min_wake_w',$min_wake_w,'n_wake_bins',$n_wake_bins,
+ 'e_max',$e_max,'min_cor',$min_cor,
+ 'max_shdev',$max_shdev,'max_shdev_sum',$max_shdev_sum,
+ 'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth),
+ 'min_hab',round($min_hab),'PPI_editing_enabled',$PPI_editing_enabled,
+ 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin,
+ 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end,
+ 'BT_bin_start',$BT_bin_start,'BT_bin_search_above',$BT_bin_search_above,
+ 'BT_max_bin_spread',$BT_max_bin_spread,'BT_max_w_difference',$BT_max_w_difference,
+);
+if (defined($BT_min_depth)) {
+ &antsAddParams('BT_min_depth',$BT_min_depth,'BT_max_depth',$BT_max_depth);
+} else {
+ &antsAddParams('BT_max_depth_error',$BT_max_depth_error);
+}
+$fullParams = $antsCurParams;
$LADCP_start = 1 if ($LADCP_start == 0); # ensure that there is previous ensemble
@@ -462,57 +523,44 @@
}
#----------------------------------------------------------------------
-# Step 10: Write Output Files
+# Step 10: Write Output Profiles if requested
#----------------------------------------------------------------------
-print(STDERR "Writing shear profiles...");
+if (defined($opt_p)) {
+
+ print(STDERR "Writing shear profiles...");
+
+ @antsNewLayout = ('depth','dc_nshear','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
+ 'uc_nshear','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
+ 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n');
+
+ &antsOut('EOF');
+ close(STDOUT);
+ open(STDOUT,">$opt_p") || croak("$opt_p: $!\n");
+
+ $antsCurParams = $fullParams;
-@antsNewLayout = ('depth','dc_nshear','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
- 'uc_nshear','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
- 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n');
-
-$commonParams = $antsCurParams;
-&antsAddParams('ubin_start',$ubin_start,'ubin_end',$ubin_end, # record processing params
- 'wbin_start',$wbin_start,'wbin_end',$wbin_end,
- 'shbin_start',$shbin_start,'shbin_end',$shbin_end,
- 'w_ref_bin',$w_ref_bin,'w_dif',$w_dif,
- 'wake_hd_dif',$wake_hd_dif,'wake_ang_min',$wake_ang_min,
- 'min_wake_w',$min_wake_w,'n_wake_bins',$n_wake_bins,
- 'e_max',$e_max,'min_cor',$min_cor,
- 'max_shdev',$max_shdev,'max_shdev_sum',$max_shdev_sum,
- 'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth),
- 'min_hab',round($min_hab),
- 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin,
- 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end,
- 'BT_bin_start',$BT_bin_start,'BT_bin_search_above',$BT_bin_search_above,
- 'BT_max_bin_spread',$BT_max_bin_spread,'BT_max_w_difference',$BT_max_w_difference,
-);
-if (defined($BT_min_depth)) {
- &antsAddParams('BT_min_depth',$BT_min_depth,'BT_max_depth',$BT_max_depth);
-} else {
- &antsAddParams('BT_max_depth_error',$BT_max_depth_error);
-}
-
-for (my($gi)=0; $gi<@ush_mu; $gi++) {
- &antsOut(depthOfGI($gi), # depth in center of bin
- numberp($dc_sh_n[$gi])?$dc_sh_n[$gi]:0, # downcast
- $dc_ush_mu[$gi],$dc_ush_sig[$gi],
- $dc_vsh_mu[$gi],$dc_vsh_sig[$gi],
- $dc_wsh_mu[$gi],$dc_wsh_sig[$gi],
- numberp($uc_sh_n[$gi])?$uc_sh_n[$gi]:0, # upcast
- $uc_ush_mu[$gi],$uc_ush_sig[$gi],
- $uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
- $uc_wsh_mu[$gi],$uc_wsh_sig[$gi],
- $sh_n[$gi], # combined
- $ush_mu[$gi],$ush_sig[$gi],
- $vsh_mu[$gi],$vsh_sig[$gi],
- $wsh_mu[$gi],$wsh_sig[$gi],
- $nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan,
- $nSv_prof[$gi],
- );
-}
-
-print(STDERR "\n");
+ for (my($gi)=0; $gi<@ush_mu; $gi++) {
+ &antsOut(depthOfGI($gi), # depth in center of bin
+ numberp($dc_sh_n[$gi])?$dc_sh_n[$gi]:0, # downcast
+ $dc_ush_mu[$gi],$dc_ush_sig[$gi],
+ $dc_vsh_mu[$gi],$dc_vsh_sig[$gi],
+ $dc_wsh_mu[$gi],$dc_wsh_sig[$gi],
+ numberp($uc_sh_n[$gi])?$uc_sh_n[$gi]:0, # upcast
+ $uc_ush_mu[$gi],$uc_ush_sig[$gi],
+ $uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
+ $uc_wsh_mu[$gi],$uc_wsh_sig[$gi],
+ $sh_n[$gi], # combined
+ $ush_mu[$gi],$ush_sig[$gi],
+ $vsh_mu[$gi],$vsh_sig[$gi],
+ $wsh_mu[$gi],$wsh_sig[$gi],
+ $nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan,
+ $nSv_prof[$gi],
+ );
+ }
+
+ print(STDERR "\n");
+} # if -p
#---------------------------------------
# Acoustic backscatter depth-time-series
--- a/LADCPproc.UHcode
+++ b/LADCPproc.UHcode
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . U H C O D E
# doc: Fri Sep 17 20:27:53 2010
-# dlm: Mon Jan 10 13:23:13 2011
+# dlm: Tue Jul 12 21:59:07 2011
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 295 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 392 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# PERLified functions from Eric's [merge.c]; with mods
@@ -32,6 +32,8 @@
# Oct 25, 2010: - added W_CTD_BIT, renamed W_BIT to W_OUTLIER_BIT
# Oct 26, 2010: - added TILT_BIT
# Dec 10, 2010: - modified assertion to allow processing of UL data
+# Jul 10, 2011: - added outTDseries() call
+# Jul 12, 2011: - replaced -p by $PPI_editing_enabled flag
#======================================================================
# VELOCITY EDITING
@@ -265,7 +267,7 @@
set_wake_flag($ens,$De);
set_misc_flags($ens,$De);
set_PPI_flags($ens,$De)
- if $opt_p && ($clip_margin > 0); # PPI editing is off by default
+ if $PPI_editing_enabled && ($clip_margin > 0); # PPI editing is off by default
}
}
@@ -304,6 +306,7 @@
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$V]) / ($d1-$d0);
$wsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$W] -
$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$W]) / ($d1-$d0);
+ $ens[$i] = $ens;
}
return $nvel;
@@ -362,13 +365,13 @@
my($start,$end,$grid_dz,$edit_shear) = @_;
my($De) = $start<$end ? 1 : -1; # downcast: De = 1; upcast: De = -1
- my(@ush_vals,@vsh_vals,@wsh_vals);
+ local(@ush_vals,@vsh_vals,@wsh_vals,@ens_vals);
my($nvel,$nsh) = (0,0);
for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) { # loop over all ens from start to end
local(@sh_i0,@sh_i1);
- local(@dsh,@ush,@vsh,@wsh);
+ local(@dsh,@ush,@vsh,@wsh,@ens);
uv_to_shear($ens);
if ($edit_shear) {
@@ -382,9 +385,12 @@
push(@{$ush_vals[$gi]},$ush[$i]);
push(@{$vsh_vals[$gi]},$vsh[$i]);
push(@{$wsh_vals[$gi]},$wsh[$i]);
+ push(@{$ens_vals[$gi]},$ens[$i]);
}
} # $ens loop
+ outTDseries($De==1) if ($edit_shear); # output depth-time time series
+
@ush_mu = @vsh_mu = @wsh_mu = ();
@ush_sig = @vsh_sig = @wsh_sig = ();
--- a/LADCPproc.defaults
+++ b/LADCPproc.defaults
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . D E F A U L T S
# doc: Fri Sep 17 09:44:21 2010
-# dlm: Fri Jul 8 00:06:07 2011
+# dlm: Tue Jul 12 21:11:56 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 91 16 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 26 46 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# default parameters for [mkShearProf]
@@ -23,6 +23,7 @@
# Dec 9, 2010: - added doc for ASCII CTD file support
# Jun 15, 2011: - added Svbin_start, Svbin_end
# Jul 7, 2011: - added BT processing parameters
+# Jul 12, 2011: - added $PPI_editing_enabled
#----------------------------------------------------------------------
# ASCII CTD file support
@@ -251,6 +252,8 @@
## interference is expanded before clipping. A reasonable value would
## be something like 32 (2 depth bins).
+undef($PPI_editing_enabled); # set to 1 to enable PPI editing
+
# NOTES:
# - default value (90m) taken from comment in merge control files
# - clipping is disabled if water_depth cannot be determined