--- a/LADCP_VKE Mon Apr 20 14:04:04 2015 +0000
+++ b/LADCP_VKE Sun Jul 26 20:04:48 2015 +0000
@@ -2,13 +2,13 @@
#======================================================================
# L A D C P _ V K E
# doc: Tue Oct 14 11:05:16 2014
-# dlm: Thu Apr 16 10:29:45 2015
+# dlm: Wed Jun 17 11:40:44 2015
# (c) 2012 A.M. Thurnherr
-# uE-Info: 188 54 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 399 18 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';
-$antsMinLibVersion = 6.0;
+$antsMinLibVersion = 6.1;
# HISTORY:
# Oct 14, 2014: - created from [LADCPfs]
@@ -18,38 +18,127 @@
# Nov 7, 2014: - changed parameterization constant $c to 0.0215
# Apr 16, 2015: - disabled output activation unless ANTS tools are available
# - removed superfluous $ANTSBIN definition
-
-# NOTES:
-# - requires power densities
-# - output spectra (-s) have ADCP-related corrections applied unless -c is set
+# May 18, 2015: - added -p)ulse <length>
+# Jun 11, 2015: - removed w_z code (yfname param requirement)
+# Jun 12, 2015: - adapted to &antsParam()
+# - BUG: %k.0 and %lambda.0 had been required erroneously (.1 are first ones to be used)
+# - made finescale limits optional
+# - renamed -b=>-l, -c=>-o, -s=>-i
+# Jun 14, 2015: - renamed -p=>-a
+# - added -p
+# - removed weird evals
+# Jun 15, 2015: - added plot & other mods
+# Jun 16, 2015: - define default outputs when STDOUT is a tty
+# - adapted to re-added pwrdens.0 in LADCP_wspec output
+# - modified plot label
+# Jun 17, 2015: - added eps.w to plot
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
$ANTS_TOOLS_AVAILABLE = (`which list` ne '');
require "$ANTSLIB/ants.pl";
require "$ANTSLIB/libLADCP.pl";
+require "$ANTSLIB/libGMT.pl";
+
+use FileHandle;
+use IPC::Open2;
#----------------------------------------------------------------------
# Usage
#----------------------------------------------------------------------
-my($c) = 0.0215;
+my($c) = 0.0215; # Thurnherr et al. (GRL 2015)
+
+&antsUsage('bc:de:g:i:l:mno:p:r:s:tuw:',0,
+ '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [apply cosine-t)aper]',
+ '[-d)own/-u)pcast-only] [exclude -b)ottom window]', # LADCP_wspec options
+ '[-s)urface <layer depth to exclude[150m]>',
+ '[-g)ap <max depth layer to fill with interpolation[40m]>]',
+ '[-w)indow <power-of-two input-records[32]>]',
+ '[shortwave -c)utoff <kz or lambda[100m]>]', # LADCP_VKE options
+ '[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]',
+ '[-l) override ADCP bin <length>] [-a) override pulse <length>]',
+ "[-e)ps-parameterization <scale[$c]>",
+ '[output -i)ndividual spectra <basename>]',
+ '[output -p)lot <ps-file>',
+ '[file]');
+
+#----------------------------------------------------------------------
+# Calculate VKE spectra with [LADCP_wspec] if input is a w_ocean file
+#----------------------------------------------------------------------
-&antsUsage('cne:r:s:',2,
- '[suppress spectral -c)orrection] [tilt-correction -r)ange <max[0m]>]',
- '[output -s)pectra <basename>]',
- "[-e)ps-parameterization <scale[$c]>",
- '<lambda.min> <lambda.max> [file...]');
+if (defined(fnrNoErr('dc_w'))) { # pre-process with LADCP_wspec when handed vertical-velocity input
+ my($opts);
+ $opts .= ' -d' if ($opt_d);
+ $opts .= ' -u' if ($opt_u);
+ $opts .= ' -b' if ($opt_b);
+ $opts .= ' -t' if ($opt_t);
+ $opts .= " -s $opt_s" if defined($opt_s);
+ $opts .= " -g $opt_g" if defined($opt_g);
+ $opts .= " -w $opt_w" if defined($opt_w);
+ $opts .= " -o $opt_o" if defined($opt_o);
+ open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # spawn sub-process
+ croak("LADCP_wspec $opts: $!\n");
+ open(STDIN,"<&FROMCLD") || croak("dup(FROMCLD): $!\n");
+ close(FROMCLD);
+ print(TOCLD $antsOldHeaders); undef($antsOldHeaders); # feed already gobbled header & first record to child
+ print(TOCLD $antsPeekBuffer); undef($antsPeekBuffer);
+ undef(%P); # shouldn't matter, because we'll get the same %PARAMs back
+ undef(@antsLayout); # shouldn't matter, because it will get overwritten
+ while (<>) { print(TOCLD $_); } # feed remaining data
+ close(TOCLD);
+} elsif (defined(fnrNoErr('pwrdens.1'))) {
+ croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n")
+ if ($opt_d || $opt_u || $opt_b || defined($opt_w) ||
+ defined($opt_s) || defined($opt_g));
+} else {
+ if ($ARGV[0]) {
+ croak("$ARGV[0]: no such file or directory\n");
+ } else {
+ croak("$0: empyt input\n");
+ }
+}
-$lmin = &antsFloatArg(); # wavelength limits
-$lmax = &antsFloatArg();
+#----------------------------------------------------------------------
+# Handle LADCP_VKE usage & read data
+#----------------------------------------------------------------------
+
+&antsAddParams('ADCP_bin_length',$opt_l) # override bin length
+ if defined($opt_l);
+&antsAddParams('ADCP_pulse_length',$opt_a) # override pulse length
+ if defined($opt_a);
&antsFloatOpt(\$opt_e,$c); # default parameterization
+if (defined($opt_c)) { # shortwave cutoff supplied
+ $lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c;
+ &antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.w is calculated below
+} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied
+ $lmin = 2*$PI/antsParam('shortwave_cutoff');
+} else { # use 100m default cutoff
+ $lmin = 100;
+ &antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.w is calculated below
+}
+$lmax = 9e99; # no longwave cutoff implemented yet
+
+&antsInstallBufFull(0); # load entire file
+&antsIn();
+my($Hbuf) = $antsOldHeaders; # save for later
+
+&antsRequireParam('profile_id');
+&antsRequireParam('lambda.1');
+&antsRequireParam('k.1');
+&antsRequireParam('resolution_bandwidth');
+&antsRequireParam('input_depth_resolution');
+&antsRequireParam('output_depth_resolution');
+&antsRequireParam('ADCP_bin_length');
+&antsAddParams('ADCP_pulse_length',antsParam('ADCP_bin_length'))
+ unless defined(antsParam('ADCP_pulse_length'));
+
$imin = 0; # find frequency bin limits
-for ($nfreq=1; defined($P{"binpgrams::lambda.$nfreq"}); $nfreq++) {
- $imin = $nfreq if ($imin==0 && $P{"binpgrams::lambda.$nfreq"}<=$lmax);
- $imax = $nfreq if ($P{"binpgrams::lambda.$nfreq"} >= $lmin);
+for ($nfreq=1; defined(antsParam("lambda.$nfreq")); $nfreq++) {
+ $imin = $nfreq if ($imin==0 && antsParam("lambda.$nfreq")<=$lmax);
+ $imax = $nfreq if (antsParam("lambda.$nfreq") >= $lmin);
}
croak("$0: <lambda.min=$lmin> < min(lambda)")
unless defined($imax);
@@ -58,26 +147,49 @@
$fs_fmin = $pg_fmin + $imin; # first power field in finescale range
$fs_fmax = $pg_fmin + $imax; # last power field in finescale range
+$widxf = fnr('widx');
$mindf = fnr('depth.min');
$maxdf = fnr('depth.max');
#----------------------------------------------------------------------
-# Calculate Finescale Quantities (depending on input quantity)
+# Redirect STDOUT & create plot if STDOUT is a tty
#----------------------------------------------------------------------
-sub integrate_fs_power($) # integrate fs spectrum (always done)
+if (-t STDOUT) {
+ my($id) = &antsRequireParam('profile_id');
+ $opt_p = sprintf('%03d_VKE.ps',$id)
+ unless defined($opt_p);
+ $outfile = sprintf('%03d.VKE',$id);
+ open(STDOUT,">$outfile") || die("$outfile: $!\n");
+}
+
+#----------------------------------------------------------------------
+# Library
+#----------------------------------------------------------------------
+
+sub integrate_fs_power($) # integrate fs spectrum
{
my($r) = @_;
$ants_[$r][$fspwrf] = 0;
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
- $ants_[$r][$fspwrf] += $ants_[0][$f];
+ $ants_[$r][$fspwrf] += $ants_[$r][$f];
}
- $ants_[$r][$fspwrf] *= $P{'binpgrams::resolution_bandwidth'};
+ $ants_[$r][$fspwrf] *= antsParam('resolution_bandwidth');
}
-sub fit_universal_w_spec($) # vertical velocity => p0
+sub normalize_spectral_power($) # scale pwrdens by p0
+{
+ my($r) = @_;
+
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
+ $ants_[$r][$f] /= $ants_[$r][$p0f];
+ }
+}
+
+
+sub fit_universal_w_spec($) # vertical velocity => p0
{
my($r) = @_;
my($nsamp) = $fs_fmax - $fs_fmin + 1;
@@ -86,27 +198,27 @@
# fit slope-2 line in log-log space (main estimator)
#---------------------------------------------------
- if ($nsamp >= 2) { # require min 2 wavenumber samples
+ if ($nsamp >= 2) { # require min 2 wavenumber samples
my($sumd,$sumx,$sumy) = (0,0,0);
- for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { # loop over wavenumbers
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { # loop over wavenumbers
my($i) = $f - $pg_fmin;
- $sumd += log10($ants_[$r][$f]) + 2*log10($P{"binpgrams::k.$i"});
- $sumx += log10($P{"binpgrams::k.$i"});
+ $sumd += log10($ants_[$r][$f]) + 2*log10(antsRequireParam("k.$i"));
+ $sumx += log10(antsParam("k.$i"));
$sumy += log10($ants_[$r][$f]);
}
my($p0) = $sumd/$nsamp;
$ants_[$r][$p0f] = 10**$p0;
- my($avgx) = $sumx/$nsamp; # avg for r calc
+ my($avgx) = $sumx/$nsamp; # avg for r calc
my($avgy) = $sumy/$nsamp;
- my($sumsq,$sxx,$syy,$sxy) = (0,0,0,0); # calc rms misfit & correlation coeff
+ my($sumsq,$sxx,$syy,$sxy) = (0,0,0,0); # calc rms misfit & correlation coeff
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
my($i) = $f - $pg_fmin;
- my($xt) = log10($P{"binpgrams::k.$i"}) - $avgx;
+ my($xt) = log10(antsParam("k.$i")) - $avgx;
my($yt) = log10($ants_[$r][$f]) - $avgy;
- $sumsq += ($p0 - 2*log10($P{"binpgrams::k.$i"}) - log10($ants_[$r][$f]))**2;
+ $sumsq += ($p0 - 2*log10(antsParam("k.$i")) - log10($ants_[$r][$f]))**2;
$sxx += $xt * $xt; $syy += $yt * $yt; $sxy += $xt * $yt;
}
$ants_[$r][$rmsf] = sqrt($sumsq/$nsamp);
@@ -124,14 +236,11 @@
$fspwrf = &antsNewField('fs_pwr'); # derived fields
-if ($P{binpgrams::yfname} eq 'dc_w' || # w-finestructure parameterization
- $P{binpgrams::yfname} eq 'uc_w') {
- $p0f = &antsNewField('p0');
- $rmsf = &antsNewField('p0fit.rms');
- $sampf = &antsNewField('p0fit.nsamp');
- $rf = &antsNewField('p0fit.r');
- $wepsf = &antsNewField('eps.w');
-}
+$p0f = &antsNewField('p0'); # VKE parameterization
+$rmsf = &antsNewField('p0fit.rms');
+$sampf = &antsNewField('p0fit.nsamp');
+$rf = &antsNewField('p0fit.r');
+$wepsf = &antsNewField('eps.w');
my(@outLayout) = @antsNewLayout; # save for later
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input
@@ -140,11 +249,24 @@
}
croak("$0: cannot find fields 'pwr.tot' or 'pwr.tot.nsamp' in input\n")
unless defined($totf) || defined($tnsf);
-$lsf = defined($tnsf) ? $tnsf : $totf;
+$lsf = defined($tnsf) ? $tnsf : $totf;
-&antsInstallBufFull(0); # load entire file
-&antsIn();
-my($Hbuf) = $antsOldHeaders; # save for later
+if ($opt_p) { # begin plot
+ $xmin = (&antsParam('output_depth_resolution')>=350) ? 0.012 : 0.018;
+ $xmax = 0.2; $ymin = 1; $ymax = 1e4;
+ $plotsize = 13;
+ GMT_begin($opt_p,"-JX${plotsize}l","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4'); # init plot
+ GMT_psxy('-G25/255/25 -L');
+ printf(GMT "%f %f\n",$xmin,$xmin**(-2)*sqrt(2));
+ printf(GMT "%f %f\n",$xmax,$xmax**(-2)*sqrt(2));
+ printf(GMT "%f %f\n",$xmax,$xmax**(-2)/sqrt(2));
+ printf(GMT "%f %f\n",$xmin,$xmin**(-2)/sqrt(2));
+}
+
+my(@sumPwr,@nPwr);
+my(@sumGoodPwr,@nGoodPwr);
+my($min_depth) = 9e99;
+my($max_depth) = -9e99;
for (my($r)=0; $r<@ants_; $r++) { # loop over all windows
@@ -152,17 +274,12 @@
# apply spectral correction
#--------------------------
- unless ($opt_c) {
+ unless ($opt_m) {
for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers
- if ($P{binpgrams::yfname} eq 'dc_w' || # vertical velocity
- $P{binpgrams::yfname} eq 'uc_w') {
- $ants_[$r][$i+$pg_fmin] *= T_w($P{"binpgrams::k.$i"},$P{ADCP_bin_length},$P{binpgrams::depth_resolution},$opt_r);
- } elsif ($P{binpgrams::yfname} eq 'dc_w_depth'|| # dw/dz
- $P{binpgrams::yfname} eq 'uc_w_depth') {
- $ants_[$r][$i+$pg_fmin] *= T_w_z($P{"binpgrams::k.$i"},$P{ADCP_bin_length},$P{binpgrams::depth_resolution},$opt_r);
- } else {
- croak("$0: cannot calculate finescale $P{'binpgrams::yfname'} spectra (implementation restriction)\n");
- }
+ $ants_[$r][$i+$pg_fmin] *=
+ T_w(antsParam("k.$i"),antsParam('ADCP_bin_length'),
+ antsParam('ADCP_pulse_length'),antsParam('input_depth_resolution'),
+ $opt_r);
}
}
@@ -170,19 +287,54 @@
# calculate fs quantities
#------------------------
- integrate_fs_power($r); # always integrate spectra
- fit_universal_w_spec($r)
- if (($P{binpgrams::yfname} eq 'dc_w') || ($P{binpgrams::yfname} eq 'uc_w'));
- $ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2;
+ integrate_fs_power($r); # calculate total finescale power
+ fit_universal_w_spec($r); # fit kz^-2 spectrum
+# normalize_spectral_power($r); # normalize by p0 (even w/o shortwave cutoff)
+
+ $ants_[$r][$wepsf] = $ants_[$r][$rmsf] <= 0.4 ? # do not calc eps unless good fit
+ ($ants_[$r][$p0f] / $opt_e)**2 : nan;
+
+ if (numberp($ants_[$r][$p0f])) { # update min/max depth
+ $min_depth = $ants_[$r][$mindf] if ($ants_[$r][$mindf] < $min_depth);
+ $max_depth = $ants_[$r][$maxdf] if ($ants_[$r][$maxdf] > $max_depth);
+ }
#---------------
# produce output
#---------------
- if (defined($opt_s)) { # output individual spectra
+ if (defined($opt_p)) { # plot current spectrum on -p & calc mean
+ my($R) = 0; # RGB map
+ my($G) = int(200*(1-$r/@ants_));
+ my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_));
+ GMT_psxy("-W8,$R/$G/$B");
+ for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { # avg & plot high-quality spectra only
+ next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]); # - omit zero wavenumber ($pg_fmin)
+ next if ($ants_[$r][$rmsf] > 0.4);
+ my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
+ printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]);
+ $nGoodPwr[$f]++;
+ $sumGoodPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f];
+ }
+ GMT_psxy("-W4,$R/$G/$B,."); # avg & plot all spectra with dots
+ for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) {
+ next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]);
+ my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
+ printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]);
+ $nPwr[$f]++;
+ $sumPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f];
+ }
+ }
+
+
+ if (defined($opt_i)) { # output current spectrum on -i
open(STDOUT_DUP,">&",STDOUT) || croak("$0: cannot dup STDOUT\n");
@antsNewLayout = ('k','lambda','pwrdens','finescale','pwrdens_fit');
- $ofname = sprintf("${opt_s}_%02d.wspec",$r+1);
+ if ($ants_[$r][$widxf] > 0) {
+ $ofname = sprintf("${opt_i}_%02d.wspec",$r+1);
+ } else {
+ $ofname = sprintf("${opt_i}_btm.wspec");
+ }
open(STDOUT,">$ofname") || croak("$ofname: $!\n");
undef($antsOldHeaders);
&antsActivateOut() if ($ANTS_TOOLS_AVAILABLE);
@@ -194,9 +346,9 @@
&antsAddParams($outLayout[$f],$ants_[$r][$f]);
}
- for (my($f)=$pg_fmin; $f<$pg_fmin+$nfreq; $f++) {
- my($k) = eval(sprintf('$P{"binpgrams::k.%d"}',$f-$pg_fmin));
- my($l) = eval(sprintf('$P{"binpgrams::lambda.%d"}',$f-$pg_fmin));
+ for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) {
+ my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
+ my($l) = antsParam(sprintf("lambda.%d",$f-$pg_fmin));
&antsOut($k,$l,$ants_[$r][$f],
($f>=$fs_fmin && $f<=$fs_fmax),
numberp($ants_[$r][$p0f]) ? $ants_[$r][$p0f] * $k**(-2) : nan);
@@ -209,6 +361,59 @@
}
}
+if (defined($opt_p)) { # complete plot
+ if (@nPwr) {
+ GMT_psxy('-W10,255/100/100,.'); # plot mean spectrum; dotted; entire range
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
+ my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
+ printf(GMT "$k %g\n",$sumPwr[$f]/$nPwr[$f]);
+ }
+ }
+
+ if (@nGoodPwr) {
+ GMT_psxy('-W4/255/100/100'); # plot mean spectrum; entire range
+ for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) {
+ my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
+ printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]);
+ }
+ GMT_psxy('-W16/255/100/100'); # plot mean fs spectrum (heavy)
+ for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {
+ my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin));
+ printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]);
+ }
+ }
+
+ GMT_psbasemap('-Bf3a2:"Vertical Wavenumber [rad/m]":/' . # annotate axes
+ 'f3a1:"VKE/p@-0@- [m@+3@+s@+-1/2@+/(rad m@+-1@+)]":WeS');
+ GMT_setJ("-JX-${plotsize}l");
+ GMT_setR(sprintf('-R%f/%f/%f/%f',2*$PI/$xmax,2*$PI/$xmin,$ymin,$ymax));
+ GMT_psbasemap('-B3:"Vertical Wavelength [m]:N"');
+
+ GMT_unitcoords(); # print profile number
+ GMT_pstext();
+ if (defined($outfile)) { printf(GMT "0.98 0.98 14 0 0 TR $outfile %s\n",antsParam('input_data')); }
+ else { printf(GMT "0.98 0.98 14 0 0 TR %03d %s\n",antsParam('profile_id'),antsParam('input_data')); }
+
+ GMT_set('ANNOT_FONT_SIZE_PRIMARY 10','ANNOT_OFFSET_PRIMARY 0.01c', # eps profile inset
+ 'LABEL_FONT_SIZE 10','LABEL_OFFSET 0');
+ $min_depth = 400 if ($min_depth > 1e4);
+ $max_depth = 1500 if ($max_depth < 0);
+ GMT_setR(sprintf("-R5e-12/1e-7/%d/%d",round($min_depth-250,500),round($max_depth+250,500)));
+ GMT_setJ(sprintf('-JX%fl/-%f',$plotsize/2.2,$plotsize/2.2));
+ GMT_psbasemap('-X0.3 -Y0.3 -Bg1a1f3p:"@~e@~@-VKE@- [W kg@+-1@+]":/g1000a500f100:"Depth [m]":wEsN');
+ GMT_psxy('-M');
+ for (my($r)=0; $r<@ants_; $r++) {
+ next unless numberp($ants_[$r][$wepsf]);
+ my($R) = 0; my($G) = int(200*(1-$r/@ants_));
+ my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_));
+ print(GMT "> -W12,$R/$G/$B\n");
+ print(GMT "$ants_[$r][$wepsf] $ants_[$r][$mindf]\n");
+ print(GMT "$ants_[$r][$wepsf] $ants_[$r][$maxdf]\n");
+ }
+
+ GMT_end(); # finish plot
+}
+
@antsNewLayout = @outLayout;
$antsOldHeaders = $Hbuf;
$antsHeadersPrinted = 0;
--- a/LADCP_w Mon Apr 20 14:04:04 2015 +0000
+++ b/LADCP_w Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P _ W
# doc: Fri Dec 17 18:11:13 2010
-# dlm: Mon Apr 20 14:01:10 2015
+# dlm: Sun Jul 26 17:11:01 2015
# (c) 2010 A.M. Thurnherr
-# uE-Info: 167 78 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 1465 90 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# TODO:
@@ -166,6 +166,20 @@
# Apr 20, 2015: - improved comments
# - improved diagnostic messages (warning on missing CTD temp)
# - added support for empirical backscatter correction
+# Apr 21, 2015: - improved screen log output
+# Apr 22, 2015: - BUG: $realLastGoodEns could be undefined, breaking plotting routines
+# Apr 24, 2015: - removed @ARGS
+# - added %profile_id
+# May 13, 2015: - loosened "insufficient valid velocities" test for short casts
+# May 15, 2015: - added $min_valid_vels
+# - BUG: LADCP_atbottom could be less than firstGoodEns
+# May 18, 2015: - added %LADCP_pulse_length, %dnXX
+# May 20, 2015: - added $PROF as newer alternative to $STN
+# - replaced "require ProcessingParams" by "do ProcessingParams"
+# Jun 15, 2015: - clarified warning message
+# Jul 26, 2015: - added %output_grid_dz %output_grid_minsamp for use by [LADCP_w_regrid]
+# - began work on support for [libGMT.pl]
+# - -v usage message had wrong default
# CTD REQUIREMENTS
# - elapsed elapsed seconds; see note below
@@ -251,11 +265,11 @@
# Usage
#------
-@ARGS = @ARGV; # save opts & arguments
+#@ARGS = @ARGV; # save opts & arguments
$antsParseHeader = 0;
&antsUsage('3:4a:b:c:e:g:h:i:k:m:n:o:p:qs:t:uv:w:x:',1,
- '[-v)erbosity <level[2]>]',
+ '[-v)erbosity <level[1]>]',
'[-q)uick (no single-ping denoising)]',
'[require -4)-beam solutions] [apply beamvel-m)ask <file> if it exists]',
'[valid LADCP -b)ins <bin[2],bin[*]>',
@@ -269,17 +283,17 @@
'[pressure-sensor -a)cceleration correction <residual/CTD_w_tt>]',
'[-o)utput bin <resolution[10m]>] [-k) require <min[20]> samples]',
'[e-x)ecute <perl-expr>]',
- '<station-id> [run-label]');
+ '<profile-id> [run-label]');
&antsUsageError() if ($opt_u && !defined($opt_i));
-&antsCardOpt(\$opt_s,0); # skip # initial ensembles
-$opt_p = '+' unless defined($opt_p); # separate uc/dc time lagging
+&antsCardOpt(\$opt_s,0); # skip # initial ensembles
+$opt_p = '+' unless defined($opt_p); # separate uc/dc time lagging
-&antsFloatOpt(\$opt_a,1); # pressure acceleration correction
+&antsFloatOpt(\$opt_a,1); # pressure acceleration correction
-$STN = &antsCardArg(); # station id
-if (@ARGV) { # run label
+$PROF = $STN = &antsCardArg(); # station id ($STN for compatibility)
+if (@ARGV) { # run label
$RUN = $ARGV[0];
shift;
} else {
@@ -291,7 +305,7 @@
#-----------------------------
require "$WCALC/defaults.pl"; # load default/option parameters
-require "$processing_param_file"; # load processing parameters
+do "$processing_param_file"; # load processing parameters
eval($opt_x) if defined($opt_x); # eval cmd-line expression to override anything
$RDI_Coords::minValidVels = 4 if ($opt_4); # decode cmd-line arguments
@@ -310,7 +324,7 @@
croak("$0: \$out_basename undefined\n") # all plotting routines use this
unless defined($out_basename);
&antsAddParams('out_basename',$out_basename);
-&antsAddParams('run_label',$RUN);
+&antsAddParams('profile_id',$PROF,'run_label',$RUN);
#----------------------------------------------------------------------
# Screen Logging
@@ -333,8 +347,7 @@
sub info(@)
{
print(LOGF "\t"),printf(LOGF @_) if defined($out_log);
- print(STDERR ($opt_v > 1) ? "\t" : "$LADCP_file: ");
- printf(STDERR @_) if ($opt_v > 0);
+ printf(STDERR @_) if ($opt_v > 1);
}
sub warning(@)
@@ -342,15 +355,20 @@
my($lvl,@msg) = @_;
if (defined($out_log)) {
- print(LOGF "\nWARNING (level $lvl): ");
+ print(LOGF "\nWARNING (L$lvl): ");
printf(LOGF @msg);
print(LOGF "\n");
}
return if ($opt_v == 0);
- print(STDERR "\n") if ($opt_v > 1);
- print(STDERR "WARNING (level $lvl): ");
- printf(STDERR @msg);
- print(STDERR "\n") if ($opt_v > 1);
+
+ if ($opt_v == 1) {
+ print(STDERR "$LADCP_file: WARNING (L$lvl): ");
+ printf(STDERR @msg);
+ } else {
+ print(STDERR "\n-------------\nWARNING (L$lvl): ");
+ printf(STDERR @msg);
+ print(STDERR "-------------\n\n")
+ }
}
sub error($)
@@ -398,6 +416,7 @@
if ($LADCP{BLANKING_DISTANCE}<$LADCP{BIN_LENGTH} && $refLr_firstBin==1);
&antsAddParams('ADCP_bin_length',$LADCP{BIN_LENGTH},
+ 'ADCP_pulse_length',$LADCP{TRANSMITTED_PULSE_LENGTH},
'ADCP_frequency',$LADCP{BEAM_FREQUENCY},
'ADCP_blanking_distance',$LADCP{BLANKING_DISTANCE});
@@ -534,7 +553,7 @@
progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");
}
-error("$LADCP_file: insufficient valid velocities\n") unless ($nvw > 1000);
+error("$LADCP_file: insufficient valid velocities\n") unless ($nvw >= $min_valid_vels);
#----------------------------------------------
# STEP: Edit earth-coordinate -velocity data
@@ -567,14 +586,17 @@
($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =
calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$opt_g);
-error("$LADCP_file: no good ensembles\n")
- unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0);
+error("$LADCP_file: insufficient valid data\n")
+ unless defined($firstGoodEns) && ($lastGoodEns>$firstGoodEns) && ($LADCP_atbottom>=$firstGoodEns);
my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} -
$LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED};
error("$0: implausibly short cast ($cast_duration seconds)\n")
unless ($cast_duration > 600);
+my($year) = $LADCP{ENSEMBLE}[$firstGoodEns]->{YEAR} % 100;
+&antsAddParams("dn$year",$LADCP{ENSEMBLE}[$firstGoodEns]->{DAYNO});
+
$LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1);
progress("\tStart of cast : %s (#%5d)\n",
@@ -658,7 +680,7 @@
($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w');
$CTD_temp = &fnrNoErr('temp');
-warning(2,"no CTD temperature --- using ADCP temperature instead => Sv degraded!\n",$s)
+warning(0,"no CTD temperature --- using ADCP temperature instead => Sv degraded!\n",$s)
unless defined($CTD_temp);
$CTD_maxdepth = -1;
@@ -792,9 +814,10 @@
while (@splits > 1) {
push(@CTD_tl_fromEns,$splits[0]);
push(@CTD_tl_toEns,$splits[1]);
+ debugmsg("lag($splits[0],$splits[1]) = ");
my($lag) = calc_lag($number_of_timelag_windows[1],$length_of_timelag_windows[1],
1,$splits[0],$splits[1]);
- debugmsg("lag($splits[0],$splits[1]) = $lag\n");
+ debugmsg("$lag\n");
if (defined($lag)) {
progress("\tcast-piece: elapsed(CTD) = elapsed(LADCP) + %.2fs\n",$lag);
push(@CTD_time_lag,$lag);
@@ -896,12 +919,10 @@
progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n",
100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI));
if (sqrt($sumWsqI/$nWsqI) > 0.05) {
- warning(1,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
+ warning(1,"%.2f cm/s (large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
} elsif (sqrt($sumWsqI/$nWsqI) > 0.1) {
- warning(2,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
+ warning(2,"%.2f cm/s (very large) reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI));
}
-# error("$0: rms reference-layer w_ocean is too large\n")
-# unless (sqrt($sumWsqI/$nWsqI) < 0.07);
} elsif ($nWsq > 0) {
&antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan);
progress("\t%.2f cm/s rms reference-layer w_ocean\n",100*sqrt($sumWsq/$nWsq));
@@ -1066,9 +1087,11 @@
#---------------------------------------------------------------------------
progress("Creating binned profiles at ${opt_o}m resolution...\n");
+&antsAddParams('output_grid_dz',$opt_o,'output_grid_minsamp',$opt_k); # will be used by LADCP_w_regrid
my($min_depth) = 9e99;
my($max_depth) = 0;
+my($realLastGoodEns);
progress("\tdowncast...\n");
for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast
@@ -1076,6 +1099,7 @@
$firstGoodEns++ if ($ens == $firstGoodEns); # start has been edited away
next;
}
+ $realLastGoodEns = $ens;
my(@bindepth) = calc_binDepths($ens);
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);
@@ -1125,7 +1149,6 @@
progress("\tupcast...\n"); # upcast
-my($realLastGoodEns);
for ($ens=$LADCP_atbottom; $ens<=$lastGoodEns; $ens++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH});
$realLastGoodEns = $ens;
@@ -1438,7 +1461,15 @@
foreach my $of (@out_profile) {
progress("<$of> ");
- $of = ">$of" unless ($of =~ /^$|^\s*\|/);
+
+ my($plot,$out) = ($of =~ /^([^\(]+)\(([^\)]+)\)$/); # plot_sub(out_file)
+ if (defined($out)) {
+ require "$WCALC/${plot}.pl";
+ &{$plot}($out);
+ next;
+ }
+
+ $of = ">$of" unless ($of =~ /^$|^\s*\|/); # pipe or file output
open(STDOUT,$of) || error("$of: $!\n");
undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/LADCP_w_CTD Sun Jul 26 20:04:48 2015 +0000
@@ -0,0 +1,520 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P _ W _ C T D
+# doc: Mon Nov 3 17:34:19 2014
+# dlm: Mon Jul 20 09:28:16 2015
+# (c) 2014 A.M. Thurnherr
+# uE-Info: 38 25 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
+$antsMinLibVersion = 6.0;
+
+# HISTORY:
+# Nov 3, 2014: - created
+# Nov 4, 2014: - improved
+# Nov 6, 2014: - BUG: sound speed was not calculated correctly
+# - added -a
+# - added conductivity & temperature editing
+# Nov 7, 2014: - loosened outlier editing
+# - added no-valid-data error message
+# - modified binning criterion to allow any sampling
+# frequency (not just divisors of 24)
+# Apr 17, 2015: - added in-situ temperature to output
+# Apr 23, 2015: - had to losen unrealistic-soundspeed criterion
+# because of DUK2#45
+# May 13, 2015: - BUG: did not work with casts beginning at depth
+# - added errors when editing edits too much
+# May 14, 2015: - BUG: depth was wrong for partial-depth casts
+# - BUG: $badval had not been handled correctly
+# - BUG: editing errors were too tight
+# Jun 16, 2015: - added STDOUT redirection (default output file) on tty stdout
+# - added -p)lot option
+# Jun 17, 2015: - added cond to std output
+# - added conductivity spike filter (different from pressure)
+# Jun 18, 2015: - renamed from SBE_w
+# Jun 19, 2015: - added press & salin to make output suitable for u/v processing
+# - removed a couple of data assertions
+# Jul 20, 2015: - allow for non-numeric -i
+
+($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
+
+require "$ANTS/ants.pl";
+require "$ANTS/fft.pl";
+require "$ANTS/libGMT.pl";
+require "$ANTS/libSBE.pl";
+require "$ANTS/libEOS83.pl";
+
+$antsParseHeader = 0; # usage
+&antsUsage('ai:l:rp:s:v:',1,
+ '[-v)erbosity <level[0]>]',
+ '[use -a)lternate sensor pair]',
+ '[-r)etain all data (no editing)]',
+ '[-s)ampling <rate[6Hz]>]',
+ '[-l)owpass w <cutoff[2s]>]',
+ '[profile -i)d <id>]',
+ '[-p)lot_basenames <[%03d_wpkg.ps],[%03d_sspd.ps]>]',
+ '<SBE CNV file>');
+
+&antsFloatOpt(\$opt_l,2); # defaults
+&antsCardOpt(\$opt_s,6);
+&antsCardOpt(\$opt_v,0);
+
+$CNVfile = $ARGV[0]; # open CNV file
+open(F,&antsFileArg());
+&antsActivateOut(); # activate ANTS file
+
+#----------------------------------------------------------------------
+# Read Data
+#----------------------------------------------------------------------
+
+print(STDERR "Reading $CNVfile...") if ($opt_v);
+
+($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header
+ SBE_parseHeader(F,0,0); # SBE field names, no time check
+
+croak("$CNVfile: unexpected sampling interval $sampint\n")
+ unless (abs($sampint-1/24) < 1e-5);
+croak("$CNVfile: unknown latitude\n")
+ unless numberp($lat);
+
+$pressF = fnr('prDM');
+
+if ($opt_a) { # temp/cond alternate sensor pair
+ $tempF = fnr('t190C');
+ $condF = fnrNoErr('c1S/m');
+ if (defined($condF)) {
+ $condHistRes = 10; # 0.1 bins
+ } else {
+ $condF = fnr('c1mS/cm');
+ $condHistRes = 1; # 1.0 bins
+ }
+} else { # primary sensor pair
+ $tempF = fnr('t090C');
+ $condF = fnrNoErr('c0S/m');
+ if (defined($condF)) {
+ $condHistRes = 10;
+ } else {
+ $condF = fnr('c0mS/cm');
+ $condHistRes = 1;
+ }
+}
+
+&antsInstallBufFull(0); # read entire CNV file
+&SBEin(F,$ftype,$nfields,$nscans,$badval);
+
+printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1);
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Redirect STDOUT to %.6Hz & create %_wpkg.ps,%_sspd.ps if STDOUT is a tty
+#----------------------------------------------------------------------
+
+$id = defined($opt_i) ? $opt_i : &antsParam('station');
+croak("$CNVfile: no station information in header => -i required\n")
+ unless defined($id);
+
+if (-t STDOUT) {
+ if (numberp($id)) {
+ $opt_p = '%03d_wpkg.ps,%03d_sspd.ps'
+ unless defined($opt_p);
+ $outfile = sprintf('%03d.%dHz',$id,$opt_s);
+ } else {
+ $opt_p = '%s_wpkg.ps,%s_sspd.ps'
+ unless defined($opt_p);
+ $outfile = sprintf('%s.%dHz',$id,$opt_s);
+ }
+ open(STDOUT,">$outfile") || die("$outfile: $!\n");
+}
+
+#----------------------------------------------------------------------
+# Edit Data
+# - pressure outliers & spikes
+# - conductivity outliers & spikes
+#----------------------------------------------------------------------
+
+unless ($opt_r) {
+ print(STDERR "Editing Data...") if ($opt_v);
+
+ #----------------------------------------
+ # trim initial records with nan pressures
+ #----------------------------------------
+ my($trimmed) = 0; # trim leading nan pressures
+ shift(@ants_),$trimmed++
+ until numberp($ants_[0][$pressF]) && numberp($ants_[0][$condF]);
+ printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1);
+ my($lvp) = $ants_[0][$pressF];
+ my($lvc) = $ants_[0][$condF];
+
+ #------------------------------------------------
+ # edit pressure outliers outside contiguous range
+ #------------------------------------------------
+ my($outliers) = 0; my($min) = 9e99; my($max);
+ for (my($s)=0; $s<$nscans; $s++) {
+ next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
+ $phist[$ants_[$s][$pressF]+1000]++;
+ $min = $ants_[$s][$pressF] if ($ants_[$s][$pressF] < $min);
+ }
+ for ($max=$min+25; $phist[$max+1000]||$phist[$max+1001]; $max++) {} # outliers after 2 gaps
+ for ($min=$max; $phist[$min+1000]||$phist[$min+ 999]; $min--) {}
+ for (my($s)=0; $s<$nscans; $s++) {
+ if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
+ if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
+ }
+ &antsAddParams("pressure_outliers",sprintf("%d",$outliers));
+ printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
+ croak("$CNVfile: pressure editing removed too many 'outliers'\n")
+ unless ($outliers < 100);
+
+ #----------------------------------------------------
+ # edit conductivity outliers outside contiguous range
+ #----------------------------------------------------
+ $outliers = 0;
+ my($modeSamp)=0;
+ undef(@phist);
+ for (my($s)=0; $s<$nscans; $s++) {
+ next unless ($ants_[$s][$condF] > 0);
+ my($b) = $ants_[$s][$condF]*$condHistRes; # 1/10 S/m histogram resolution (1 mS/cm)
+ $phist[$b]++;
+ next unless ($phist[$b] > $modeSamp);
+ $modeSamp = $phist[$b]; $modeBin = $b;
+ }
+ for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {} # outliers after 2 gaps
+ for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}
+ $max /= $condHistRes; $min /= $condHistRes;
+ for (my($s)=0; $s<$nscans; $s++) {
+ if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }
+ if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }
+ }
+ &antsAddParams("conductivity_outliers",sprintf("%d",$outliers));
+ printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
+ croak("$CNVfile: conductivity editing removed too many 'outliers'\n")
+ unless ($outliers/$nscans < 0.4);
+
+ #----------------------------------------------------
+ # edit temperature outliers outside contiguous range
+ #----------------------------------------------------
+ $outliers = 0;
+ my($modeSamp)=0;
+ undef(@phist);
+ for (my($s)=0; $s<$nscans; $s++) {
+ next unless ($ants_[$s][$tempF] > 0);
+ my($b) = $ants_[$s][$tempF]*10; # 10th of a degree histogram resolution
+ $phist[$b]++;
+ next unless ($phist[$b] > $modeSamp);
+ $modeSamp = $phist[$b]; $modeBin = $b;
+ }
+ for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {} # outliers after 2 gaps
+ for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}
+ $max /= 10; $min /= 10;
+ for (my($s)=0; $s<$nscans; $s++) {
+ if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }
+ if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }
+ }
+ &antsAddParams("temperature_outliers",sprintf("%d",$outliers));
+ printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
+ croak("$CNVfile: temperature editing removed too many 'outliers'\n")
+ unless ($outliers/$nscans < 0.4);
+
+ #----------------------------------------
+ # edit pressure spikes based on gradients
+ #----------------------------------------
+
+ for (my($s)=1; $s<$nscans; $s++) { # calculate pressure gradients (across gaps)
+ if (numberp($ants_[$s][$pressF])) {
+ $dp[$s-1] = $ants_[$s][$pressF] - $lvp;
+ $lvp = $ants_[$s][$pressF];
+ } else {
+ $dp[$s-1] = nan;
+ }
+ }
+
+ my($ns1,$ns2) = (0,0);
+ for (my($s)=0; $s<$nscans-2; $s++) { # consecutive large pressure gradients of opposite sign
+ if (($dp[$s]*$dp[$s+1] < 0) && # tests return false if either of the dps is not defined
+ (abs($dp[$s]) > 10) &&
+ (abs($dp[$s+1]) > 10)) {
+ $ants_[$s+1][$pressF] = nan;
+ $dp[$s] = $dp[$s+1] = undef;
+ $ns1++;
+ }
+ }
+ for (my($s)=0; $s<$nscans-3; $s++) { # 3 consecutive large pressure gradients of opposite sign
+ if (($dp[$s]>2 && $dp[$s+1]<-4 && $dp[$s+2]>2) ||
+ ($dp[$s]<-2 && $dp[$s+1]>4 && $dp[$s+2]<-2)) {
+ $ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan;
+ $dp[$s] = $dp[$s+1] = $dp[$s+2] = undef;
+ $ns2+=2;
+ }
+ }
+ &antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));
+ printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1);
+
+ #--------------------------------------------------
+ # edit conductivity spikes based on large gradients
+ # - $lvc = $ants_[0] is guaranteed numeric here
+ #--------------------------------------------------
+
+ my($nv) = my($ns) = my($last_dcond) = my($lvcs) = 0;
+ my($run_start) = my($run_dcond) = 0;
+ for (my($s)=1; $s<$nscans; $s++) { # calculate conductivity gradients (across gaps)
+ next unless numberp($ants_[$s][$condF]);
+
+ my($dcond) = $ants_[$s][$condF] - $lvc; # integrate gradient across runs
+# print(STDERR "ldc/dc/rdc: $last_dcond/$dcond/$run_dcond\n");
+ if ($last_dcond*$dcond >= 0) { # run is continuing (gradient does not change sign)
+ $run_dcond += $dcond;
+ } else { # run has ended
+# print(STDERR "new run at $lvcs (run_dcond = $last_dcond)\n");
+ $run_start = $lvcs;
+ $run_dcond = $dcond;
+ }
+
+ if (abs($run_dcond) <= 0.02) { # small integrated gradient => okay
+ $lvc = $ants_[$s][$condF]; # update stored previous values
+ $lvcs = $s;
+ $last_dcond = $dcond;
+ next; # process next scan
+ }
+
+# print(STDERR "large gradient ($ants_[$run_start][$condF]-$ants_[$s][$condF], $run_dcond) run $run_start-$s\n");
+ my($i);
+ my($max_spike_length) = 24;
+ for ($i=$s; $i<=$run_start+$max_spike_length && $i<$nscans-1; $i++) { # large gradient => check whether values return within 10s
+ next unless (numberp($ants_[$i][$condF]) &&
+ numberp($ants_[$i+1][$condF]));
+ last if ((abs($ants_[$i][$condF]-$ants_[$run_start][$condF]) <= 0.005) && # 2 vals to avoid large gradients
+ (abs($ants_[$i+1][$condF]-$ants_[$run_start][$condF]) <= 0.005));
+ }
+
+ if ($i>$run_start+$max_spike_length || $i==$nscans-1) { # values don't return => leave data alone
+# print(STDERR "values don't return => new run at $s\n");
+ $run_start = $s; $run_dcond = 0;
+ $lvc = $ants_[$run_start][$condF]; # start new run
+ $lvcs = $run_start;
+ $last_dcond = 0;
+ next; # process next scan
+ }
+
+ $ns++;
+ for (my($j)=$run_start+1; $j<$i; $j++) { # values return => remove bad spike
+ $ants_[$j][$condF] = nan;
+ $nv++;
+ }
+# print(STDERR "values return at $i ($ants_[$i][$condF]) => deleted $run_start+1-$i; new run\n");
+ $run_start = $i; $run_dcond = 0;
+ $lvc = $ants_[$run_start][$condF];
+ $lvcs = $run_start;
+ $last_dcond = 0;
+ }
+ &antsAddParams("conductivity_spikes_removed",$ns);
+ printf(STDERR "\n\t%d conductivity values removed from %d spikes",$nv,$ns) if ($opt_v > 1);
+
+ printf(STDERR "\n") if ($opt_v);
+
+} # if $opt_r
+
+#----------------------------------------------------------------------
+# Correcting for pressure bias
+#----------------------------------------------------------------------
+
+print(STDERR "Correcting for pressure bias...") if ($opt_v);
+
+my($minP) = 9e99;
+for (my($s)=0; $s<$nscans; $s++) {
+ $minP = $ants_[$s][$pressF]
+ if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);
+}
+croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
+ unless ($minP < 9e99);
+
+if ($minP < 25) {
+ &antsAddParams('pressure_bias',$minP);
+ printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1);
+ for (my($s)=0; $s<$nscans; $s++) {
+ $ants_[$s][$pressF] -= $minP
+ if numberp($ants_[$s][$pressF]);
+ }
+} else {
+ printf(STDERR "\n\tpartial-depth cast below %.1f dbar (no correction applied)",$minP) if ($opt_v > 1);
+}
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Binning data
+#----------------------------------------------------------------------
+
+my($sps) = round(1 / $sampint / $opt_s);
+print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);
+&antsAddParams('sampling_frequency',1/$opt_s);
+&antsAddParams('sampling_rate',$opt_s);
+
+my(@press,@temp,@cond);
+my($sp,$np,$st,$nt,$sc,$nc);
+
+$sp = $st = $sc = $np = $nt = $nc = 0;
+for (my($rec)=1,my($s)=0; $s<$nscans; $s++) {
+ if ($s*$sampint > $rec/$opt_s) {
+ $rec++;
+ push(@press,$np>0?$sp/$np:nan);
+ push(@temp, $nt>0?$st/$nt:nan);
+ push(@cond, $nc>0?$sc/$nc:nan);
+ $sp = $st = $sc = $np = $nt = $nc = 0;
+ }
+ $sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]);
+ $st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]);
+ $sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]);
+}
+
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Calculating derived quantities
+#----------------------------------------------------------------------
+
+print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v);
+
+my($maxP) = -9e99; my($atBtm);
+my($min_sspd) = 9e99; my($max_sspd) = -9e99;
+for (my($r)=0; $r<@press; $r++) {
+ $maxP=$press[$r],$atBtm=$r if ($press[$r] > $maxP);
+ $elapsed[$r] = $r/$opt_s;
+ $depth[$r] = &depth($press[$r],$lat);
+ $salin[$r] = &salin($cond[$r],$temp[$r],$press[$r]);
+ $sspd[$r] = &sVel($salin[$r],$temp[$r],$press[$r]);
+ $min_sspd = $sspd[$r] if ($sspd[$r] < $min_sspd);
+ $max_sspd = $sspd[$r] if ($sspd[$r] > $max_sspd);
+}
+
+$w[0] = nan;
+for (my($r)=1; $r<@depth-1; $r++) {
+ $w[$r] = numbersp($depth[$r-1],$depth[$r+1])
+ ? ($depth[$r+1] - $depth[$r-1]) * $opt_s
+ : nan;
+}
+push(@w,nan);
+
+printf(STDERR "\n") if ($opt_v);
+
+#----------------------------------------------------------------------
+# Low-pass filter velocity data
+# - interpolate missing vertical velocities first
+#----------------------------------------------------------------------
+
+if ($opt_l > 0) {
+ print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);
+ &antsAddParams('w_lowpass_cutoff',$opt_l);
+
+ my($trimmed) = 0;
+ shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++
+ until numberp($w[0]);
+ my($interpolated) = 0;
+ for ($r=1; $r<@w; $r++) {
+ next if numberp($w[$r]);
+ my($lv) = $r-1;
+ for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {}
+ if ($nv < @depth) {
+ while ($r < $nv) {
+ $w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]);
+ $interpolated++;
+ $r++;
+ }
+
+ } else {
+ $trimmed += @w-$r;
+ splice(@w,$r); splice(@depth,$r);
+ splice(@elapsed,$r); splice(@sspd,$r);
+ }
+ }
+ &antsAddParams('w_interpolated',$interpolated);
+ printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1);
+
+
+#--------------------
+# Zero Pad Data
+#--------------------
+
+ for ($pot=1; $pot<@w; $pot<<=1) {} # determine power of two
+
+ for ($r=0; $r<@w; $r++) { # copy data
+ $fftbuf[2*$r] = $w[$r];
+ $fftbuf[2*$r+1] = 0;
+ }
+ printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1);
+ while ($r < $pot) { # pad with zeroes
+ $fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;
+ $r++;
+ }
+
+#--------------------
+# Low-Pass Filter
+#--------------------
+
+ @fco = &FOUR1(-1,@fftbuf); # forward FFT
+ $n = @fco/2;
+ for (my($ip)=2; $ip<=$n; $ip+=2) { # +ve freq fco
+ my($in) = 2*$n-$ip; # -ve freq fco
+ my($f) = $ip/2/$n*$opt_s; # frequency
+ $fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0
+ if ($f > 1/$opt_l); # low-pass filter
+ }
+ @w_lp = &FOUR1(1,@fco); # inverse FFT
+
+ printf(STDERR "\n") if ($opt_v);
+} else {
+ @w_lp = @w;
+}
+
+#----------------------------------------------------------------------
+
+if (defined($opt_p)) {
+ print(STDERR "Plotting data...\n") if ($opt_v);
+ my(@pfmt) = split(',',$opt_p);
+ croak("$0: cannot decode -p $opt_p\n")
+ unless (@pfmt == 2);
+
+ my($xmin) = $elapsed[0]/60;
+ my($xmax) = $elapsed[$#elapsed]/60;
+ my($ymin) = -3; my($ymax) = 3;
+ my($plotsize) = '13c';
+ GMT_begin(sprintf($pfmt[0],$id),"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4');
+ GMT_psxy('-W2,255/127/80');
+ for ($r=0; $r<@w; $r++) {
+ printf(GMT "%f %f\n",$elapsed[$r]/60,$w_lp[2*$r]/@w_lp);
+ GMT_psxy('-W2,46/139/87') if ($r == $atBtm);
+ }
+ GMT_psbasemap('-Bg60a30f5:"Elapsed Time [min]":/g1a1f0.1:"Vertical Package Velocity [ms@+-1@+]":WeSn');
+ GMT_unitcoords();
+ GMT_pstext();
+ if (defined($outfile)) { printf(GMT "0.98 0.98 14 0 0 TR $outfile\n",$id); }
+ else { printf(GMT "0.98 0.98 14 0 0 TR %03d\n",$id); }
+ GMT_end();
+
+ my($xmin) = round($min_sspd-3,5);
+ my($xmax) = round($max_sspd+3,5);
+ my($ymin) = 0; my($ymax) = round($depth[$atBtm]+70,100);
+ my($plotsize) = '13c';
+ GMT_begin(sprintf($pfmt[1],$id),"-JX${plotsize}/-${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4');
+ GMT_psbasemap('-Bg10a10f1:"Speed of Sound [m/s]":/g1000a500f100:"Depth [m]":WeSn');
+ GMT_psxy('-W8,255/127/80');
+ for ($r=0; $r<@w; $r++) {
+ printf(GMT "%f %f\n",$sspd[$r],$depth[$r]);
+ GMT_psxy('-W6,46/139/87') if ($r == $atBtm);
+ }
+ GMT_unitcoords();
+ GMT_pstext();
+ if (defined($outfile)) { printf(GMT "0.98 0.98 14 0 0 TR $outfile\n",$id); }
+ else { printf(GMT "0.98 0.98 14 0 0 TR %03d\n",$id); }
+ GMT_end();
+}
+
+#----------------------------------------------------------------------
+
+print(STDERR "Writing output...\n") if ($opt_v);
+
+@antsNewLayout = ('elapsed','press','temp','cond','depth','salin','sspd','w.raw','w');
+for ($r=0; $r<@w; $r++) {
+ &antsOut($elapsed[$r],$press[$r],$temp[$r],$cond[$r],$depth[$r],$salin[$r],$sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp);
+}
+
+exit(0); # don't flush @ants_
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/LADCP_w_regrid Sun Jul 26 20:04:48 2015 +0000
@@ -0,0 +1,314 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P _ W _ R E G R I D
+# doc: Fri Apr 24 17:15:59 2015
+# dlm: Sun Jul 26 10:26:54 2015
+# (c) 2015 A.M. Thurnherr
+# uE-Info: 234 43 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
+$antsMinLibVersion = 6.1;
+
+# HISTORY:
+# Apr 24, 2015: - created
+# Apr 25, 2015: - maded gridding work
+# Apr 26, 2015: - made editing work
+# Apr 27, 2015: - added -p
+# May 5, 2015: - modified Editfile syntax to use parens()
+# May 7, 2015: - allow leading whitespace before Editfile labels
+# May 17, 2015: - removed warning about missing ./Editfile
+# May 18, 2015: - added important %PARAMs for dual-head data
+# - updated to libV6.1
+# May 19, 2015: - added hab to output
+# - allow setting %PARAMS in Editfile
+# May 20, 2015: - Editfile => EditParams
+# Jun 18, 2015: - added -i
+# - implemented default output on -t stdout
+# - changed to libGMT.pl
+# - removed dead -d option
+# Jul 26, 2015: - adapted for %output_grid_dz
+
+($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
+
+die("$0: ANTSlib required but not found (bad \$PATH?)\n")
+ unless ($ANTS ne '');
+
+require "$ANTS/ants.pl";
+require "$ANTS/libstats.pl";
+require "$ANTS/libGMT.pl";
+
+&antsUsage('i:k:o:p:',1,
+ '[profile -i)d <id>]',
+ '[-o)utput bin <resolution[10m]>] [-k) require <min> samples]',
+ '[-p)lot <[%03d_wprof.eps]>]',
+ '<.samp file> [.samp file]');
+
+$dual_head = (@ARGV==1);
+
+$opt_o_override = defined($opt_o);
+$opt_o = $P{output_grid_dz};
+&antsCardOpt(\$opt_o,10);
+
+$opt_k_override = defined($opt_k);
+$opt_k = 2 * $P{ouput_grid_minsamp}
+ if defined($P{ouput_grid_minsamp});
+&antsCardOpt(\$opt_k,$dual_head?40:20);
+
+
+#----------------------------------------------------------------------
+# Redirect STDOUT to %.wprof & create %_wprof.ps if STDOUT is a tty
+#----------------------------------------------------------------------
+
+$id = defined($opt_i) ? $opt_i : &antsParam('profile_id');
+croak("$0: no profile_id in first file => -i required\n")
+ unless defined($id);
+
+if (-t STDOUT) {
+ $opt_p = '%03d_wprof.ps' unless defined($opt_p);
+ $outfile = sprintf('%03d.wprof',$id);
+ open(STDOUT,">$outfile") || die("$outfile: $!\n");
+}
+
+#----------------------------------------------------------------------
+# Library
+#----------------------------------------------------------------------
+
+my(@dscMin,@dscMax,@dscDUc);
+
+sub discard($$$) # define bad-data filter
+{
+ my($fnm,$min,$max) = @_;
+ my($fnr) = fnr($fnm);
+ $dscMin[$fnr] = numberp($min) ? $min : -9e99;
+ $dscMax[$fnr] = numberp($max) ? $max : 9e99;
+ $dscDUc[$fnr] = 2;
+}
+
+sub discard_dc($$$) # define bad-data filter
+{
+ my($fnm,$min,$max) = @_;
+ my($fnr) = fnr($fnm);
+ $dscMin[$fnr] = numberp($min) ? $min : -9e99;
+ $dscMax[$fnr] = numberp($max) ? $max : 9e99;
+ $dscDUc[$fnr] = 1;
+}
+
+sub discard_uc($$$) # define bad-data filter
+{
+ my($fnm,$min,$max) = @_;
+ my($fnr) = fnr($fnm);
+ $dscMin[$fnr] = numberp($min) ? $min : -9e99;
+ $dscMax[$fnr] = numberp($max) ? $max : 9e99;
+ $dscDUc[$fnr] = 0;
+}
+
+sub isBad()
+{
+ for (my($f)=0; $f<$antsBufNFields; $f++) {
+ next unless defined($dscDUc[$f]);
+ next unless ($dscDUc[$f]==2 || $dscDUc[$f]==$ants_[0][$dcF]);
+ return 1 unless ($ants_[0][$f] < $dscMin[$f]) ||
+ ($ants_[0][$f] > $dscMax[$f]);
+ }
+ return 0;
+}
+
+#----------------------------------------------------------------------
+# Edit Data
+#----------------------------------------------------------------------
+
+$wF = &fnr('w');
+$eF = &fnr('elapsed');
+$dF = &fnr('depth');
+$dcF = &fnr('downcast');
+
+$first_label = &antsRequireParam('run_label');
+$bin_length = &antsRequireParam('ADCP_bin_length');
+$pulse_length = &antsRequireParam('ADCP_pulse_length');
+$blanking_dist = &antsRequireParam('ADCP_blanking_distance');
+($dayNoP,$dn) = &antsFindParam('dn\d\d');
+croak("$0: cannot determine day number\n")
+ unless defined($dayNoP);
+
+my($R,$R2);
+if ($opt_p) { # begin plot
+ $xmin = -0.1; $x2min = -200;
+ $xmax = 0.35; $x2max = 200;
+ $ymin = 0;
+ $ymax = antsParam('water_depth');
+ $ymax = antsRequireParam('max_depth') unless defined($ymax);
+ $plotsize = 13;
+ $R = "-R$xmin/$xmax/$ymin/$ymax";
+ $R2 = "-R$x2min/$x2max/$ymin/$ymax";
+ GMT_begin(sprintf($opt_p,$id),"-JX$plotsize/-$plotsize",$R,'-X6 -Y4');
+ GMT_psxy('-W1');
+ print(GMT "0 $ymin\n0 $ymax");
+ GMT_psxy('-L -G200');
+ print(GMT "0.07 $ymin\n0.07 $ymax\n0.18 $ymax\n0.18 $ymin\n");
+ GMT_setR($R2);
+ GMT_psxy('-M -W1');
+ print(GMT ">\n50 $ymin\n50 $ymax\n");
+ print(GMT ">\n100 $ymin\n100 $ymax\n");
+ print(GMT ">\n150 $ymin\n150 $ymax\n");
+ GMT_setR($R);
+}
+
+$min_depth = 9e99; # sentinels
+$max_depth = -9e99;
+
+$curF = ''; # current input file (sentinel)
+$filt = 0;
+for (my($curF),$r=0; &antsIn(); $r++) {
+ if ($P{PATHNAME} ne $curF) { # new file
+ $curF = $P{PATHNAME};
+
+ &antsInfo("WARNING: inconsistent %output_grid_dz in input files")
+ if (defined($P{output_grid_dz}) && $P{output_grid_dz}!=$opt_o &&!$opt_o_override);
+ &antsInfo("WARNING: inconsistent %output_grid_minsamp in input files")
+ if (defined($P{output_grid_minsamp}) && $P{output_grid_minsamp}*2!=$opt_k &&!$opt_k_override);
+
+ $bin_length = sprintf('%d & %d',round($bin_length),($P{ADCP_bin_length}))
+ unless (round($bin_length) == round($P{ADCP_bin_length}));
+ $pulse_length = sprintf('%d & %d',round($pulse_length),round($P{ADCP_pulse_length}))
+ unless (round($pulse_length) == round($P{ADCP_pulse_length}));
+ $blanking_dist = sprintf('%d & %d',round($blanking_dist),($P{ADCP_blanking_distance}))
+ unless (round($blanking_dist) == round($P{ADCP_blanking_distance}));
+
+ $PROF = $STN = $id; $RUN = antsRequireParam('run_label');
+ undef(@dscMin); undef(@dscMax);
+ do "./EditParams";
+
+ if (@dcw1) { # 2nd file in dual-head profile => plot 1st
+ GMT_psxy('-W2,255/127/80,-');
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",median(@{$dcw1[$bi]}),($bi+0.5)*$opt_o);
+ }
+ GMT_psxy('-W2,46/139/87,-');
+ for (my($bi)=0; $bi<=$#ucw1; $bi++) {
+ printf(GMT "%f %f\n",median(@{$ucw1[$bi]}),($bi+0.5)*$opt_o);
+ }
+ undef(@dcw1); undef(@ucw1);
+ }
+ }
+ $filt++,next if &isBad();
+ $min_depth = $ants_[0][$dF] if ($ants_[0][$dF] < $min_depth);
+ $max_depth = $ants_[0][$dF] if ($ants_[0][$dF] > $max_depth);
+ my($bi) = $ants_[0][$dF]/$opt_o;
+ if ($ants_[0][$dcF]) { # downcast
+ push(@{$dcw[$bi]},$ants_[0][$wF]); # vertical velocity
+ push(@{$dcw1[$bi]},$ants_[0][$wF]) if ($dual_head);
+ push(@{$dce[$bi]},$ants_[0][$eF]); # elapsed time
+ } else { # upcast
+ push(@{$ucw[$bi]},$ants_[0][$wF]);
+ push(@{$ucw1[$bi]},$ants_[0][$wF]) if ($dual_head);
+ push(@{$uce[$bi]},$ants_[0][$eF]);
+ }
+}
+
+if (@dcw1) { # 2nd file in dual-head profile in buffer => plot it
+ GMT_psxy('-W2,255/127/80,.');
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",median(@{$dcw1[$bi]}),($bi+0.5)*$opt_o);
+ }
+ GMT_psxy('-W2,46/139/87,.');
+ for (my($bi)=0; $bi<=$#ucw1; $bi++) {
+ printf(GMT "%f %f\n",median(@{$ucw1[$bi]}),($bi+0.5)*$opt_o);
+ }
+}
+
+&antsInfo("%d measurements edited (%d%% of total)",$filt,round(100*$filt/$r))
+ if ($filt > 0);
+
+#----------------------------------------------------------------------
+# Average and Output Profile, Add to Plot
+#----------------------------------------------------------------------
+
+@antsNewLayout = ('depth','hab',
+ 'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
+ 'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp');
+
+if ($dual_head) { # selected %PARAMs only
+ &antsAddParams('profile_id',$id,'lat',$P{lat},'lon',$P{lon});
+ &antsAddParams($dayNoP,$dn,'run_label',"$first_label & $P{run_label}");
+ &antsAddParams('output_grid_dz',$opt_o,'output_grid_minsamp',$opt_k);
+ &antsAddParams('ADCP_bin_length',round($bin_length),'ADCP_pulse_length',round($pulse_length))
+ &antsAddParams('ADCP_blanking_distance',round($blanking_dist));
+ &antsAddParams('min_depth',round($min_depth),'max_depth',round($max_depth));
+ &antsAddParams('water_depth',$P{water_depth},'water_depth.sig',$P{water_depth.sig});
+ undef($antsOldHeaders);
+}
+
+&antsInfo("WARNING: unknown water depth (no height-above-bottom)")
+ unless numberp($P{water_depth});
+
+my(@dcwm,@ucwm,@dcns,@ucns,@dcwmad,@ucwmad);
+for (my($bi)=0; $bi<=max($#dcw,$#ucw); $bi++) {
+ $dcwm[$bi] = median(@{$dcw[$bi]});
+ $ucwm[$bi] = median(@{$ucw[$bi]});
+ $dcns[$bi] = @{$dcw[$bi]};
+ $ucns[$bi] = @{$ucw[$bi]};
+ $dcwmad[$bi] = mad2($dcwm[$bi],@{$dcw[$bi]});
+ $ucwmad[$bi] = mad2($ucwm[$bi],@{$ucw[$bi]});
+ push(@{$out[$bi]},
+ ($bi+0.5)*$opt_o,
+ (numberp($P{water_depth}) ? $P{water_depth}-($bi+0.5)*$opt_o : nan),
+ avg(@{$dce[$bi]}),
+ (($dcns[$bi]>=$opt_k)?$dcwm[$bi]:nan),(($dcns[$bi]>=$opt_k)?$dcmad[$bi]:nan),
+ scalar(@{$dcw[$bi]}),
+ avg(@{$uce[$bi]}),
+ (($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan),(($ucns[$bi]>=$opt_k)?$ucmad[$bi]:nan),
+ scalar(@{$ucw[$bi]}));
+ &antsOut(@{$out[$bi]});
+}
+
+if ($opt_p) { # complete plot
+ GMT_setR($R);
+
+ GMT_psxy('-W4,255/127/80'); # median profiles
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",$dcwm[$bi],($bi+0.5)*$opt_o);
+ }
+ GMT_psxy('-W4,46/139/87');
+ for (my($bi)=0; $bi<=$#ucw1; $bi++) {
+ printf(GMT "%f %f\n",$ucwm[$bi],($bi+0.5)*$opt_o);
+ }
+
+ GMT_psxy('-Sc0.1 -G255/127/80'); # m.a.d. profiles
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",$dcwmad[$bi],($bi+0.5)*$opt_o);
+ }
+ GMT_psxy('-Sc0.1 -G46/139/87');
+ for (my($bi)=0; $bi<=$#ucw1; $bi++) {
+ printf(GMT "%f %f\n",$ucwmad[$bi],($bi+0.5)*$opt_o);
+ }
+
+ GMT_setR($R2);
+ GMT_psxy('-Mn -W1,255/127/80');
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",$dcns[$bi],($bi+0.5)*$opt_o);
+ }
+ GMT_psxy('-Mn -W1,46/139/87');
+ for (my($bi)=0; $bi<=$#dcw1; $bi++) {
+ printf(GMT "%f %f\n",$ucns[$bi],($bi+0.5)*$opt_o);
+ }
+
+ GMT_psbasemap('-Bf10a1000-950:" # of Samples":N');
+ GMT_psbasemap('-Ba1000-900N'); GMT_psbasemap('-Ba1000-850N');
+
+ $depth_tics = ($ymax-$ymin< 1000) ? 'f10a100' : 'f100a500';
+ GMT_setR($R);
+ GMT_psbasemap('-Bf0.01a10-10.05:"Vertical Velocity [m/s] ":/' .
+ $depth_tics . ':"Depth [m]":WeS');
+ GMT_psbasemap('-Ba10-9.95S'); GMT_psbasemap('-Ba10-9.85S');
+
+ GMT_setR('-R0/1/0/1');
+ GMT_pstext('-Gblue -N');
+ if (defined($outfile)) { print(GMT "0.01 -0.06 14 0 0 TL $outfile\n"); }
+ else { printf(GMT "0.01 -0.06 14 0 0 TL %03d\n",$id); }
+ GMT_pstext();
+ print(GMT '0.62 0.98 12 0 0 MR m.a.d.');
+ GMT_end();
+}
+
+&antsExit(0);
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/LADCP_wspec Sun Jul 26 20:04:48 2015 +0000
@@ -0,0 +1,354 @@
+#!/usr/bin/perl
+#======================================================================
+# L A D C P _ W S P E C
+# doc: Thu Jun 11 12:02:49 2015
+# dlm: Tue Jun 16 13:20:57 2015
+# (c) 2012 A.M. Thurnherr
+# uE-Info: 334 19 NIL 0 0 72 10 2 4 NIL ofnI
+#======================================================================
+
+$antsSummary = 'calculate VKE window spectra from LADCP profiles';
+$antsMinLibVersion = 6.1;
+
+# HISTORY:
+# Jun 11, 2015: - adapted from [binpgrams]
+# Jun 12, 2015: - renamed %PARAM prefixes
+# Jun 15, 2015: - BUG: de-meaning did not respect _gap variables
+# - added %output_depth_resolution
+# - reversed semantics of -d/-u
+# Jun 16, 2015: - reversed semantics of -t
+# - re-added pwrdens.0 to make output consistent with [binpgrams]
+
+($ANTSBIN) = (`which list` =~ m{^(.*)/[^/]*$});
+($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
+
+require "$ANTSLIB/ants.pl";
+require "$ANTSLIB/antsfilters.pl";
+require "$ANTSLIB/libstats.pl";
+require "$ANTSLIB/fft.pl";
+require "$ANTSLIB/lfit.pl";
+require "$ANTSLIB/nrutil.pl";
+require "$ANTSBIN/.lsfit.poly";
+require "$ANTSBIN/.nminterp.linear";
+
+#----------------------------------------------------------------------
+# Usage
+#----------------------------------------------------------------------
+
+&antsUsage('bc:dg:o:s:tuw:',0,
+ '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]',
+ '[-d)own/-u)pcast-only] [exclude -b)ottom window]',
+ '[shortwave -c)utoff <kz or lambda>]',
+ '[-s)urface <layer depth to exclude[150m]>',
+ '[-g)ap <max depth layer to fill with interpolation[40m]>]',
+ '[-w)indow <power-of-two input-records[32]>]',
+ '[LADCP-profile(s)]');
+
+&antsIntOpt(\$opt_o,0); # polynomial order to remove
+if ($opt_o >= 0) { # init model
+ &modelUsage();
+ matrix(\@covar,1,$modelNFit,1,$modelNFit);
+ vector(\@afunc,1,$modelNFit);
+ &antsAddParams('wspec::demean_poly_order',$opt_o);
+}
+
+croak("$0: cannot ignore both down- and upcast\n")
+ unless ($opt_d+$opt_u < 2);
+if ($opt_d) {
+ &antsAddParams('wspec::input_data','dc');
+} elsif ($opt_u) {
+ &antsAddParams('wspec::input_data','uc');
+} else {
+ &antsAddParams('wspec::input_data','dc/uc');
+}
+&antsAddParams('wspec::cos_taper_applied',$opt_t ? 'no' : 'yes');
+&antsAddParams('wspec::btm_window_included',$opt_b ? 'no' : 'yes');
+
+if (defined($opt_c)) { # shortwave cutoff
+ $kzlim = ($opt_c < 1) ? $opt_c : 2*$PI/$opt_c;
+ &antsAddParams('wspec::shortwave_cutoff',$kzlim);
+} else {
+ $kzlim = 9e99;
+}
+
+&antsCardOpt($opt_w); # window size
+
+&antsCardOpt(\$opt_g,40); # gap length [m]
+&antsAddParams('wspec::min_gap_thickness',$opt_g);
+
+&antsCardOpt(\$opt_s,150); # surface layer
+&antsAddParams('wspec::surface_layer',$opt_s);
+
+&ISUsage; # interpolation model
+
+#----------------------------------------------------------------------
+# Read Data & Define Layout
+#----------------------------------------------------------------------
+
+croak("LADCP_VKE: spectral-input mode must be selected manually (-f)\n")
+ if defined(fnrNoErr('pwrdens.0')) && !defined(fnrNoErr('dc_w'));
+
+$zfnr = fnr('depth'); # required fields
+unless (defined(fnrNoErr('dc_w'))) { #
+}
+$dcwfnr = fnr('dc_w');
+$ucwfnr = fnr('uc_w');
+$habfnr = fnrNoErr('hab'); # optional fields
+
+&antsInstallBufFull(0); # read entire file
+&antsIn();
+
+while ($ants_[0][$zfnr] < $opt_s) { # remove surface layer
+ shift(@ants_);
+}
+
+for ($trimmed=0; !numberp($ants_[$#ants_][$dcwfnr]) && !numberp($ants_[$#ants_][$ucwfnr]); $trimmed++) {
+ pop(@ants_);
+}
+&antsInfo("$trimmed trailing non-numeric records trimmed")
+ if ($trimmed);
+
+$dz = &antsXCheck($zfnr,0,$#ants_,1.01); # calc dT; 1% jitter
+&antsAddParams("wspec::input_depth_resolution",$dz);
+
+$opt_g = int(($opt_g - 1) / $dz); # [m] -> [records]
+
+unless (defined($opt_w)) { # window size
+ for ($opt_w=32; $opt_w*$dz>600; $opt_w/=2) {}
+ &antsInfo("%d-m windows ($opt_w records)",$opt_w*$dz);
+}
+&antsAddParams('wspec::window_size',$opt_w,'wspec::output_depth_resolution',$dz*$opt_w);
+
+croak(sprintf("$0: insufficient data (%d records found, %d required)\n",scalar(@ants_),$opt_w))
+ unless (@ants_ >= $opt_w);
+
+$zrange = $opt_w * $dz; # NB: not equal to max-min!!!
+$resolution_bandwidth = 1 / $zrange;
+$resolution_bandwidth *= 2*$PI;
+&antsAddParams('wspec::resolution_bandwidth',$resolution_bandwidth);
+
+push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','nsamp');
+for (my($i)=0; $i<$opt_w/2+1; $i++) {
+ my($kz) = 2*$PI*$i/$zrange;
+ last if ($kz > $kzlim);
+ &antsAddParams(sprintf('wspec::k.%d',$i),$kz);
+ &antsAddParams(sprintf('wspec::lambda.%d',$i),$i ? $zrange/$i : inf);
+ push(@antsNewLayout,sprintf('pwrdens.%d',$i));
+}
+push(@antsNewLayout,'pwr.tot');
+
+&antsActivateOut();
+
+#----------------------------------------------------------------------
+# interpolate short gaps
+#----------------------------------------------------------------------
+
+&ISInit($dcwfnr,$zfnr);
+&ISInit($ucwfnr,$zfnr);
+
+my($dcLastValid,$ucLastValid,$interp);
+for (my($r)=0; $r<=$#ants_; $r++) {
+ if (numberp($ants_[$r][$dcwfnr])) { # number => no interp.
+ $dcLastValid = $r;
+ } elsif (defined($dcLastValid)) {
+ $ants_[$r][$dcwfnr] = &interpolate($zfnr,$opt_g,$dcwfnr,$dcLastValid,$r);
+ $interp++ if numberp($ants_[$r][$dcwfnr]);
+ }
+ if (numberp($ants_[$r][$ucwfnr])) { # number => no interp.
+ $ucLastValid = $r;
+ } elsif (defined($ucLastValid)) {
+ $ants_[$r][$ucwfnr] = &interpolate($zfnr,$opt_g,$ucwfnr,$ucLastValid,$r);
+ $interp++ if numberp($ants_[$r][$ucwfnr]);
+ }
+}
+
+&antsInfo("$interp non-numeric values interpolated")
+ if ($interp);
+
+#----------------------------------------------------------------------
+# loop over windows
+#----------------------------------------------------------------------
+
+sub avgF($$) # average field over window
+{
+ my($f,$r) = @_;
+ my(@vals);
+
+ push(@vals,$ants_[$r++][$f])
+ while (@vals < $opt_w);
+ return avg(@vals);
+}
+
+unless ($opt_t) { # compile taper function only if needed
+ sub cosTaperWeight($)
+ {
+ my($z) = $_[0] - $ants_[$fromR][$zfnr]; # elapsed time
+ return ($z<0.1*$zrange || $z>0.9*$zrange) ?
+ 0.5*(1+cos(10*$PI*$z/$zrange-$PI)) : 1;
+ }
+}
+
+WINDOW: for (my($widx)=1; 1; $widx++) {
+ undef(@out);
+ $out[0] = $widx;
+ local($fromR) = round(($widx-1)*($opt_w/2)); # local, cuz it's used in cosTaperWeight
+
+ #----------------------
+ # partial bottom window
+ #----------------------
+
+ if ($fromR+$opt_w-1 > $#ants_) {
+ last if ($opt_b);
+ $out[0] = -1;
+ $fromR = @ants_ - $opt_w;
+ $opt_b = 1;
+ }
+
+ #-----------------------------------
+ # output nan on non-numeric y values
+ #-----------------------------------
+
+ my($dc_gap) = $opt_u; # exclude dc with -d, uc with -u
+ my($uc_gap) = $opt_d;
+ for (my($r)=$fromR; $r<$fromR+$opt_w; $r++) {
+ $dc_gap |= !numberp($ants_[$r][$dcwfnr]);
+ $uc_gap |= !numberp($ants_[$r][$ucwfnr]);
+ }
+
+ if ($dc_gap && $uc_gap) {
+ push(@out,$ants_[$fromR+$opt_w/2][$zfnr]);
+ if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) {
+ push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);
+ push(@out,$ants_[$fromR][$zfnr]);
+ } else {
+ push(@out,$ants_[$fromR][$zfnr]);
+ push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);
+ }
+ push(@out,defined($habfnr) ?
+ avgF($habfnr,$fromR) : nan);
+ push(@out,nan);
+ for ($i=0; $i<=$opt_w/2; $i++) {
+ push(@out,nan);
+ }
+ &antsOut(@out); # output nan record and go to next window
+ next WINDOW;
+ }
+
+ #--------------------
+ # save current values
+ #--------------------
+
+ for (my($i)=0; $i<$opt_w; $i++) {
+ $dcwbuf[$i] = $ants_[$fromR+$i][$dcwfnr];
+ $ucwbuf[$i] = $ants_[$fromR+$i][$ucwfnr];
+ }
+
+ #------------------------
+ # polynomial de-"mean"ing
+ #------------------------
+
+ if ($opt_o >= 0) {
+ my($calc);
+ unless ($dc_gap) {
+ &modelInit(); # dc
+ for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; }
+ &lfit($zfnr,$dcwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1);
+ &modelCleanup();
+ for (my($i)=0; $i<$opt_w; $i++) {
+ modelEvaluate($i,$zfnr,\@afunc);
+ for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) {
+ $calc += $A[$p] * $afunc[$p];
+ }
+ $ants_[$fromR+$i][$dcwfnr] -= $calc;
+ }
+ }
+ unless ($uc_gap) {
+ &modelInit(); # uc
+ for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; }
+ &lfit($zfnr,$ucwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1);
+ &modelCleanup();
+ for (my($i)=0; $i<$opt_w; $i++) {
+ modelEvaluate($i,$zfnr,\@afunc);
+ for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) {
+ $calc += $A[$p] * $afunc[$p];
+ }
+ $ants_[$fromR+$i][$ucwfnr] -= $calc;
+ }
+ }
+ }
+
+ #-----------
+ # taper data
+ #-----------
+
+ unless ($opt_t) {
+ for (my($i)=0; $i<$opt_w; $i++) {
+ $ants_[$fromR+$i][$dcwfnr] *= &cosTaperWeight($ants_[$fromR+$i][$zfnr]);
+ $ants_[$fromR+$i][$ucwfnr] *= &cosTaperWeight($ants_[$fromR+$i][$zfnr]);
+ }
+ $taper_correction = 1/0.875;
+ } else {
+ $taper_correction = 1;
+ }
+
+ #-------------
+ # PSD Estimate
+ #-------------
+
+# for (my($r)=0; $r<$opt_w; $r++) {
+# print(STDERR "$ants_[$fromR+$r][$dcwfnr], $ants_[$fromR+$r][$ucwfnr]\n");
+# }
+
+ @dc_coeff = &cFFT($dcwfnr,nan,$opt_w,$fromR) # FFT
+ unless ($dc_gap);
+ @uc_coeff = &cFFT($ucwfnr,nan,$opt_w,$fromR)
+ unless ($uc_gap);
+ croak("$0: -n $opt_w not a power-of-two\n")
+ unless (@dc_coeff/2==$opt_w || @uc_coeff/2==$opt_w);
+
+ @dc_pwr = &pgram_onesided($opt_w,@dc_coeff) # total power
+ unless ($dc_gap);
+ @uc_pwr = &pgram_onesided($opt_w,@uc_coeff)
+ unless ($uc_gap);
+
+ push(@out,$ants_[$fromR+$opt_w/2][$zfnr]); # middle z
+ if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) { # input descending
+ push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); # z.min
+ push(@out,$ants_[$fromR][$zfnr]); # z.max
+ } else { # input ascending
+ push(@out,$ants_[$fromR][$zfnr]); # z.min
+ push(@out,$ants_[$fromR+$opt_w-1][$zfnr]); # z.max
+ }
+ push(@out,defined($habfnr) ? # hab
+ avgF($habfnr,$fromR) : nan);
+ my($nsamp) = !$dc_gap + !$uc_gap; # nsamp
+ push(@out,$nsamp);
+
+ my($totP) = 0; # power
+ my($i);
+ for ($i=0; $i<$opt_w/2+1; $i++) {
+ my($sumP) = 0;
+ $sumP += $dc_pwr[$i] * $taper_correction unless ($dc_gap);
+ $sumP += $uc_pwr[$i] * $taper_correction unless ($uc_gap);
+ push(@out,$sumP/$nsamp/$resolution_bandwidth)
+ unless (antsParam("k.$i") > $kzlim);
+ $totP += $sumP;
+ }
+ push(@out,$totP); # total power
+
+ &antsOut(@out);
+
+ #--------------------------------
+ # undo tapering and/or de-meaning
+ #--------------------------------
+
+ for ($i=0; $i<$opt_w; $i++) {
+ $ants_[$fromR+$i][$dcwfnr] = $dcwbuf[$i];
+ $ants_[$fromR+$i][$ucwfnr] = $ucwbuf[$i];
+ }
+}
+
+&antsInfo("$skipped windows with non-numeric values skipped")
+ if ($skipped);
+
+antsExit();
--- a/LWplot_BR Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_BR Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ B R
# doc: Mon Oct 17 10:57:12 2011
-# dlm: Thu Apr 16 10:40:35 2015
+# dlm: Wed May 20 20:43:04 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 19 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 20 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -16,6 +16,7 @@
# Apr 5, 2015: - made fixbb optional
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -30,13 +31,14 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$max_bin" ] && max_bin=`expr -- "$line" : '#ANTS#PARAMS#.*max_bin{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
--- a/LWplot_CAE Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_CAE Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ C A E
# doc: Wed May 15 19:35:58 2013
-# dlm: Thu Apr 16 10:40:42 2015
+# dlm: Wed May 20 20:43:33 2015
# (c) 2013 A.M. Thurnherr
-# uE-Info: 18 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 40 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -15,6 +15,7 @@
# Apr 5, 2015: - made fixbb optional
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -29,12 +30,13 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
--- a/LWplot_Sv Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_Sv Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ S V
# doc: Sat Oct 15 13:42:50 2011
-# dlm: Thu Apr 16 11:18:10 2015
+# dlm: Wed May 20 20:43:55 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 93 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 52 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -19,6 +19,8 @@
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
# - added seabed if available
+# Apr 20, 2015: - removed _skel dependence
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -34,9 +36,9 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
@@ -46,6 +48,7 @@
[ -z "$max_ens" ] && max_ens=`expr -- "$line" : '#ANTS#PARAMS#.*max_ens{\([^}]*\)}`
[ -z "$ADCP_bin_length" ] && ADCP_bin_length=`expr -- "$line" : '#ANTS#PARAMS#.*ADCP_bin_length{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
@@ -102,9 +105,10 @@
ens_tics=f500a1000
fi
-cat `which LWplot_Sv | sed 's@LWplot_Sv$@Sv_scale.skel@'` >> "$eps_file"
+psbasemap -O -K $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
-psbasemap -O $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
+gmtset ANNOT_FONT_SIZE_PRIMARY 7
+psscale -O -E -D8/2/3/0.4 $C -B/:S@-v@-: >> "$eps_file"
rm $TMPFILE
[ -n "`which fixbb`" ] && fixbb "$eps_file"
--- a/LWplot_TL Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_TL Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ T L
# doc: Thu Oct 13 10:51:49 2011
-# dlm: Thu Apr 16 10:40:56 2015
+# dlm: Wed May 20 20:44:17 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 27 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 46 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -24,15 +24,16 @@
# - BUG: ps file had not been closed properly (fixbb took care of that)
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
+# May 20, 2015: - made it quit on EOF
USAGE="Usage: $0 <eps-file> [in-file]"
[ $# -eq 2 ] && exec <"$2" "$0" "$1"
[ $# -ne 1 ] && { echo $USAGE >&2; exit 1; }
eps_file="$1"
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_elapsed" ] && min_elapsed=`expr -- "$line" : '#ANTS#PARAMS#.*elapsed.min{\([^}]*\)}`
@@ -41,6 +42,7 @@
[ -z "$best_scan_offsets" ] && best_scan_offsets=`expr -- "$line" : '#ANTS#PARAMS#.*best_scan_offsets{\([^}]*\)}`
[ -z "$to_elapsed_limits" ] && to_elapsed_limits=`expr -- "$line" : '#ANTS#PARAMS#.*to_elapsed_limits{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
--- a/LWplot_corr Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_corr Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ C O R R
# doc: Sat Oct 15 13:42:50 2011
-# dlm: Thu Apr 16 10:41:02 2015
+# dlm: Wed May 20 20:44:44 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 20 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 48 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -17,6 +17,8 @@
# Apr 5, 2015: - made fixbb optional
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
+# Apr 20, 2015: - removed _skel dependence
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -31,9 +33,9 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
@@ -42,6 +44,7 @@
[ -z "$max_ens" ] && max_ens=`expr -- "$line" : '#ANTS#PARAMS#.*max_ens{\([^}]*\)}`
[ -z "$ADCP_bin_length" ] && ADCP_bin_length=`expr -- "$line" : '#ANTS#PARAMS#.*ADCP_bin_length{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
@@ -90,9 +93,10 @@
ens_tics=f500a1000
fi
-cat `which LWplot_corr | sed 's@LWplot_corr$@corr_scale.skel@'` >> "$eps_file"
+psbasemap -O -K $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
-psbasemap -O $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
+gmtset ANNOT_FONT_SIZE_PRIMARY 7
+psscale -O -E -D8/2/3/0.4 $C -B/:corr: >> "$eps_file"
rm $TMPFILE
[ -n "`which fixbb`" ] && fixbb "$eps_file"
--- a/LWplot_prof Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_prof Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ P R O F
# doc: Fri Oct 14 09:42:36 2011
-# dlm: Thu Apr 16 10:41:09 2015
+# dlm: Wed May 20 20:45:55 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 24 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 56 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -21,6 +21,12 @@
# Apr 5, 2015: - made fixbb optional
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
+# Apr 26, 2015: - adapted to LADCP_w_regrid output
+# - removed -N
+# - moved run_label above plot area
+# May 19, 2015: - adapted to new layout (hab)
+# - BUG: erroneously let .VKE profiles pass
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -35,21 +41,24 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$done" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
[ -z "$max_depth" ] && max_depth=`expr -- "$line" : '#ANTS#PARAMS#.*max_depth{\([^}]*\)}`
- [ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ [ -z "$water_depth" ] && water_depth=`expr -- "$line" : '#ANTS#PARAMS#.*water_depth{\([^}]*\)}`
+ newfields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ [ -z "$newfields" ] || fields=$newfields
+ done=`expr -- "$line" : '^[^#]\(.*\)'`
+ read line
done
set -- $fields
[ "$1" = depth -a "$4" = dc_w -a "$5" = dc_w.mad -a "$6" = dc_w.nsamp -a \
- "${11}" = uc_w -a "${12}" = uc_w.mad -a "${13}" = uc_w.nsamp -a \
- "${20}" = BT_w -a "${21}" = BT_w.mad -a "${22}" = BT_w.nsamp ] || {
- echo "$0: file layout error" >&2
+ "$8" = uc_w -a "$9" = uc_w.mad -a "${10}" = uc_w.nsamp ] || {
+ echo "$0: file layout error ($1, $4, $5, $6, $8, $9, ${10})" >&2
exit 1
}
@@ -62,7 +71,10 @@
cd /tmp/$$
TMPFILE=/tmp/$$.LWplot_prof
-cat > $TMPFILE
+{
+ echo $line
+ sed s/#.*//
+} > $TMPFILE
[ -f .gmtdefaults4 ] ||
gmtset PAPER_MEDIA letter+ \
@@ -70,8 +82,16 @@
WANT_EURO_FONT true \
PLOT_DEGREE_FORMAT ddd:mm:ssF
-R=-R-0.07/0.35/$min_depth/$max_depth
-R2=-R-200/200/$min_depth/$max_depth
+[ -n "$water_depth" ] && blim=`echo "scale=0;$water_depth/1+25"|bc` \
+ || blim=`echo "scale=0;$max_depth/1+25"|bc`
+
+XTICS="-0.05 0.05 0.15"
+XMIN=-0.1
+
+[ -r ./.LWplot_prof ] && . ./.LWplot_prof
+
+R=-R$XMIN/0.35/0/$blim
+R2=-R-200/200/0/$blim
U=-R0/1/0/1
X=-JX10/-10
@@ -83,22 +103,27 @@
{ echo 150 $min_depth; echo 150 $max_depth; } | psxy -O -K $R2 $X >> "$eps_file"
# VERTICAL VELOCITIES
-awk '{print $4, $1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4/coral >> "$eps_file"
-awk '{print $11,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4/SeaGreen >> "$eps_file"
-awk '{print $20,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4/black >> "$eps_file"
+awk '{print $4,$1}' $TMPFILE | psxy -O -K -Mn $R $X -W4/coral >> "$eps_file"
+awk '{print $8,$1}' $TMPFILE | psxy -O -K -Mn $R $X -W4/SeaGreen >> "$eps_file"
-# MEAN ABSOLUTE DEVIATIONS (COMBINED SOLUTION)
-awk '{print $5,$1, $4}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gcoral >> "$eps_file"
-awk '{print $12,$1,$11}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -GSeaGreen >> "$eps_file"
-awk '{print $21,$1,$20}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gblack >> "$eps_file"
+# MEAN ABSOLUTE DEVIATIONS
+awk '{print $5,$1,$4}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -Sc0.1c -Gcoral >> "$eps_file"
+awk '{print $9,$1,$8}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -Sc0.1c -GSeaGreen >> "$eps_file"
+
+# NUMBER OF SAMPLES
+awk '{print $6,$1,$5}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -Mn $R2 $X -W1/coral >> "$eps_file"
+awk '{print $10,$1,$9}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -Mn $R2 $X -W1/SeaGreen >> "$eps_file"
-# NUMBER OF SAMPLES (COMBINED SOLUTION)
-awk '{print $6,$1, $4}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/coral >> "$eps_file"
-awk '{print $13,$1,$11}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/SeaGreen >> "$eps_file"
-awk '{print $22,$1,$20}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/black >> "$eps_file"
+# SEABED
+[ -n "$water_depth" ] && {
+ echo $XMIN $blim;
+ echo 0.35 $blim;
+ echo 0.35 $water_depth;
+ echo $XMIN $water_depth;
+} | psxy -O -K $R $X -G204/153/102 >> "$eps_file"
# LABELS
-echo 0.02 0.02 12 0 0 TL $out_basename $run_label | pstext -O -K $U $X >> "$eps_file"
+echo 0.02 -0.05 10 0 0 TL $out_basename $run_label | pstext -O -K -Gblue -N $U $X >> "$eps_file"
echo 0.6 0.98 12 0 0 BR m.a.d. | pstext -O -K $U $X >> "$eps_file"
if [ 0 -eq `echo "($max_depth-$min_depth)>1000"|bc` ]
--- a/LWplot_prof_2beam Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_prof_2beam Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ P R O F _ 2 B E A M
# doc: Fri Oct 14 09:42:36 2011
-# dlm: Thu Apr 16 13:59:55 2015
+# dlm: Wed May 20 20:42:08 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 1 9 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 21 39 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -16,6 +16,16 @@
# - changed shell from /bin/ksh to /bin/sh
# - slightly increased w range
# - added seabed if available
+# Apr 23, 2015: - added support for .LWplot_prof_2beam
+# - removed all -N
+# May 20, 2015: - made it quit on EOF
+
+# NOTES:
+# - In order to extend the x-axis range, create a file
+# called [.LWplot_prof_2beam] in the processing directory and
+# add something like the following lines:
+# XMIN=-0.27
+# XTICS="-0.25 -0.15 -0.05 0.05"
#--------------------------------------------------
# Usage
@@ -30,15 +40,16 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
[ -z "$max_depth" ] && max_depth=`expr -- "$line" : '#ANTS#PARAMS#.*max_depth{\([^}]*\)}`
[ -z "$water_depth" ] && water_depth=`expr -- "$line" : '#ANTS#PARAMS#.*water_depth{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
@@ -53,6 +64,19 @@
# Plot Data
#--------------------------------------------------
+[ -n "$water_depth" ] && blim=`echo "scale=0;$water_depth/1+25"|bc` \
+ || blim=`echo "scale=0;$max_depth/1+25"|bc`
+
+XTICS="-0.05 0.05 0.15"
+XMIN=-0.1
+
+[ -r ./.LWplot_prof_2beam ] && . ./.LWplot_prof_2beam
+
+R=-R$XMIN/0.35/0/$blim
+R2=-R-200/200/0/$blim
+U=-R0/1/0/1
+X=-JX10/-10
+
eps_file="$PWD/$eps_file" # make outfile name absolute (hopefully, it is not already...)
mkdir /tmp/$$ # GMT makes tmpfiles and is not reentrant
cd /tmp/$$
@@ -66,13 +90,6 @@
WANT_EURO_FONT true \
PLOT_DEGREE_FORMAT ddd:mm:ssF
-[ -n "$water_depth" ] && blim=`echo "scale=0;$water_depth/1+25"|bc` \
- || blim=`echo "scale=0;$max_depth/1+25"|bc`
-R=-R-0.1/0.35/0/$blim
-R2=-R-200/200/0/$blim
-U=-R0/1/0/1
-X=-JX10/-10
-
# FRAME
{ echo 0 0; echo 0 $blim; } | psxy -P -K $R $X > "$eps_file"
{ echo 0.07 0; echo 0.07 $blim; echo 0.18 $blim; echo 0.18 0; } | psxy -O -K $R $X -L -G200 >> "$eps_file"
@@ -81,28 +98,28 @@
{ echo 150 0; echo 150 $blim; } | psxy -O -K $R2 $X >> "$eps_file"
# VERTICAL VELOCITIES (2-BEAM SOLUTIONS)
-awk '{print $7, $1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,coral,6_2:0 >> "$eps_file"
-awk '{print $8, $1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,coral,4_6:0 >> "$eps_file"
-awk '{print $14,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,SeaGreen,6_2:0 >> "$eps_file"
-awk '{print $15,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,SeaGreen,4_6:0 >> "$eps_file"
-awk '{print $16,$1}' $TMPFILE | psxy -O -K -N -Mn $R $X -W4,black >> "$eps_file"
+awk '{print $7, $1}' $TMPFILE | psxy -O -K -Mn $R $X -W4,coral,6_2:0 >> "$eps_file"
+awk '{print $8, $1}' $TMPFILE | psxy -O -K -Mn $R $X -W4,coral,4_6:0 >> "$eps_file"
+awk '{print $14,$1}' $TMPFILE | psxy -O -K -Mn $R $X -W4,SeaGreen,6_2:0 >> "$eps_file"
+awk '{print $15,$1}' $TMPFILE | psxy -O -K -Mn $R $X -W4,SeaGreen,4_6:0 >> "$eps_file"
+awk '{print $16,$1}' $TMPFILE | psxy -O -K -Mn $R $X -W4,black >> "$eps_file"
# MEAN ABSOLUTE DEVIATIONS (COMBINED SOLUTION)
-awk '{print $5,$1, $4}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gcoral >> "$eps_file"
-awk '{print $12,$1,$11}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -GSeaGreen >> "$eps_file"
-awk '{print $17,$1,$20}' $TMPFILE | grep -vi nan | psxy -O -K -N $R $X -Sc0.1c -Gblack >> "$eps_file"
+awk '{print $5,$1, $4}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -Sc0.1c -Gcoral >> "$eps_file"
+awk '{print $12,$1,$11}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -Sc0.1c -GSeaGreen >> "$eps_file"
+awk '{print $17,$1,$20}' $TMPFILE | grep -vi nan | psxy -O -K $R $X -Sc0.1c -Gblack >> "$eps_file"
# NUMBER OF SAMPLES (COMBINED SOLUTION)
-awk '{print $6,$1, $4}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/coral >> "$eps_file"
-awk '{print $13,$1,$11}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/SeaGreen >> "$eps_file"
-awk '{print $18,$1,$20}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -N -Mn $R2 $X -W1/black >> "$eps_file"
+awk '{print $6,$1, $4}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -Mn $R2 $X -W1/coral >> "$eps_file"
+awk '{print $13,$1,$11}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -Mn $R2 $X -W1/SeaGreen >> "$eps_file"
+awk '{print $18,$1,$20}' $TMPFILE | sed '/nan/s/.*/nan/' | psxy -O -K -Mn $R2 $X -W1/black >> "$eps_file"
# SEABED
[ -n "$water_depth" ] && {
- echo -0.1 $blim;
+ echo $XMIN $blim;
echo 0.35 $blim;
echo 0.35 $water_depth;
- echo -0.1 $water_depth;
+ echo $XMIN $water_depth;
} | psxy -O -K $R $X -G204/153/102 >> "$eps_file"
# LABELS
@@ -113,9 +130,11 @@
[ $blim -lt 1000 ] && depth_tics=f10a100\
|| depth_tics=f100a500
-psbasemap -O -K $R $X -Bf0.01a10-10.05:"Vertical Velocity [m/s] ":/$depth_tics:"Depth [m]":WeS >> "$eps_file"
-psbasemap -O -K $R $X -Ba10-9.95S >> "$eps_file"
-psbasemap -O -K $R $X -Ba10-9.85S >> "$eps_file"
+psbasemap -O -K $R $X -Bf0.01:"Vertical Velocity [m/s] ":/$depth_tics:"Depth [m]":WeS >> "$eps_file"
+for t in $XTICS
+do
+ psbasemap -O -K $R $X -Ba10-`echo "10-($t)"|bc`S >> "$eps_file"
+done
psbasemap -O -K $R2 $X -Bf10a1000-950:" # of Samples":N >> "$eps_file"
psbasemap -O -K $R2 $X -Ba1000-900N >> "$eps_file"
--- a/LWplot_residuals Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_residuals Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ R E S I D U A L S
# doc: Sat Oct 15 13:42:50 2011
-# dlm: Thu Apr 16 11:20:19 2015
+# dlm: Wed May 20 20:46:41 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 91 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 50 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -19,6 +19,7 @@
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
# - added seabed if available
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -33,9 +34,9 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
@@ -45,6 +46,7 @@
[ -z "$max_ens" ] && max_ens=`expr -- "$line" : '#ANTS#PARAMS#.*max_ens{\([^}]*\)}`
[ -z "$ADCP_bin_length" ] && ADCP_bin_length=`expr -- "$line" : '#ANTS#PARAMS#.*ADCP_bin_length{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
@@ -100,9 +102,10 @@
ens_tics=f500a1000
fi
-cat `which LWplot_residuals | sed 's@LWplot_residuals$@residuals_scale.skel@'` >> "$eps_file"
+psbasemap -O -K $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
-psbasemap -O $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
+gmtset ANNOT_FONT_SIZE_PRIMARY 7
+psscale -O -E -D8/2/3/0.4 $C -B/:residuals: >> "$eps_file"
rm $TMPFILE
[ -n "`which fixbb`" ] && fixbb "$eps_file"
--- a/LWplot_spec Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_spec Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ S P E C
# doc: Thu Sep 5 18:23:41 2013
-# dlm: Thu Apr 16 10:41:27 2015
+# dlm: Wed May 20 20:47:00 2015
# (c) 2013 A.M. Thurnherr
-# uE-Info: 20 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 44 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# NB: THIS PLOTTING ROUTINE REQUIRES NON-PUBLIC ANTS TOOLS
@@ -17,6 +17,7 @@
# Apr 5, 2015: - made fixbb optional
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -31,14 +32,15 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
[ -z "$max_depth" ] && max_depth=`expr -- "$line" : '#ANTS#PARAMS#.*max_depth{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
--- a/LWplot_w Mon Apr 20 14:04:04 2015 +0000
+++ b/LWplot_w Sun Jul 26 20:04:48 2015 +0000
@@ -2,9 +2,9 @@
#======================================================================
# L W P L O T _ W
# doc: Sat Oct 15 13:42:50 2011
-# dlm: Thu Apr 16 11:17:23 2015
+# dlm: Wed May 20 20:47:23 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 76 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 52 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -20,6 +20,8 @@
# Apr 16, 2015: - removed copy of input on stdout
# - changed shell from /bin/ksh to /bin/sh
# - added seabed if available
+# Apr 20, 2015: - removed _skel dependence
+# May 20, 2015: - made it quit on EOF
#--------------------------------------------------
# Usage
@@ -34,9 +36,9 @@
# Read Header Data
#--------------------------------------------------
-while [ -z "$fields" ]
+read line
+while [ -z "$fields" -a -n "$line" ]
do
- read line
[ -z "$out_basename" ] && out_basename=`expr -- "$line" : '#ANTS#PARAMS#.*out_basename{\([^}]*\)}`
[ -z "$run_label" ] && run_label=`expr -- "$line" : '#ANTS#PARAMS#.*run_label{\([^}]*\)}`
[ -z "$min_depth" ] && min_depth=`expr -- "$line" : '#ANTS#PARAMS#.*min_depth{\([^}]*\)}`
@@ -46,6 +48,7 @@
[ -z "$max_ens" ] && max_ens=`expr -- "$line" : '#ANTS#PARAMS#.*max_ens{\([^}]*\)}`
[ -z "$ADCP_bin_length" ] && ADCP_bin_length=`expr -- "$line" : '#ANTS#PARAMS#.*ADCP_bin_length{\([^}]*\)}`
[ -z "$fields" ] && fields=`expr -- "$line" : '#ANTS#FIELDS# \(.*\)' | sed -e s/{//g -e s/}//g`
+ read line
done
set -- $fields
@@ -77,12 +80,13 @@
U=-R0/1/0/1
X=-JX10/-10
+C=-C`which LWplot_w | sed 's@LWplot_w$@w.cpt@'`
ens_width=`echo "scale=5;10/($max_ens-$min_ens+1)"|bc`
bin_length=`echo "scale=5;10*$ADCP_bin_length/($max_depth-$min_depth+$ADCP_bin_length)"|bc`
awk "{print \$1, \$4, \$7, $ens_width, $bin_length}" $TMPFILE \
- | psxy -P -K $R $X -C`which LWplot_w | sed 's@LWplot_w$@w.cpt@'` -Sr > "$eps_file"
+ | psxy -P -K $R $X $C -Sr > "$eps_file"
[ -n "$water_depth" ] && {
echo $min_ens $blim;
echo $max_ens $blim;
@@ -101,9 +105,10 @@
ens_tics=f500a1000
fi
-cat `which LWplot_w | sed 's@LWplot_w$@w_scale.skel@'` >> "$eps_file"
+psbasemap -O -K $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
-psbasemap -O $R $X -B$ens_tics:"Ensemble":/$depth_tics:"Depth [m]":WeSn >> "$eps_file"
+gmtset ANNOT_FONT_SIZE_PRIMARY 7
+psscale -O -E -D8/2/3/0.4 $C -B/:w: >> "$eps_file"
rm $TMPFILE
[ -n "`which fixbb`" ] && fixbb "$eps_file"
--- a/Makefile Mon Apr 20 14:04:04 2015 +0000
+++ b/Makefile Sun Jul 26 20:04:48 2015 +0000
@@ -1,9 +1,9 @@
#======================================================================
# M A K E F I L E
# doc: Mon Oct 17 13:29:27 2011
-# dlm: Fri Apr 17 11:26:13 2015
+# dlm: Tue Apr 21 20:55:33 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 20 36 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 19 0 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
MAKE_DIR = /Data/Makefiles
@@ -16,12 +16,8 @@
mkCPT -oc polar -- -0.03 -0.02 -0.01 -0.005 0.005 0.01 0.02 0.03 > $@
Sv.cpt:
-# mkCPT -m255/255/255 -o -- \#-90--60:2 > $@
- mkCPT -m255/255/255 -o -- \#-100--60:2 > $@
+ mkCPT -m255/255/255 -o -- \#-90--60:2 > $@
+# mkCPT -m255/255/255 -o -- \#-100--60:2 > $@
corr.cpt:
mkCPT -no -- \#70-120:5 \#120-130:0.5 > $@
-
-%_scale.skel: %.cpt
- gmtset ANNOT_FONT_SIZE_PRIMARY 7
- psscale -O -K -E -D8/2/3/0.4 -C$< -B/:"`echo $@ | sed s/_scale.skel//`": > $@
--- a/SBE_w Mon Apr 20 14:04:04 2015 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,365 +0,0 @@
-#!/usr/bin/perl
-#======================================================================
-# S B E _ W
-# doc: Mon Nov 3 17:34:19 2014
-# dlm: Fri Apr 17 10:06:21 2015
-# (c) 2014 A.M. Thurnherr
-# uE-Info: 23 55 NIL 0 0 72 2 2 4 NIL ofnI
-#======================================================================
-
-$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
-$antsMinLibVersion = 6.0;
-
-# HISTORY:
-# Nov 3, 2014: - created
-# Nov 4, 2014: - improved
-# Nov 6, 2014: - BUG: sound speed was not calculated correctly
-# - added -a
-# - added conductivity & temperature editing
-# Nov 7, 2014: - loosened outlier editing
-# - added no-valid-data error message
-# - modified binning criterion to allow any sampling
-# frequency (not just divisors of 24)
-# Apr 17, 2015: - added in-situ temperature to output
-
-($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
-
-require "$ANTS/ants.pl";
-require "$ANTS/fft.pl";
-require "$ANTS/libSBE.pl";
-require "$ANTS/libEOS83.pl";
-
-$antsParseHeader = 0; # usage
-&antsUsage('al:rs:v:',1,
- '[-v)erbosity <level[0]>]',
- '[use -a)lternate sensor pair]',
- '[-r)etain all data (no editing)]',
- '[-s)ampling <rate[6Hz]>]',
- '[-l)owpass w <cutoff[2s]>]',
- '<SBE CNV file>');
-
-&antsFloatOpt(\$opt_l,2); # defaults
-&antsCardOpt(\$opt_s,6);
-&antsCardOpt(\$opt_v,0);
-
-$CNVfile = $ARGV[0]; # open CNV file
-open(F,&antsFileArg());
-&antsActivateOut(); # activate ANTS file
-
-#----------------------------------------------------------------------
-# Read Data
-#----------------------------------------------------------------------
-
-print(STDERR "Reading $CNVfile...") if ($opt_v);
-
-($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) = # decode SBE header
- SBE_parseHeader(F,0,0); # SBE field names, no time check
-
-croak("$CNVfile: unexpected sampling interval $sampint\n")
- unless (abs($sampint-1/24) < 1e-5);
-croak("$CNVfile: unknown latitude\n")
- unless numberp($lat);
-
-$pressF = fnr('prDM');
-
-if ($opt_a) { # temp/cond alternate sensor pair
- $tempF = fnr('t190C');
- $condF = fnrNoErr('c1S/m');
- if (defined($condF)) {
- $condHistRes = 10; # 0.1 bins
- } else {
- $condF = fnr('c1mS/cm');
- $condHistRes = 1; # 1.0 bins
- }
-} else { # primary sensor pair
- $tempF = fnr('t090C');
- $condF = fnrNoErr('c0S/m');
- if (defined($condF)) {
- $condHistRes = 10;
- } else {
- $condF = fnr('c0mS/cm');
- $condHistRes = 1;
- }
-}
-
-&antsInstallBufFull(0); # read entire CNV file
-&SBEin(F,$ftype,$nfields,$nscans,$bad);
-
-printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1);
-printf(STDERR "\n") if ($opt_v);
-
-#----------------------------------------------------------------------
-# Edit Data
-# - pressure outliers & spikes
-# - conductivity outliers & spikes
-#----------------------------------------------------------------------
-
-unless ($opt_r) {
- print(STDERR "Editing Data...") if ($opt_v);
-
- #----------------------------------------
- # trim initial records with nan pressures
- #----------------------------------------
- my($trimmed) = 0; # trim leading nan pressures
- shift(@ants_),$trimmed++
- until numberp($ants_[0][$pressF]);
- printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1);
- my($lvp) = $ants_[0][$pressF];
-
- #------------------------------------------------
- # edit pressure outliers outside contiguous range
- #------------------------------------------------
- my($outliers) = 0; my($min,$max);
- for (my($s)=0; $s<$nscans; $s++) {
- $phist[$ants_[$s][$pressF]+1000]++
- if ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
- }
- for ($max=25; $phist[$max+1000]||$phist[$max+1001]; $max++) {} # outliers after 2 gaps
- for ($min=$max; $phist[$min+1000]||$phist[$min+ 999]; $min--) {}
- for (my($s)=0; $s<$nscans; $s++) {
- if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
- if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
- }
- &antsAddParams("pressure_outliers",sprintf("%d",$outliers));
- printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
-
- #----------------------------------------------------
- # edit conductivity outliers outside contiguous range
- #----------------------------------------------------
- $outliers = 0;
- my($modeSamp)=0;
- undef(@phist);
- for (my($s)=0; $s<$nscans; $s++) {
- next unless ($ants_[$s][$condF] > 0);
- my($b) = $ants_[$s][$condF]*$condHistRes; # 1/10 S/m histogram resolution (1 mS/cm)
- $phist[$b]++;
- next unless ($phist[$b] > $modeSamp);
- $modeSamp = $phist[$b]; $modeBin = $b;
- }
- for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {} # outliers after 2 gaps
- for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}
- $max /= $condHistRes; $min /= $condHistRes;
- for (my($s)=0; $s<$nscans; $s++) {
- if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }
- if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }
- }
- &antsAddParams("conductivity_outliers",sprintf("%d",$outliers));
- printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
-
- #----------------------------------------------------
- # edit temperature outliers outside contiguous range
- #----------------------------------------------------
- $outliers = 0;
- my($modeSamp)=0;
- undef(@phist);
- for (my($s)=0; $s<$nscans; $s++) {
- next unless ($ants_[$s][$tempF] > 0);
- my($b) = $ants_[$s][$tempF]*10; # 10th of a degree histogram resolution
- $phist[$b]++;
- next unless ($phist[$b] > $modeSamp);
- $modeSamp = $phist[$b]; $modeBin = $b;
- }
- for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {} # outliers after 2 gaps
- for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}
- $max /= 10; $min /= 10;
- for (my($s)=0; $s<$nscans; $s++) {
- if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }
- if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }
- }
- &antsAddParams("temperature_outliers",sprintf("%d",$outliers));
- printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);
-
- #----------------------------------------
- # edit pressure spikes based on gradients
- #----------------------------------------
-
- for (my($s)=1; $s<$nscans; $s++) { # calculated pressure gradients (across gaps)
- if (numberp($ants_[$s][$pressF])) {
- $dp[$s-1] = $ants_[$s][$pressF] - $lvp;
- $lvp = $ants_[$s][$pressF];
- } else {
- $dp[$s-1] = nan;
- }
- }
-
- my($ns1,$ns2) = (0,0);
- for (my($s)=0; $s<$nscans-2; $s++) { # consecutive large pressure gradients of opposite sign
- if (($dp[$s]*$dp[$s+1] < 0) && # tests return false if either of the dps is not defined
- (abs($dp[$s]) > 10) &&
- (abs($dp[$s+1]) > 10)) {
- $ants_[$s+1][$pressF] = nan;
- $dp[$s] = $dp[$s+1] = undef;
- $ns1++;
- }
- }
- for (my($s)=0; $s<$nscans-3; $s++) { # 3 consecutive large pressure gradients of opposite sign
- if (($dp[$s]>2 && $dp[$s+1]<-4 && $dp[$s+2]>2) ||
- ($dp[$s]<-2 && $dp[$s+1]>4 && $dp[$s+2]<-2)) {
- $ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan;
- $dp[$s] = $dp[$s+1] = $dp[$s+2] = undef;
- $ns2+=2;
- }
- }
- &antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));
- printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1);
-
- printf(STDERR "\n") if ($opt_v);
-
-} # if $opt_r
-
-#----------------------------------------------------------------------
-# Correcting for pressure bias
-#----------------------------------------------------------------------
-
-print(STDERR "Correcting for pressure bias...") if ($opt_v);
-
-my($minP) = 9e99;
-for (my($s)=0; $s<$nscans; $s++) {
- $minP = $ants_[$s][$pressF]
- if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);
-}
-croak("$CNVfile: no valid CTD pressure data below 25dbar\n")
- unless ($minP < 9e99);
-&antsAddParams('pressure_bias',$minP);
-printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1);
-for (my($s)=0; $s<$nscans; $s++) {
- $ants_[$s][$pressF] -= $minP
- if numberp($ants_[$s][$pressF]);
-}
-printf(STDERR "\n") if ($opt_v);
-
-#----------------------------------------------------------------------
-# Binning data
-#----------------------------------------------------------------------
-
-my($sps) = round(1 / $sampint / $opt_s);
-print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);
-&antsAddParams('sampling_frequency',1/$opt_s);
-&antsAddParams('sampling_rate',$opt_s);
-
-my(@press,@temp,@cond);
-my($sp,$np,$st,$nt,$sc,$nc);
-
-$sp = $st = $sc = $np = $nt = $nc = 0;
-for (my($rec)=1,my($s)=0; $s<$nscans; $s++) {
- if ($s*$sampint > $rec/$opt_s) {
- $rec++;
- push(@press,$np>0?$sp/$np:nan);
- push(@temp, $nt>0?$st/$nt:nan);
- push(@cond, $nc>0?$sc/$nc:nan);
- $sp = $st = $sc = $np = $nt = $nc = 0;
- }
- $sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]);
- $st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]);
- $sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]);
-}
-
-printf(STDERR "\n") if ($opt_v);
-
-#----------------------------------------------------------------------
-# Calculating derived quantities
-#----------------------------------------------------------------------
-
-print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v);
-
-for (my($r)=0; $r<@press; $r++) {
- $elapsed[$r] = $r/$opt_s;
- $depth[$r] = &depth($press[$r],$lat);
- croak(sprintf("$CNVfile: unrealistic depth %d m at elapsed=%.1f s (r=$r)\n",
- $depth[$r],$elapsed[$r]))
- if numberp($depth[$r]) && ($depth[$r]<0 || $depth[$r]>6100);
- $sspd[$r] = &sVel(&salin($cond[$r],$temp[$r],$press[$r]),$temp[$r],$press[$r]);
- croak(sprintf("$CNVfile: unrealistic soundspeed %dm/s at elapsed=%.1fs & depth=%.1fm ($cond[$r],$temp[$r],$press[$r])\n",
- $sspd[$r],$elapsed[$r],$depth[$r]))
- if numberp($sspd[$r]) && ($sspd[$r]<1400 || $sspd[$r]>1600);
-}
-
-$w[0] = nan;
-for (my($r)=1; $r<@depth-1; $r++) {
- $w[$r] = numbersp($depth[$r-1],$depth[$r+1])
- ? ($depth[$r+1] - $depth[$r-1]) * $opt_s
- : nan;
-}
-push(@w,nan);
-
-printf(STDERR "\n") if ($opt_v);
-
-#----------------------------------------------------------------------
-# Low-pass filter velocity data
-# - interpolate missing vertical velocities first
-#----------------------------------------------------------------------
-
-if ($opt_l > 0) {
- print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);
- &antsAddParams('w_lowpass_cutoff',$opt_l);
-
- my($trimmed) = 0;
- shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++
- until numberp($w[0]);
- my($interpolated) = 0;
- for ($r=1; $r<@w; $r++) {
- next if numberp($w[$r]);
- my($lv) = $r-1;
- for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {}
- if ($nv < @depth) {
- while ($r < $nv) {
- $w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]);
- $interpolated++;
- $r++;
- }
-
- } else {
- $trimmed += @w-$r;
- splice(@w,$r); splice(@depth,$r);
- splice(@elapsed,$r); splice(@sspd,$r);
- }
- }
- &antsAddParams('w_interpolated',$interpolated);
- printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1);
-
-
-#--------------------
-# Zero Pad Data
-#--------------------
-
- for ($pot=1; $pot<@w; $pot<<=1) {} # determine power of two
-
- for ($r=0; $r<@w; $r++) { # copy data
- $fftbuf[2*$r] = $w[$r];
- $fftbuf[2*$r+1] = 0;
- }
- printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1);
- while ($r < $pot) { # pad with zeroes
- $fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;
- $r++;
- }
-
-#--------------------
-# Low-Pass Filter
-#--------------------
-
- @fco = &FOUR1(-1,@fftbuf); # forward FFT
- $n = @fco/2;
- for (my($ip)=2; $ip<=$n; $ip+=2) { # +ve freq fco
- my($in) = 2*$n-$ip; # -ve freq fco
- my($f) = $ip/2/$n*$opt_s; # frequency
- $fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0
- if ($f > 1/$opt_l); # low-pass filter
- }
- @w_lp = &FOUR1(1,@fco); # inverse FFT
-
- printf(STDERR "\n") if ($opt_v);
-} else {
- @w_lp = @w;
-}
-
-#----------------------------------------------------------------------
-
-print(STDERR "Writing output...\n") if ($opt_v);
-
-@antsNewLayout = ('elapsed','depth','temp','sspd','w.raw','w');
-for ($r=0; $r<@w; $r++) {
- &antsOut($elapsed[$r],$depth[$r],$temp[$r],$sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp);
-}
-
-exit(0); # don't flush @ants_
--- a/Sv.cpt Mon Apr 20 14:04:04 2015 +0000
+++ b/Sv.cpt Sun Jul 26 20:04:48 2015 +0000
@@ -1,25 +1,20 @@
-#ANTS# [04/17/15 ...vertical_velocity] mkCPT '-m255/255/255' '-o' '--' '#-100--60:2'
+#ANTS# [04/21/15 ...vertical_velocity] mkCPT '-m255/255/255' '-o' '--' '#-90--60:2'
#ANTS#FIELDS# {from_z} {from_R} {from_G} {from_B} {to_z} {to_R} {to_G} {to_B}
--100 194 0 255 -98 194 0 255
--98 134 0 255 -96 134 0 255
--96 73 0 255 -94 73 0 255
--94 12 0 255 -92 12 0 255
--92 0 49 255 -90 0 49 255
--90 0 109 255 -88 0 109 255
--88 0 170 255 -86 0 170 255
--86 0 231 255 -84 0 231 255
--84 0 255 219 -82 0 255 219
--82 0 255 158 -80 0 255 158
--80 0 255 97 -78 0 255 97
--78 0 255 36 -76 0 255 36
--76 24 255 0 -74 24 255 0
--74 85 255 0 -72 85 255 0
--72 146 255 0 -70 146 255 0
--70 206 255 0 -68 206 255 0
--68 255 243 0 -66 255 243 0
--66 255 182 0 -64 255 182 0
--64 255 121 0 -62 255 121 0
--62 255 61 0 -60 255 61 0
+-90 175 0 255 -88 175 0 255
+-88 96 0 255 -86 96 0 255
+-86 16 0 255 -84 16 0 255
+-84 0 64 255 -82 0 64 255
+-82 0 143 255 -80 0 143 255
+-80 0 223 255 -78 0 223 255
+-78 0 255 207 -76 0 255 207
+-76 0 255 128 -74 0 255 128
+-74 0 255 48 -72 0 255 48
+-72 32 255 0 -70 32 255 0
+-70 112 255 0 -68 112 255 0
+-68 191 255 0 -66 191 255 0
+-66 255 239 0 -64 255 239 0
+-64 255 159 0 -62 255 159 0
+-62 255 80 0 -60 255 80 0
B 255 0 255
F 255 0 0
N 255 255 255
--- a/Sv_scale.skel Mon Apr 20 14:04:04 2015 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,144 +0,0 @@
-PSL_font_encode 0 get 0 eq {ISOLatin1+_Encoding /Helvetica /Helvetica PSL_reencode PSL_font_encode 0 1 put} if
-0 setlinecap
-0 setlinejoin
-3.32551 setmiterlimit
-944.882 59.0551 T
-
-%% PostScript produced by:
-%%GMT: psscale -O -K -E -D8/2/3/0.4 -CSv.cpt -B/:Sv:
-%%PROJ: xy -100.00000000 -60.00000000 0.00000000 0.15748031 -100.000 -60.000 0.000 0.157 +xy +a=6378137.000 +b=6356752.314245
-0 A
-5 W
-47.2441 0 T 90 R
-V N 0 0 T 354 47 scale /DeviceRGB setcolorspace
-<< /ImageType 1 /Decode [0 1 0 1 0 1] /Width 20 /Height 1 /BitsPerComponent 8
- /ImageMatrix [20 0 0 -1 0 1] /DataSource currentfile /ASCII85Decode filter /LZWDecode filter
->> image
-J1uKaoJ?u?!sV$+!uL_4gC7.M>k$ABD?c\`QsaAg`&+M(8OHY.8":$*"$FB>TEkf~>
-U
-{1 0 1 C} FS
-O1
-0 47 M
-0 -47 D
--24 24 D
-FO
-{1 0 0 C} FS
-354 47 M
-0 -47 D
-24 24 D
-FO
-0 0 M 354 0 D S
-0 47 M 354 0 D S
-0 0 M 0 47 D S
-354 0 M 0 47 D S
-1 W
-0 0 M 0 47 D S
-18 0 M 0 47 D S
-35 0 M 0 47 D S
-53 0 M 0 47 D S
-71 0 M 0 47 D S
-89 0 M 0 47 D S
-106 0 M 0 47 D S
-124 0 M 0 47 D S
-142 0 M 0 47 D S
-159 0 M 0 47 D S
-177 0 M 0 47 D S
-195 0 M 0 47 D S
-213 0 M 0 47 D S
-230 0 M 0 47 D S
-248 0 M 0 47 D S
-266 0 M 0 47 D S
-283 0 M 0 47 D S
-301 0 M 0 47 D S
-319 0 M 0 47 D S
-337 0 M 0 47 D S
-354 0 M 0 47 D S
-2 W
-0 0 M
-0 -24 D S
-0 0 M 29 F0 (-100) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-0 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-100) Z U
-18 0 M
-0 -24 D S
-0 0 M 29 F0 (-98) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-18 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-98) Z U
-35 0 M
-0 -24 D S
-0 0 M 29 F0 (-96) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-35 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-96) Z U
-53 0 M
-0 -24 D S
-0 0 M 29 F0 (-94) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-53 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-94) Z U
-71 0 M
-0 -24 D S
-0 0 M 29 F0 (-92) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-71 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-92) Z U
-89 0 M
-0 -24 D S
-0 0 M 29 F0 (-90) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-89 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-90) Z U
-106 0 M
-0 -24 D S
-0 0 M 29 F0 (-88) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-106 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-88) Z U
-124 0 M
-0 -24 D S
-0 0 M 29 F0 (-86) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-124 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-86) Z U
-142 0 M
-0 -24 D S
-0 0 M 29 F0 (-84) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-142 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-84) Z U
-159 0 M
-0 -24 D S
-0 0 M 29 F0 (-82) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-159 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-82) Z U
-177 0 M
-0 -24 D S
-0 0 M 29 F0 (-80) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-177 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-80) Z U
-195 0 M
-0 -24 D S
-0 0 M 29 F0 (-78) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-195 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-78) Z U
-213 0 M
-0 -24 D S
-0 0 M 29 F0 (-76) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-213 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-76) Z U
-230 0 M
-0 -24 D S
-0 0 M 29 F0 (-74) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-230 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-74) Z U
-248 0 M
-0 -24 D S
-0 0 M 29 F0 (-72) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-248 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-72) Z U
-266 0 M
-0 -24 D S
-0 0 M 29 F0 (-70) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-266 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-70) Z U
-283 0 M
-0 -24 D S
-0 0 M 29 F0 (-68) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-283 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-68) Z U
-301 0 M
-0 -24 D S
-0 0 M 29 F0 (-66) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-301 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-66) Z U
-319 0 M
-0 -24 D S
-0 0 M 29 F0 (-64) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-319 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-64) Z U
-337 0 M
-0 -24 D S
-0 0 M 29 F0 (-62) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-337 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-62) Z U
-354 0 M
-0 -24 D S
-0 0 M 29 F0 (-60) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-354 -110 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-60) Z U
-0 0 M 29 F0 (Sv) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-402 24 M V -90 R PSL_dim_w 2 div neg 0 G 29 F0 (Sv) Z U
--90 R -47.2441 0 T
--944.882 -59.0551 T 0 A
--- a/acoustic_backscatter.pl Mon Apr 20 14:04:04 2015 +0000
+++ b/acoustic_backscatter.pl Sun Jul 26 20:04:48 2015 +0000
@@ -1,9 +1,9 @@
#======================================================================
# A C O U S T I C _ B A C K S C A T T E R . P L
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Mon Apr 20 13:56:56 2015
+# dlm: Thu Jun 18 13:04:26 2015
# (c) 2010 A.M. Thurnherr
-# uE-Info: 24 71 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 27 82 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -22,6 +22,9 @@
# Apr 20, 2015: - added comments
# - removed SS_{min,max}_allowed_range from calc_backscatter_profs()
# - added correct_backscatter() & linterp() from laptop
+# Apr 21, 2015: - added debug statements
+# May 14, 2015: - BUG: code did not work for partial-depth casts
+# Jun 18, 2015: - removed assertion marked by ##???, which bombed on P16N1#41 DL
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -135,9 +138,6 @@
return $miy + ($x-$mix)/($max-$mix)*($may-$miy);
}
-$Sv_ref_bin = 1; # bin 2 is slightly better than bin 5 => use closest bin as reference as originally intended
- # default setting choses first bin with data; do not set to values below 1
-
sub correct_backscatter($$)
{
my($LADCP_start,$LADCP_end) = @_;
@@ -161,8 +161,18 @@
push(@refSvProf,$refSvProf[$#refSvProf]);
push(@refSvSamp,$refSvSamp[$#refSvSamp]);
}
+ for ($i=$#refSvProf-1; $i>=0; $i--) { # extrapolate upward to surface
+ next if ($refSvSamp[$i] > 0);
+ $refSvProf[$i] = $refSvProf[$i+1];
+ $refSvSamp[$i] = $refSvSamp[$i+1];
+ }
info("\tusing bin %d as reference\n",$Sv_ref_bin);
+# for ($i=0; $i<@refSvProf; $i++) {
+# print(STDERR "$i $refSvProf[$i] $refSvSamp[$i]\n");
+# }
+# die;
+
my(@dSvProf); # create profiles for all bins
for ($bin=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
my(@dSvSamp);
@@ -171,11 +181,23 @@
$dSvProf[int($depth/100)][$bin] += $sSv[$depth][$bin] / $nSv[$depth][$bin];
$dSvSamp[int($depth/100)]++;
}
+# print(STDERR "dSvProf[bin$bin] = ");
for ($i=0; $i<@dSvSamp; $i++) {
- next unless ($refSvSamp[$i] > 0) && ($dSvSamp[$i] > 0);
+ next unless ($dSvSamp[$i] > 0);
+ die("assertion failed (refSvSamp[$i] = $refSvSamp[$i])") unless ($refSvSamp[$i] > 0);
$dSvProf[$i][$bin] = $dSvProf[$i][$bin]/$dSvSamp[$i] - $refSvProf[$i];
+ $dSvProf[$i][$bin] = $dSvProf[$i-1][$bin]
+ if (abs($dSvProf[$i][$bin]) > 10);
+# printf(STDERR "%.1f ",$dSvProf[$i][$bin]);
}
- $dSvProf[$i][$bin] = $dSvProf[$i-1][$bin]; # extrapolate 100m
+ $i--; # trim deepest, often anomalous, value
+# print(STDERR "[delete] ");
+ while ($i <= int(scalar(@nSv)/100)+1) {
+ $dSvProf[$i][$bin] = $dSvProf[$i-1][$bin]; # extrapolate 100m
+# printf(STDERR "[%.1f] ",$dSvProf[$i][$bin]);
+ $i++;
+ }
+# print(STDERR "\n");
}
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) { # correct Sv data
@@ -184,9 +206,15 @@
for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {
next unless numberp($LADCP{ENSEMBLE}[$ens]->{SV}[$bin]);
my($depth) = int($bd[$bin]);
- $LADCP{ENSEMBLE}[$ens]->{SV}[$bin] -= $dSvProf[int($depth/100)][$bin];
- linterp($depth,100*int($depth/100),100*int($depth/100)+100,
- $dSvProf[int($depth/100)][$bin],$dSvProf[int($depth/100)+1][$bin]);
+ if (numberp($dSvProf[int($depth/100)][$bin]) && numberp($dSvProf[int($depth/100)+1][$bin])) {
+# print(STDERR "\n$LADCP{ENSEMBLE}[$ens]->{SV}[$bin]? $dSvProf[int($depth/100)][$bin] $dSvProf[int($depth/100)+1][$bin] int($depth/100)+1[$bin]");
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin] -= # $dSvProf[int($depth/100)][$bin];
+ linterp($depth,100*int($depth/100),100*int($depth/100)+100,
+ $dSvProf[int($depth/100)][$bin],$dSvProf[int($depth/100)+1][$bin]);
+##??? die unless ($LADCP{ENSEMBLE}[$ens]->{SV}[$bin] < 0);
+ } else {
+ $LADCP{ENSEMBLE}[$ens]->{SV}[$bin] = nan;
+ }
}
}
}
--- a/corr_scale.skel Mon Apr 20 14:04:04 2015 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,249 +0,0 @@
-PSL_font_encode 0 get 0 eq {ISOLatin1+_Encoding /Helvetica /Helvetica PSL_reencode PSL_font_encode 0 1 put} if
-0 setlinecap
-0 setlinejoin
-10 setmiterlimit
-944.882 59.0551 T
-
-%% PostScript produced by:
-%%GMT: psscale -O -K -E -D8/2/3/0.4 -Ccorr.cpt -B/:corr:
-%%PROJ: xy 70.00000000 130.00000000 0.00000000 0.15748031 70.000 130.000 0.000 0.157 +xy +a=6378137.000 +b=6356752.314245
-0 A
-5 W
-47.2441 0 T 90 R
-{0.839 0 1 C} FS
-O0
-47 30 0 0 SB
-{0.678 0 1 C} FS
-47 29 30 0 SB
-{0.518 0 1 C} FS
-47 30 59 0 SB
-{0.353 0 1 C} FS
-47 29 89 0 SB
-{0.192 0 1 C} FS
-47 30 118 0 SB
-{0.0314 0 1 C} FS
-47 29 148 0 SB
-{0 0.129 1 C} FS
-47 30 177 0 SB
-{0 0.29 1 C} FS
-47 29 207 0 SB
-{0 0.451 1 C} FS
-47 30 236 0 SB
-{0 0.612 1 C} FS
-47 29 266 0 SB
-{0 0.773 1 C} FS
-47 3 295 0 SB
-{0 0.937 1 C} FS
-47 3 298 0 SB
-{0 1 0.902 C} FS
-47 3 301 0 SB
-{0 1 0.741 C} FS
-47 3 304 0 SB
-{0 1 0.58 C} FS
-47 3 307 0 SB
-{0 1 0.42 C} FS
-47 3 310 0 SB
-{0 1 0.259 C} FS
-47 3 313 0 SB
-{0 1 0.098 C} FS
-47 3 316 0 SB
-{0.0627 1 0 C} FS
-47 3 319 0 SB
-{0.227 1 0 C} FS
-47 3 322 0 SB
-{0.388 1 0 C} FS
-47 3 325 0 SB
-{0.549 1 0 C} FS
-47 3 328 0 SB
-{0.71 1 0 C} FS
-47 3 331 0 SB
-{0.871 1 0 C} FS
-47 3 334 0 SB
-{1 0.969 0 C} FS
-47 3 337 0 SB
-{1 0.808 0 C} FS
-47 3 340 0 SB
-{1 0.647 0 C} FS
-47 2 343 0 SB
-{1 0.482 0 C} FS
-47 3 345 0 SB
-{1 0.322 0 C} FS
-47 3 348 0 SB
-{1 0.161 0 C} FS
-47 3 351 0 SB
-{1 0 1 C} FS
-O1
-0 47 M
-0 -47 D
--24 24 D
-FO
-{1 0 0 C} FS
-354 47 M
-0 -47 D
-24 24 D
-FO
-0 0 M 354 0 D S
-0 47 M 354 0 D S
-0 0 M 0 47 D S
-354 0 M 0 47 D S
-1 W
-0 0 M 0 47 D S
-30 0 M 0 47 D S
-59 0 M 0 47 D S
-89 0 M 0 47 D S
-118 0 M 0 47 D S
-148 0 M 0 47 D S
-177 0 M 0 47 D S
-207 0 M 0 47 D S
-236 0 M 0 47 D S
-266 0 M 0 47 D S
-295 0 M 0 47 D S
-298 0 M 0 47 D S
-301 0 M 0 47 D S
-304 0 M 0 47 D S
-307 0 M 0 47 D S
-310 0 M 0 47 D S
-313 0 M 0 47 D S
-316 0 M 0 47 D S
-319 0 M 0 47 D S
-322 0 M 0 47 D S
-325 0 M 0 47 D S
-328 0 M 0 47 D S
-331 0 M 0 47 D S
-334 0 M 0 47 D S
-337 0 M 0 47 D S
-340 0 M 0 47 D S
-343 0 M 0 47 D S
-345 0 M 0 47 D S
-348 0 M 0 47 D S
-351 0 M 0 47 D S
-354 0 M 0 47 D S
-2 W
-0 0 M
-0 -24 D S
-0 0 M 29 F0 (70.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-0 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (70.0) Z U
-30 0 M
-0 -24 D S
-0 0 M 29 F0 (75.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-30 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (75.0) Z U
-59 0 M
-0 -24 D S
-0 0 M 29 F0 (80.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-59 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (80.0) Z U
-89 0 M
-0 -24 D S
-0 0 M 29 F0 (85.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-89 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (85.0) Z U
-118 0 M
-0 -24 D S
-0 0 M 29 F0 (90.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-118 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (90.0) Z U
-148 0 M
-0 -24 D S
-0 0 M 29 F0 (95.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-148 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (95.0) Z U
-177 0 M
-0 -24 D S
-0 0 M 29 F0 (100.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-177 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (100.0) Z U
-207 0 M
-0 -24 D S
-0 0 M 29 F0 (105.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-207 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (105.0) Z U
-236 0 M
-0 -24 D S
-0 0 M 29 F0 (110.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-236 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (110.0) Z U
-266 0 M
-0 -24 D S
-0 0 M 29 F0 (115.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-266 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (115.0) Z U
-295 0 M
-0 -24 D S
-0 0 M 29 F0 (120.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-295 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (120.0) Z U
-298 0 M
-0 -24 D S
-0 0 M 29 F0 (120.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-298 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (120.5) Z U
-301 0 M
-0 -24 D S
-0 0 M 29 F0 (121.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-301 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (121.0) Z U
-304 0 M
-0 -24 D S
-0 0 M 29 F0 (121.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-304 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (121.5) Z U
-307 0 M
-0 -24 D S
-0 0 M 29 F0 (122.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-307 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (122.0) Z U
-310 0 M
-0 -24 D S
-0 0 M 29 F0 (122.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-310 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (122.5) Z U
-313 0 M
-0 -24 D S
-0 0 M 29 F0 (123.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-313 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (123.0) Z U
-316 0 M
-0 -24 D S
-0 0 M 29 F0 (123.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-316 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (123.5) Z U
-319 0 M
-0 -24 D S
-0 0 M 29 F0 (124.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-319 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (124.0) Z U
-322 0 M
-0 -24 D S
-0 0 M 29 F0 (124.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-322 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (124.5) Z U
-325 0 M
-0 -24 D S
-0 0 M 29 F0 (125.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-325 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (125.0) Z U
-328 0 M
-0 -24 D S
-0 0 M 29 F0 (125.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-328 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (125.5) Z U
-331 0 M
-0 -24 D S
-0 0 M 29 F0 (126.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-331 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (126.0) Z U
-334 0 M
-0 -24 D S
-0 0 M 29 F0 (126.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-334 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (126.5) Z U
-337 0 M
-0 -24 D S
-0 0 M 29 F0 (127.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-337 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (127.0) Z U
-340 0 M
-0 -24 D S
-0 0 M 29 F0 (127.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-340 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (127.5) Z U
-343 0 M
-0 -24 D S
-0 0 M 29 F0 (128.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-343 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (128.0) Z U
-345 0 M
-0 -24 D S
-0 0 M 29 F0 (128.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-345 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (128.5) Z U
-348 0 M
-0 -24 D S
-0 0 M 29 F0 (129.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-348 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (129.0) Z U
-351 0 M
-0 -24 D S
-0 0 M 29 F0 (129.5) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-351 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (129.5) Z U
-354 0 M
-0 -24 D S
-0 0 M 29 F0 (130.0) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-354 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (130.0) Z U
-0 0 M 29 F0 (corr) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-402 24 M V -90 R PSL_dim_w 2 div neg 0 G 29 F0 (corr) Z U
--90 R -47.2441 0 T
--944.882 -59.0551 T 0 A
--- a/defaults.pl Mon Apr 20 14:04:04 2015 +0000
+++ b/defaults.pl Sun Jul 26 20:04:48 2015 +0000
@@ -1,9 +1,9 @@
#======================================================================
# D E F A U L T S . P L
# doc: Tue Oct 11 17:11:21 2011
-# dlm: Mon Apr 20 13:51:17 2015
+# dlm: Sun Jul 26 17:12:51 2015
# (c) 2011 A.M. Thurnherr
-# uE-Info: 283 24 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 379 12 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -46,6 +46,15 @@
# plotting sub-system)
# - croak -> error
# - added $SS_use_BT, $SS_min_signal, $SS_min_samp
+# Apr 20, 2015: - reduced value of $SS_min_allowed_range
+# - added $Sv_ref_bin
+# Apr 21: 2015: - BUG: typo in $Sv_ref_bin
+# - decreased default verbosity
+# May 15, 2015: - added $min_valid_vels
+# May 20, 2015: - STN -> PROF
+# Jul 23, 2015: - began adaptation to libGMT.pl
+# - changed .prof output .wprof
+# - -v docu was wrong
#======================================================================
# Data Input
@@ -70,6 +79,12 @@
$pitch_bias = $roll_bias = $heading_bias = 0;
+# minimum valid velocities to require in a LADCP file
+# throws an error if not enough valid vels in a file
+
+$min_valid_vels = 50;
+
+
# bins to use in w calculations
# - set with -b
# - defaults to 2-last
@@ -81,13 +96,13 @@
#======================================================================
# - there are 4 verbosity levels, selected by -v
-# 0: only print errors
-# 1: UNIX-like (warnings and info messages that are not produced for every cast)
-# 2: (default) progress messages and useful information
+# 0 : errors
+# 1*: UNIX-like (warnings and info messages that are not produced for every cast; *DEFAULT)
+# 2 : progress messages and useful information
# >2: debug messges
# - the most useful ones of these are 1 & 2
-&antsCardOpt(\$opt_v,2);
+&antsCardOpt(\$opt_v,1);
# output bin size in meters
@@ -102,7 +117,7 @@
# output base name
-$out_basename = sprintf('%03d',$STN);
+$out_basename = sprintf('%03d',$PROF);
# output subdirectories
@@ -242,9 +257,19 @@
#======================================================================
-# Seabed Search
+# Acoustic Backscatter and Seabed Search
#======================================================================
+# After applying the method of Deines (1999), an empirical correction
+# for Sv is applied to the data. The following variable determines which
+# bin is chosen to construct a reference profile for Sv. The bin number
+# is automatically increased if the selected bin does not contain valid
+# data, i.e. the default value of 1 ensures that the closest valid bin
+# is used to construct the reference profile.
+
+$Sv_ref_bin = 1;
+
+
# Set to folloing variable to 1 to use ADCP BT data to detect seabed
# instead of default code based on Sv (echo amplitude). I do not know
# which code is better.
@@ -349,16 +374,13 @@
#----------------------------------------------------------------------
# Vertical-velocity profile output and plots:
# Data:
-# *.prof vertical velocity profiles
+# *.wprof vertical velocity profiles
# Standard Plots:
-# *_prof.eps vertical velocity profiles (main output plot)
-# Optional Plots:
-# *_wspec.eps vertical-velocity wavenumber spectra
+# *_wprof.ps vertical velocity profiles (main output plot)
#----------------------------------------------------------------------
-@out_profile = ("| LWplot_prof_2beam $plot_subdir/${out_basename}_prof.eps",
-# "| LWplot_spec $plot_subdir/${out_basename}_spec.eps",
- "$data_subdir/$out_basename.prof");
+@out_profile = ("plot_wprof($plot_subdir/${out_basename}_wprof.ps)",
+ "$data_subdir/$out_basename.wprof");
#----------------------------------------------------------------------
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/plot_wprof.pl Sun Jul 26 20:04:48 2015 +0000
@@ -0,0 +1,121 @@
+#======================================================================
+# P L O T _ W P R O F . P L
+# doc: Sun Jul 26 11:08:50 2015
+# dlm: Sun Jul 26 18:14:54 2015
+# (c) 2015 A.M. Thurnherr
+# uE-Info: 105 49 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Jul 26, 2015: - created from LWplot_prof_2beam
+
+# Tweakables:
+#
+# $plot_wprof_xmin = -0.27;
+# $plot_wprof_ymax = 5000;
+# $plot_wprof_xtics = "-0.25 -0.15 -0.05 0.05";
+
+require "$ANTS/libGMT.pl";
+
+sub setR1() { GMT_setR("-R$plot_wprof_xmin/0.35/0/$plot_wprof_ymax"); }
+sub setR2() { GMT_setR("-R-200/200/0/$plot_wprof_ymax"); }
+
+sub plotDC($)
+{
+ my($f) = @_;
+ for (my($bi)=0; $bi<=$#{$DNCAST{$f}}; $bi++) {
+ printf(GMT (numberp($DNCAST{$f}[$bi]) ? "%g %g\n" : "nan nan\n"),
+ $DNCAST{$f}[$bi],($bi+0.5)*$opt_o)
+ if ($DNCAST{N_SAMP}[$bi] >= $opt_k);
+ }
+}
+
+sub plotUC($)
+{
+ my($f) = @_;
+ for (my($bi)=0; $bi<=$#{$UPCAST{$f}}; $bi++) {
+ printf(GMT (numberp($UPCAST{$f}[$bi]) ? "%g %g\n" : "nan nan\n"),
+ $UPCAST{$f}[$bi],($bi+0.5)*$opt_o)
+ if ($UPCAST{N_SAMP}[$bi] >= $opt_k);
+ }
+}
+
+sub plotBT($)
+{
+ my($f) = @_;
+ for (my($bi)=0; $bi<=$#{$BT{$f}}; $bi++) {
+ printf(GMT (numberp($BT{$f}[$bi]) ? "%g %g\n" : "nan nan\n"),
+ $BT{$f}[$bi],($bi+0.5)*$opt_o)
+ if ($BT{N_SAMP}[$bi] >= $opt_k);
+ }
+}
+
+
+sub plot_wprof($)
+{
+ my($pfn) = @_;
+
+ $plot_wprof_xmin = -0.1
+ unless defined($plot_wprof_xmin);
+ $plot_wprof_ymax = ($P{water_depth} > 0) ?
+ round($P{water_depth} + 25) :
+ round($P{max_depth} + 25)
+ unless defined($plot_wprof_ymax);
+ $plot_wprof_xtics = "-0.05 0.05 0.15"
+ unless defined($plot_wprof_xtics);
+
+ GMT_begin($pfn,'-JX10/-10',"-R$plot_wprof_xmin/0.35/0/$plot_wprof_ymax",'-P'); # START PLOT
+
+ GMT_psxy('-G200'); # MAD background
+ print(GMT "0.07 0\n 0.07 $plot_wprof_ymax\n0.18 $plot_wprof_ymax\n0.18 0\n");
+
+ if ($P{water_depth} > 0) { # SEABED
+ GMT_psxy('-G204/153/102');
+ print(GMT "$plot_wprof_xmin $plot_wprof_ymax\n0.35 $plot_wprof_ymax\n0.35 $P{water_depth}\n $plot_wprof_xmin $P{water_depth}\n");
+ }
+
+ setR1(); # FRAME
+ GMT_psxy('-W1');
+ print(GMT "0 0\n 0 $plot_wprof_ymax\n");
+ setR2();
+ GMT_psxy('-W1 -M');
+ print(GMT ">\n50 0\n 50 $plot_wprof_ymax\n");
+ print(GMT ">\n100 0\n 100 $plot_wprof_ymax\n");
+ print(GMT ">\n150 0\n 150 $plot_wprof_ymax\n");
+
+ setR1(); # VERTICAL VELOCITIES
+ GMT_psxy('-Mn -W4,coral,6_2:0'); plotDC('MEDIAN_W12');
+ GMT_psxy('-Mn -W4,coral,4_6:0'); plotDC('MEDIAN_W34');
+ GMT_psxy('-Mn -W4,SeaGreen,6_2:0'); plotUC('MEDIAN_W12');
+ GMT_psxy('-Mn -W4,SeaGreen,4_6:0'); plotUC('MEDIAN_W34');
+ GMT_psxy('-Mn -W4,black'); plotBT('MEDIAN_W');
+
+ GMT_psxy('-Sc0.1c -Gcoral'); plotDC('MAD_W'); # MEAN ABSOLUTE DEVIATIONS
+ GMT_psxy('-Sc0.1c -GSeaGreen'); plotUC('MAD_W');
+ GMT_psxy('-Sc0.1c -Gblack'); plotBT('MAD_W');
+
+ setR2(); # SAMPLES
+ GMT_psxy('-Mn -W1/coral'); plotDC('N_SAMP');
+ GMT_psxy('-Mn -W1/SeaGreen'); plotUC('N_SAMP');
+ GMT_psxy('-Mn -W1/black'); plotBT('N_SAMP');
+
+ GMT_unitcoords(); # LABELS
+ GMT_pstext();
+ print(GMT "0.02 0.02 12 0 0 TL $P{out_basename} $P{run_label}\n");
+ print(GMT "0.6 0.98 12 0 0 BR m.a.d.\n");
+
+ my($depth_tics) = ($plot_wprof_ymax < 1000 ) ? 'f10a100' : 'f100a500'; # AXES
+ setR1();
+ GMT_psbasemap("-Bf0.01:'Vertical Velocity [m/s] ':/$depth_tics:'Depth [m]':WeS");
+ foreach my $t (split('\s+',$plot_wprof_xtics)) {
+ GMT_psbasemap(sprintf('-Ba10-%fS',10-$t));
+ }
+ setR2();
+ GMT_psbasemap('-Bf10a1000-950:" # of Samples":N');
+ GMT_psbasemap('-Ba1000-900N');
+ GMT_psbasemap('-Ba1000-850N');
+
+ GMT_end(); # FINISH PLOT
+}
+
+1; # return true on require
--- a/residuals_scale.skel Mon Apr 20 14:04:04 2015 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,88 +0,0 @@
-PSL_font_encode 0 get 0 eq {ISOLatin1+_Encoding /Helvetica /Helvetica PSL_reencode PSL_font_encode 0 1 put} if
-0 setlinecap
-0 setlinejoin
-10 setmiterlimit
-944.882 59.0551 T
-
-%% PostScript produced by:
-%%GMT: psscale -O -K -E -D8/2/3/0.4 -Cresiduals.cpt -B/:residuals:
-%%PROJ: xy -0.03000000 0.03000000 0.00000000 0.15748031 -0.030 0.030 0.000 0.157 +xy +a=6378137.000 +b=6356752.314245
-0 A
-5 W
-47.2441 0 T 90 R
-{0.286 0.286 1 C} FS
-O0
-47 59 0 0 SB
-{0.573 0.573 1 C} FS
-47 59 59 0 SB
-{0.784 0.784 1 C} FS
-47 30 118 0 SB
-{1 A} FS
-47 59 148 0 SB
-{1 0.784 0.784 C} FS
-47 29 207 0 SB
-{1 0.573 0.573 C} FS
-47 59 236 0 SB
-{1 0.286 0.286 C} FS
-47 59 295 0 SB
-{0 0 1 C} FS
-O1
-0 47 M
-0 -47 D
--24 24 D
-FO
-{1 0 0 C} FS
-354 47 M
-0 -47 D
-24 24 D
-FO
-0 0 M 354 0 D S
-0 47 M 354 0 D S
-0 0 M 0 47 D S
-354 0 M 0 47 D S
-1 W
-0 0 M 0 47 D S
-59 0 M 0 47 D S
-118 0 M 0 47 D S
-148 0 M 0 47 D S
-207 0 M 0 47 D S
-236 0 M 0 47 D S
-295 0 M 0 47 D S
-354 0 M 0 47 D S
-2 W
-0 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.030) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-0 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.030) Z U
-59 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.020) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-59 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.020) Z U
-118 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.010) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-118 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.010) Z U
-148 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.005) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-148 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.005) Z U
-207 0 M
-0 -24 D S
-0 0 M 29 F0 (0.005) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-207 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.005) Z U
-236 0 M
-0 -24 D S
-0 0 M 29 F0 (0.010) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-236 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.010) Z U
-295 0 M
-0 -24 D S
-0 0 M 29 F0 (0.020) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-295 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.020) Z U
-354 0 M
-0 -24 D S
-0 0 M 29 F0 (0.030) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-354 -135 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.030) Z U
-0 0 M 29 F0 (residuals) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-402 24 M V -90 R PSL_dim_w 2 div neg 0 G 29 F0 (residuals) Z U
--90 R -47.2441 0 T
--944.882 -59.0551 T 0 A
--- a/time_lag.pl Mon Apr 20 14:04:04 2015 +0000
+++ b/time_lag.pl Sun Jul 26 20:04:48 2015 +0000
@@ -1,9 +1,9 @@
#======================================================================
# T I M E _ L A G . P L
# doc: Fri Dec 17 21:59:07 2010
-# dlm: Thu Apr 16 12:13:25 2015
+# dlm: Fri Jun 19 07:23:38 2015
# (c) 2010 A.M. Thurnherr
-# uE-Info: 276 41 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 61 10 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -57,6 +57,8 @@
# - BUG: executable flag was not set on file output
# - disabled active output when ANTS are not available
# - croak -> error
+# May 15, 2015: - fiddled with assertions
+# Jun 19, 2015: - disabled L2 warning on partial-depth time-lagging failures
# DIFFICULT STATIONS:
# NBP0901#131 this requires the search-radius doubling heuristic
@@ -106,9 +108,6 @@
sub bestLag($$$$) # find best lag in window
{
my($fe,$le,$ww,$soi) = @_; # first/last LADCP ens, window width, scan-offset increment
- die("assertion failed\n\tfe = $fe, le = $le, firstGoodEns = $firstGoodEns, lastGoodEns = $lastGoodEns")
- unless ($fe>=$firstGoodEns && $le<=$lastGoodEns);
-
my($bestso) = 0; # error at first-guess offset
my($bestmad) = mad_w($fe,$le,0);
@@ -159,13 +158,15 @@
my($n_valid_windows) = 0;
$first_ens = $approx_joint_profile_start_ens
- if ($first_ens < $approx_joint_profile_start_ens);
+ if ($first_ens < $approx_joint_profile_start_ens);
my($last_lag_piece) = ($last_ens == $lastGoodEns); # none is following
$last_ens = $approx_joint_profile_end_ens
if ($last_ens > $approx_joint_profile_end_ens);
for (my($wi)=0; $wi<$n_windows; $wi++) {
my($fe) = $first_ens + int(($last_ens-$first_ens-$window_ens)*$wi/($n_windows-1)+0.5);
+ die("assertion failed\n\tfe = $fe, first_ens = $first_ens, last_ens = $last_ens, window_ens = $window_ens, firstGoodEns = $firstGoodEns, lastGoodEns = $lastGoodEns")
+ unless ($fe>=$firstGoodEns && $fe+$window_ens<=$lastGoodEns);
my($so,$mad) = bestLag($fe,$fe+$window_ens,$search_radius,$scan_increment);
$elapsed[$wi] = $LADCP{ENSEMBLE}[$fe+int($w_size/2/$LADCP{MEAN_DT}+0.5)]->{ELAPSED};
die("assertion failed\nfe=$fe, lastGoodEns=$lastGoodEns, w_size=$w_size") unless ($elapsed[$wi]);
@@ -215,7 +216,8 @@
unless ($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]} >= $opt_3*$n_valid_windows) {
if (max(@best_lag)-min(@best_lag) > $TL_max_allowed_three_lag_spread) {
warning(2,"$0: cannot determine a valid $ctmsg lag; top 3 tags account for %d%% of total (use -3 to relax criterion)\n",
- int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_valid_windows+0.5));
+ int(100*($nBest{$best_lag[0]}+$nBest{$best_lag[1]}+$nBest{$best_lag[2]})/$n_valid_windows+0.5))
+ unless ($ctmsg == 'partial-cast');
$failed = 1;
} else {
warning(1,"top 3 tags account for only %d%% of total\n",
--- a/w_scale.skel Mon Apr 20 14:04:04 2015 +0000
+++ /dev/null Thu Jan 01 00:00:00 1970 +0000
@@ -1,116 +0,0 @@
-PSL_font_encode 0 get 0 eq {ISOLatin1+_Encoding /Helvetica /Helvetica PSL_reencode PSL_font_encode 0 1 put} if
-0 setlinecap
-0 setlinejoin
-10 setmiterlimit
-944.882 59.0551 T
-
-%% PostScript produced by:
-%%GMT: psscale -O -K -E -D8/2/3/0.4 -Cw.cpt -B/:w:
-%%PROJ: xy -0.07000000 0.07000000 0.00000000 0.15748031 -0.070 0.070 0.000 0.157 +xy +a=6378137.000 +b=6356752.314245
-0 A
-5 W
-47.2441 0 T 90 R
-{0.251 0.251 1 C} FS
-O0
-47 51 0 0 SB
-{0.439 0.439 1 C} FS
-47 25 51 0 SB
-{0.561 0.561 1 C} FS
-47 25 76 0 SB
-{0.686 0.686 1 C} FS
-47 26 101 0 SB
-{0.812 0.812 1 C} FS
-47 25 127 0 SB
-{1 A} FS
-47 50 152 0 SB
-{1 0.812 0.812 C} FS
-47 26 202 0 SB
-{1 0.686 0.686 C} FS
-47 25 228 0 SB
-{1 0.561 0.561 C} FS
-47 25 253 0 SB
-{1 0.439 0.439 C} FS
-47 26 278 0 SB
-{1 0.251 0.251 C} FS
-47 50 304 0 SB
-{0 0 1 C} FS
-O1
-0 47 M
-0 -47 D
--24 24 D
-FO
-{1 0 0 C} FS
-354 47 M
-0 -47 D
-24 24 D
-FO
-0 0 M 354 0 D S
-0 47 M 354 0 D S
-0 0 M 0 47 D S
-354 0 M 0 47 D S
-1 W
-0 0 M 0 47 D S
-51 0 M 0 47 D S
-76 0 M 0 47 D S
-101 0 M 0 47 D S
-127 0 M 0 47 D S
-152 0 M 0 47 D S
-202 0 M 0 47 D S
-228 0 M 0 47 D S
-253 0 M 0 47 D S
-278 0 M 0 47 D S
-304 0 M 0 47 D S
-354 0 M 0 47 D S
-2 W
-0 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.07) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-0 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.07) Z U
-51 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.05) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-51 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.05) Z U
-76 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.04) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-76 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.04) Z U
-101 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.03) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-101 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.03) Z U
-127 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.02) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-127 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.02) Z U
-152 0 M
-0 -24 D S
-0 0 M 29 F0 (-0.01) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-152 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (-0.01) Z U
-202 0 M
-0 -24 D S
-0 0 M 29 F0 (0.01) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-202 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.01) Z U
-228 0 M
-0 -24 D S
-0 0 M 29 F0 (0.02) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-228 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.02) Z U
-253 0 M
-0 -24 D S
-0 0 M 29 F0 (0.03) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-253 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.03) Z U
-278 0 M
-0 -24 D S
-0 0 M 29 F0 (0.04) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-278 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.04) Z U
-304 0 M
-0 -24 D S
-0 0 M 29 F0 (0.05) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-304 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.05) Z U
-354 0 M
-0 -24 D S
-0 0 M 29 F0 (0.07) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-354 -119 M V -90 R PSL_dim_w neg PSL_dim_h 2 div neg G 29 F0 (0.07) Z U
-0 0 M 29 F0 (w) E /PSL_dim_w exch def FP pathbbox N /PSL_dim_h exch def pop /PSL_dim_d exch def pop
-402 24 M V -90 R PSL_dim_w 2 div neg 0 G 29 F0 (w) Z U
--90 R -47.2441 0 T
--944.882 -59.0551 T 0 A