new file mode 100644
--- /dev/null
+++ b/.hgignore
@@ -0,0 +1,11 @@
+#======================================================================
+# . H G R C
+# doc: Tue May 15 18:38:20 2012
+# dlm: Tue May 15 18:38:20 2012
+# (c) 2012 A.M. Thurnherr
+# uE-Info: 3 11 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+syntax: glob
+
+non_public/
deleted file mode 100755
--- a/2011_07_10/LADCP2CTDmatch
+++ /dev/null
@@ -1,333 +0,0 @@
-#!/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();
-
deleted file mode 100755
--- a/2011_07_10/LADCPintsh
+++ /dev/null
@@ -1,446 +0,0 @@
-#!/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();
deleted file mode 100755
--- a/2011_07_10/LADCPproc
+++ /dev/null
@@ -1,625 +0,0 @@
-#!/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();
-
deleted file mode 100644
--- a/2011_07_10/LADCPproc.BT
+++ /dev/null
@@ -1,187 +0,0 @@
-#======================================================================
-# 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;
deleted file mode 100644
--- a/2011_07_10/LADCPproc.UHcode
+++ /dev/null
@@ -1,439 +0,0 @@
-#======================================================================
-# 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;
deleted file mode 100644
--- a/2011_07_10/LADCPproc.backscatter
+++ /dev/null
@@ -1,137 +0,0 @@
-#======================================================================
-# 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;
deleted file mode 100644
--- a/2011_07_10/LADCPproc.bestLag
+++ /dev/null
@@ -1,189 +0,0 @@
-#======================================================================
-# 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;
deleted file mode 100644
--- a/2011_07_10/LADCPproc.defaults
+++ /dev/null
@@ -1,276 +0,0 @@
-#======================================================================
-# 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;
deleted file mode 100644
--- a/2011_07_10/LADCPproc.loadCTD
+++ /dev/null
@@ -1,183 +0,0 @@
-#======================================================================
-# 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;
deleted file mode 100644
--- a/2011_07_10/README
+++ /dev/null
@@ -1,14 +0,0 @@
-======================================================================
- 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
deleted file mode 100644
--- a/2011_07_10/TODO
+++ /dev/null
@@ -1,58 +0,0 @@
-======================================================================
- 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
deleted file mode 100755
--- a/2011_07_10/splitCNV.olderVersion
+++ /dev/null
@@ -1,48 +0,0 @@
-#!/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]);
- }
-}
deleted file mode 120000
--- a/Documentation/Sv_comparison.eps
+++ /dev/null
@@ -1,1 +0,0 @@
-/Data/LADCP/Data/2011_IWISE/shearmethod/Sv_methods/Sv_comparison.eps
\ No newline at end of file
new file mode 100644
--- /dev/null
+++ b/HISTORY
@@ -0,0 +1,13 @@
+======================================================================
+ H I S T O R Y
+ doc: Tue May 15 18:04:39 2012
+ dlm: Tue May 15 19:20:18 2012
+ (c) 2012 A.M. Thurnherr
+ uE-Info: 10 12 NIL 0 0 72 3 2 8 NIL ofnI
+======================================================================
+
+May 15, 2012:
+ - V1.0beta [.hg/hgrc]
+ - began history
+ - uploaded current version to server for use with first version
+ of re-implemented shear method
--- 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: Tue Apr 17 18:14:06 2012
+# dlm: Wed Apr 18 05:23:48 2012
# (c) 2010 A.M. Thurnherr
-# uE-Info: 69 93 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 70 45 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'process LADCP data to get shear, time series';
@@ -67,6 +67,7 @@
# - BUG: code had assumed, but not ensured, that first CTD scan is valid
# Apr 13, 2012: - removed -s argument from dependencies
# Apr 17, 2012: - BAD BUG: magdec code call was bad and did not return correct value. ever.
+# Apr 18, 2012: - replaced Sv.n by Sv.nsamp
($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -593,7 +594,7 @@
@antsNewLayout = ('depth','dc_elapsed','dc_nsamp','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
'uc_elapsed','uc_nsamp','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
- 'elapsed','nsamp','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n');
+ 'elapsed','nsamp','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.nsamp');
&antsOut('EOF');
close(STDOUT);
new file mode 100644
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,16 @@
+#======================================================================
+# M A K E F I L E
+# doc: Tue May 15 18:12:31 2012
+# dlm: Tue May 15 19:19:34 2012
+# (c) 2012 A.M. Thurnherr
+# uE-Info: 16 29 NIL 0 0 72 0 2 4 NIL ofnI
+#======================================================================
+
+.PHONY: version
+version:
+ @sed -n '/^description =/s/description = //p' .hg/hgrc
+
+.PHONY: publish
+publish:
+ cd ..; \
+ scp -Cr LADCP_shearmethod miles:public_hg
deleted file mode 100644
--- a/README
+++ /dev/null
@@ -1,12 +0,0 @@
-======================================================================
- R E A D M E
- doc: Wed Oct 13 12:29:50 2010
- dlm: Wed Oct 13 12:30:04 2010
- (c) 2010 A.M. Thurnherr
- uE-Info: 13 0 NIL 0 0 72 0 2 4 NIL ofnI
-======================================================================
-
-=Differences with UH Software=
-
-- much simpler shear gridding
-- no PPI editing
new file mode 100644
--- /dev/null
+++ b/README.Install
@@ -0,0 +1,21 @@
+======================================================================
+ R E A D M E . I N S T A L L
+ doc: Tue May 15 18:42:56 2012
+ dlm: Tue May 15 19:15:36 2012
+ (c) 2012 A.M. Thurnherr
+ uE-Info: 16 37 NIL 0 0 72 3 2 4 NIL ofnI
+======================================================================
+
+- requires ADCP_tools available via link from
+ http://www.ldeo.columbia.edu/LADCP
+
+- to load output files into Matlab requires Matlab_tools available via
+ link from http://www.ldeo.columbia.edu/LADCP
+
+- requires Eric Firing's geomag code available from
+ http://currents.soest.hawaii.edu/hg
+
+- runs with perl 5.12.4 or later, but older versions may work as well
+
+- to set up, add both ADCP_tools and the current directory to path
+
new file mode 100644
--- /dev/null
+++ b/README.OutputFiles
@@ -0,0 +1,11 @@
+======================================================================
+ R E A D M E . O U T P U T F I L E S
+ doc: Tue May 15 19:17:50 2012
+ dlm: Tue May 15 19:17:50 2012
+ (c) 2012 A.M. Thurnherr
+ uE-Info: 2 0 NIL 0 0 72 3 2 8 NIL ofnI
+======================================================================
+
+All output files are in a proprietary undocumented ASCII format. To
+load the files into Matlab, use loadANTS.m (see [README.Install]).
+
new file mode 100644
--- /dev/null
+++ b/README.Run
@@ -0,0 +1,56 @@
+======================================================================
+ R E A D M E . R U N
+ doc: Tue May 15 18:49:00 2012
+ dlm: Tue May 15 19:16:34 2012
+ (c) 2012 A.M. Thurnherr
+ uE-Info: 52 0 NIL 0 0 72 3 2 8 NIL ofnI
+======================================================================
+
+NB: DO NOT RUN FROM CURRENT DIRECTORY. (Actually, you can but will have
+ to prefix each command with a ./)
+
+
+=DATA REQUIREMENTS=
+
+- raw RDI BB ADCP files from a downlooker and/or an uplooker
+- SeaBird cnv file (ASCII or binary) with lat/lon info in the header and
+ with the following fields:
+ timeS, prDM, t090C and/or t190C, sal00 and/or sal11
+- arbitrary ASCII CTD files are also acceptable but there is no
+ documentation yet
+
+
+=CALCULATE LADCP SHEAR PROFILE=
+
+assume files:
+ 001DL000.000 downlooker ADCP file
+ 001UL000.000 uplooker ADCP file
+ 001.cnv CTD file
+
+LADCPproc -p 001DL.sh -b 001.BT 001DL000.000 001.cnv > /dev/null
+ creates two files, 001DL.sh (shear profiles) and 001.BT
+ (bottom-track data)
+
+LADCPproc -p 001UL.sh 001UL000.000 001.cnv > /dev/null
+ creates one file, 001UL.sh (shear profiles)
+
+
+
+=CALCULATE LADCP VELOCITY PROFILE=
+
+LADCPintsh 001DL.sh > 001DL.bc
+ creates baroclinic velocity profile from DL data
+
+LADCPintsh -b 001.BT 001DL.sh > 001DL.vel
+ creates BT-referenced absolute velocity profile from DL data
+
+LADCPintsh -b 001.BT -u 001UL.sh 001DL.sh > 001.vel
+ creates BT-referenced absolute velocity profile from combined
+ DL/UL data
+
+
+=QUALITY CHECKS=
+
+- compare to LDEO solution
+- compare UL vs DL solutions, if available
+- compare dc vs uc u, v, and, in particular w profiles
deleted file mode 100755
--- a/splitCNV.olderVersion
+++ /dev/null
@@ -1,48 +0,0 @@
-#!/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]);
- }
-}