#!/usr/bin/perl
#======================================================================
# L A D C P _ W _ P O S T P R O C
# doc: Fri Apr 24 17:15:59 2015
# dlm: Wed Mar 9 17:20:46 2016
# (c) 2015 A.M. Thurnherr
# uE-Info: 510 1 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'edit and re-grid LADCP vertical-velocity samples';
# 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 %outgrid_*
# Jul 27, 2015: - B-t)rack <wprof>
# Sep 21, 2015: - added valid_bins()
# - replaced -t by -b
# - BUG: -k default was wrong
# - added error message output to EditParams
# Sep 22, 2015: - added output_resolution()
# Sep 24, 2015: - BUG: -k did not work any more
# - BUG: some plot commands were exectuted even without -p
# - BUG: wprof plot did not respect -k
# - BUG: typo in mad output
# Sep 26, 2015: - allow $ID as alias for $STN and $PROF
# Sep 27, 2015: - adapted to 'unknown' %water_depth
# Oct 12, 2015: - improved EditParams handling
# - BUG: single-profile plot did not work any more
# - BUG: wrong single-profile bin size warning produced
# - BUG: -o /-k code was wrong
# - BUG: plot was in landscape mode
# - added run-label(s) to plot
# - require ANTSlibs V6.2 for release
# Oct 13, 2015: - adapted to [version.pl]
# Jan 20, 2016: - added dc_w.diff, uc_w.diff for dual-head profiles
# Jan 21, 2016: - added %_wcorr.ps plot output
# - added correlation-based QC to wprof plot
# Jan 22, 2016: - added output -d)ir option
# Jan 24, 2016: - BUG: uc_/dc_corr returned single nan on failure
# Jan 25, 2016: - added software version %PARAM
# Jan 26, 2016: - added -v
# Jan 30, 2016: - added -w
# Feb 14, 2016: - BUG: all bins were off by 1! (-v, inherited limits)
# Mar 1, 2016: - added required ADCP sampling %PARAMs to dual-head
# output
# Mar 7, 2016: - BUG: correlation stats were defined/used for single-head data
# - removed good_bins() from library as -v allows more control
($ANTS) = (`which ANTSlib` =~ m{^(.*)/[^/]*$});
($WCALC) = ($0 =~ m{^(.*)/[^/]*$});
$WCALC = '.' if ($WCALC eq '');
die("$0: ANTSlib required but not found (bad \$PATH?)\n")
unless ($ANTS ne '');
require "$WCALC/version.pl";
require "$ANTS/ants.pl";
require "$ANTS/libstats.pl";
require "$ANTS/libGMT.pl";
&antsAddParams('LADCP_w_postproc',"Version $VERSION");
&antsUsage('b:d:i:k:l:o:p:v:w:',1,
'[profile -i)d <id>]',
'[-o)utput bin <resolution>] [-k) require <min> samples]',
'[-v)alid bins <DL first>,<DL last>[,<UL first>,<UL_last>]',
'[-w) <DL_dc_field>,<DL_uc_field>[,<UL_dc_field>,<UL_uc_field>]',
'[-s)urface-layer <limit[300m]>]',
'[ouptput -d)ir <name>]',
'[-p)lot <[%03d_wprof.eps]> [-b)t <wprof>]]',
'<DL.wsamp file> [UL.wsamp file] (or only <UL.wsamp file>)');
$dual_head = (@ARGV==1); # single or dual head
$id = defined($opt_i) ? $opt_i : &antsParam('profile_id'); # ensure profile id
croak("$0: no profile_id in first file => -i required\n")
unless defined($id);
if (defined($opt_d)) { # different output directory
croak("$opt_d: no such directory\n")
unless (-d "$opt_d");
}
&antsCardOpt(\$opt_s,300); # surface layer depth limit
if (defined($opt_v)) { # bin ranges with valid data to use
($fvBin,$lvBin,$UL_fvBin,$UL_lvBin) = split(/,/,$opt_v);
croak("$0: cannot decode -v $opt_v\n")
unless (defined($lvBin) && (!$dual_head || defined($UL_lvBin)));
$fvBin = &antsRequireParam('LADCP_firstBin') if ($fvBin eq '*'); # corresponding UL values set below
$lvBin = &antsRequireParam('LADCP_lastBin') if ($lvBin eq '*');
&antsAddParams('DL_first_valid_bin',$fvBin,
'DL_last_valid_bin', $lvBin);
} else {
&antsAddParams('DL_first_valid_bin',&antsRequireParam('outgrid_firstbin'),
'DL_last_valid_bin',&antsRequireParam('outgrid_lastbin'));
}
if (defined($opt_w)) { # vertical-velocity fields to use
($Ddwf,$Duwf,$Udwf,$Uuwf) = split(/,/,$opt_w); # DL dc, DL uc, ...
croak("$0: cannot decode -w $opt_w\n")
unless (defined($Duwf) && (!$dual_head || defined($Uuwf)));
&antsAddParams('DL_dc_w_field',$Ddwf,'DL_uc_w_field',$Duwf,
'UL_dc_w_field',$Udwf,'UL_uc_w_field',$Uuwf);
} else {
($Ddwf,$Duwf,$Udwf,$Uuwf) = ('w','w','w','w');
}
if (defined($opt_o)) { # output grid resolution
$opt_o_override = 1;
&antsCardOpt($opt_o);
} else {
$opt_o = &antsRequireParam('outgrid_dz');
}
if (defined($opt_k)) { # minimum number of required samples
$opt_k_override = 1;
&antsCardOpt($opt_k);
} else {
$opt_k = &antsRequireParam('outgrid_minsamp');
$opt_k *= 2 if ($dual_head);
}
#----------------------------------------------------------------------
# Redirect STDOUT to %.wprof & create plots if STDOUT is a tty
#----------------------------------------------------------------------
if (-t STDOUT) {
$opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps,$opt_d/%03d_wcorr.ps"
: '%03d_wprof.ps,%03d_wcorr.ps'
unless defined($opt_p);
$outfile = defined($opt_d) ? sprintf('%s/%03d.wprof',$opt_d,$id)
: sprintf('%03d.wprof',$id);
open(STDOUT,">$outfile") || die("$outfile: $!\n");
}
croak("$0: -b only makes sense when plots are produced\n")
if $opt_b && !defined($opt_p);
#----------------------------------------------------------------------
# EditParams Library
#
# output_resolution(dz) output_resolution(40)
# bad_range[_dc|_uc](field,min_val,max_val) bad_range_uc('depth',3500,3600)
#
#----------------------------------------------------------------------
my(@brFnr,@brMin,@brMax,@brDUc);
my(@gbFirst,@gbLast);
sub output_resolution($)
{
my($dz) = @_;
$opt_o = $dz;
$opt_o_override = 1;
}
sub bad_range($$$)
{
push(@brFnr,&fnr(shift));
push(@brMin,shift); push(@brMax,shift);
$brMin[$#brMin] = -9e99 if ($brMin[$#brMin] eq '*');
$brMax[$#brMax] = 9e99 if ($brMax[$#brMax] eq '*');
push(@brDUc,2);
}
sub bad_range_dc($$$)
{
push(@brFnr,&fnr(shift));
push(@brMin,shift); push(@brMax,shift);
$brMin[$#brMin] = -9e99 if ($brMin[$#brMin] eq '*');
$brMax[$#brMax] = 9e99 if ($brMax[$#brMax] eq '*');
push(@brDUc,1);
}
sub bad_range_uc($$$)
{
push(@brFnr,&fnr(shift));
push(@brMin,shift); push(@brMax,shift);
$brMin[$#brMin] = -9e99 if ($brMin[$#brMin] eq '*');
$brMax[$#brMax] = 9e99 if ($brMax[$#brMax] eq '*');
push(@brDUc,0);
}
#----------------------------------------------------------------------
sub isBad()
{
for (my($f)=0; $f<$antsBufNFields; $f++) {
for (my($i)=0; $i<@gbFirst; $i++) { # good bin range
return 1 if ($ants_[0][$bF] < $gbFirst[$i] ||
$ants_[0][$bF] > $gbLast[$i]);
}
for (my($i)=0; $i<@brFnr; $i++) { # bad ranges
next unless ($brFnr[$i] == $f);
next unless ($brDUc[$i]==2 || $brDUc[$i]==$ants_[0][$dcF]);
return 1 if ($ants_[0][$f] >= $brMin[$i] &&
$ants_[0][$f] <= $brMax[$i]);
}
}
return 0;
}
#----------------------------------------------------------------------
# Main Program
#----------------------------------------------------------------------
sub dc_corr()
{
my($n) = 0;
my($ax,$ay) = (0,0);
for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
$n++;
$ax += $DL_dc_median[$bi];
$ay += $UL_dc_median[$bi];
}
return (nan,nan,nan) unless ($n > 2);
$ax /= $n;
$ay /= $n;
my($syy,$sxy,$sxx) = (0,0,0);
for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
my($xt) = $DL_dc_median[$bi] - $ax;
my($yt) = $UL_dc_median[$bi] - $ay;
$sxx += $xt * $xt;
$syy += $yt * $yt;
$sxy += $xt * $yt;
}
my($R) = $sxy/(sqrt($sxx * $syy) + 1e-16);
my($esig) = sqrt(($R**2)*$sxx/$n);
my($rsig) = sqrt((1-$R**2)*$sxx/$n);
return ($R,$esig,$rsig);
}
sub uc_corr(@)
{
my($n) = 0;
my($ax,$ay) = (0,0);
for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
$n++;
$ax += $DL_uc_median[$bi];
$ay += $UL_uc_median[$bi];
}
return (nan,nan,nan) unless ($n > 2);
$ax /= $n;
$ay /= $n;
my($syy,$sxy,$sxx) = (0,0,0);
for (my($bi)=int($opt_s/$opt_o); $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
my($xt) = $DL_uc_median[$bi] - $ax;
my($yt) = $UL_uc_median[$bi] - $ay;
$sxx += $xt * $xt;
$syy += $yt * $yt;
$sxy += $xt * $yt;
}
my($R) = $sxy/(sqrt($sxx * $syy) + 1e-16);
my($esig) = sqrt(($R**2)*$sxx/$n);
my($rsig) = sqrt((1-$R**2)*$sxx/$n);
return ($R,$esig,$rsig);
}
#----------------------------------------------------------------------
$dcwF = &fnr($Ddwf); $ucwF = &fnr($Duwf);
$eF = &fnr('elapsed');
$dF = &fnr('depth');
$bF = &fnr('bin');
$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);
if (defined($opt_p)) {
($sumPF,$corrPF) = split(/,/,$opt_p);
croak("$0: cannot decode -p $opt_p\n")
unless (length($corrPF)>0);
}
my($R,$R2);
if (defined($opt_p)) { # begin summary plot
$xmin = -0.1; $x2min = -200;
$xmax = 0.35; $x2max = 200;
$ymin = 0;
$ymax = antsParam('water_depth');
$ymax = antsRequireParam('max_depth') unless numberp($ymax);
$plotsize = 13;
$R = "-R$xmin/$xmax/$ymin/$ymax";
$R2 = "-R$x2min/$x2max/$ymin/$ymax";
GMT_begin(sprintf($sumPF,$id),"-JX$plotsize/-$plotsize",$R,'-P -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);
if (defined($opt_b)) {
open(BT,$opt_b) || croak("$opt_b: $!\n");
@BTL = &antsFileLayout(BT);
my($BTwf,$BTdf,$BTmf,$BTnf);
for (my($f)=0; $f<=$#BTL; $f++) {
$BTdf = $f if ($BTL[$f] eq 'depth');
$BTwf = $f if ($BTL[$f] eq 'BT_w');
$BTmf = $f if ($BTL[$f] eq 'BT_w.mad');
$BTnf = $f if ($BTL[$f] eq 'BT_w.nsamp');
}
croak("$opt_b: file-layout error\n")
unless defined($BTdf) && defined($BTwf) &&
defined($BTmf) && defined($BTnf);
GMT_psxy('-W6');
while (@BT = &antsFileIn(BT)) {
next unless numberp($BT[$BTwf]);
printf(GMT "%f %f\n",$BT[$BTwf],$BT[$BTdf]);
}
}
}
$min_depth = 9e99; # sentinels
$max_depth = -9e99;
$curF = $P{PATHNAME}; # current input file (sentinel)
$filt = 0;
for ($r=0; &antsIn(); $r++) {
if ($P{PATHNAME} ne $curF) { # 2nd file (UL data)
$curF = $P{PATHNAME};
$dcwF = &fnr($Udwf); $ucwF = &fnr($Uuwf);
&antsInfo("WARNING: inconsistent %%outgrid_dz in input files") # consistency checks
if (defined($P{outgrid_dz}) && $P{outgrid_dz}!=$opt_o &&!$opt_o_override);
&antsInfo("WARNING: inconsistent %%outgrid_minsamp in input files")
if defined($P{outgrid_minsamp}) && !$opt_k_override &&
(( $dual_head && $P{outgrid_minsamp}*2!=$opt_k) ||
(!$dual_head && $P{outgrid_minsamp}!=$opt_k));
if (defined($opt_v)) {
$fvBin = &antsRequireParam('LADCP_firstBin') if ($UL_fvBin eq '*'); # valid bin ranges
$lvBin = &antsRequireParam('LADCP_lastBin') if ($UL_lvBin eq '*');
&antsAddParams('UL_first_valid_bin',$fvBin,
'UL_last_valid_bin', $lvBin);
} else {
&antsAddParams('UL_first_valid_bin',&antsRequireParam('outgrid_firstbin'),
'UL_last_valid_bin',&antsRequireParam('outgrid_lastbin'));
}
for (my($bi)=0; $bi<=$#dcw1; $bi++) { # calc DL median profile (before reading UL data)
$DL_dc_median[$bi] = median(@{$dcw1[$bi]});
$DL_uc_median[$bi] = median(@{$ucw1[$bi]});
}
#
# ADCP bin length, pulse length, and blanking distance for dual head casts
# with inconsistent values:
# bin length: use smaller value, which will lead to smaller spectral correction
# pulse length: same
# blanking distance: use smaller value, which is conservative e.g. for filters for ringing
#
my($warned);
unless (round($bin_length) == round($P{ADCP_bin_length})) {
unless ($warned) {
&antsInfo("WARNING: inconsistent ADCP sampling parameters --- using conservative values");
$warned = 1;
}
$bin_length = min($bin_length,$P{ADCP_bin_length});
}
unless (round($pulse_length) == round($P{ADCP_pulse_length})) {
unless ($warned) {
&antsInfo("WARNING: inconsistent ADCP sampling parameters --- using conservative values");
$warned = 1;
}
$pulse_length = min($pulse_length,$P{ADCP_pulse_length});
}
unless (round($blanking_dist) == round($P{ADCP_blanking_distance})) {
unless ($warned) {
&antsInfo("WARNING: inconsistent ADCP sampling parameters --- using conservative values");
$warned = 1;
}
$blanking_dist = min($blanking_dist,$P{ADCP_blanking_distance});
}
$PROF = $STN = $ID = $id; $RUN = antsRequireParam('run_label'); # set variables for editing
undef(@rngMin); undef(@rngMax); undef(@bins);
unless ($return = do "./EditParams") { # man perlfunc
croak("./EditParams: $@\n") if ($@);
}
if (defined($opt_p)) { # 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",$DL_dc_median[$bi],($bi+0.5)*$opt_o);
}
GMT_psxy('-W2,46/139/87,-');
for (my($bi)=0; $bi<=$#ucw1; $bi++) {
printf(GMT "%f %f\n",$DL_uc_median[$bi],($bi+0.5)*$opt_o);
}
undef(@dcw1); undef(@ucw1);
}
} # of 2nd file started
if (defined($opt_v)) { # explicit ranges of validity given
next if ($ants_[0][$bF]<$fvBin || # => apply them
$ants_[0][$bF]>$lvBin);
} else { # no range of valid bins given
next if ($ants_[0][$bF]<$P{outgrid_firstbin} || # => use values from [LADCP_w_ocean]
$ants_[0][$bF]>$P{outgrid_lastbin});
}
$filt++,next if &isBad(); # additional editing
$min_depth = $ants_[0][$dF] if ($ants_[0][$dF] < $min_depth); # update depth limits
$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][$dcwF]); # vertical velocity
push(@{$dcw1[$bi]},$ants_[0][$dcwF]) if ($dual_head); # single-instrument w
push(@{$dce[$bi]},$ants_[0][$eF]); # elapsed time
} else { # upcast
push(@{$ucw[$bi]},$ants_[0][$ucwF]);
push(@{$ucw1[$bi]},$ants_[0][$ucwF]) if ($dual_head);
push(@{$uce[$bi]},$ants_[0][$eF]);
}
} # file-read loop
if ($dual_head) {
for (my($bi)=0; $bi<=$#dcw1; $bi++) { # calc UL median & difference profiles
$UL_dc_median[$bi] = median(@{$dcw1[$bi]});
$UL_uc_median[$bi] = median(@{$ucw1[$bi]});
$dc_diff[$bi] = numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi])
? $DL_dc_median[$bi] - $UL_dc_median[$bi] : nan;
$uc_diff[$bi] = numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi])
? $DL_uc_median[$bi] - $UL_uc_median[$bi] : nan;
}
($dc_R,$dc_esig,$dc_rsig) = &dc_corr(); # correlation statistics
($uc_R,$uc_esig,$uc_rsig) = &uc_corr();
&antsAddParams('dc_R',$dc_R,'uc_R',$uc_R,
'dc_explained_stddev',$dc_esig,'uc_explained_stddev',$uc_esig,
'dc_residual_stddev',$dc_rsig,'uc_residual_stddev',$uc_rsig);
if (defined($opt_p)) { # plot 2nd-instrument profiles
GMT_psxy('-W2,255/127/80,.');
for (my($bi)=0; $bi<=$#dcw1; $bi++) {
printf(GMT "%f %f\n",$UL_dc_median[$bi],($bi+0.5)*$opt_o);
}
GMT_psxy('-W2,46/139/87,.');
for (my($bi)=0; $bi<=$#ucw1; $bi++) {
printf(GMT "%f %f\n",$UL_uc_median[$bi],($bi+0.5)*$opt_o);
}
}
}
&antsInfo("%d measurements edited (%d%% of total)",$filt,round(100*$filt/$r))
if ($filt > 0);
#----------------------------------------------------------------------
# Average and Output Profiles, Add to Summary Plot
#----------------------------------------------------------------------
if ($dual_head) { # dual-head output
@antsNewLayout = ('depth','hab',
'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp',
'dc_w.diff','uc_w.diff'); # DL-UL differences
&antsAddParams('profile_id',$id,'lat',$P{lat},'lon',$P{lon}); # selected %PARAMs
&antsAddParams($dayNoP,$dn,'run_label',"$first_label & $P{run_label}");
&antsAddParams('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k);
&antsAddParams('min_depth',round($min_depth),'max_depth',round($max_depth));
&antsAddParams('water_depth',$P{water_depth},'water_depth.sig',$P{water_depth.sig});
&antsAddParams('ADCP_bin_length',$bin_length,'ADCP_pulse_length',$pulse_length);
&antsAddParams('ADCP_blanking_distance',$blanking_dist);
undef($antsOldHeaders);
} else {
@antsNewLayout = ('depth','hab', # single-head output
'dc_elapsed','dc_w','dc_w.mad','dc_w.nsamp',
'uc_elapsed','uc_w','uc_w.mad','uc_w.nsamp');
}
#&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)?$dcwmad[$bi]:nan),
scalar(@{$dcw[$bi]}),
avg(@{$uce[$bi]}),
(($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan),(($ucns[$bi]>=$opt_k)?$ucwmad[$bi]:nan),
scalar(@{$ucw[$bi]}));
push(@{$out[$bi]},$dc_diff[$bi],$uc_diff[$bi])
if ($dual_head);
&antsOut(@{$out[$bi]});
}
if (defined($opt_p)) { # complete summary plot
if ($dual_head) {
GMT_psxy('-W2/100/100/255'); # surface layer limit
print(GMT "-0.1 $opt_s\n0.07 $opt_s\n");
if ($dc_R < 0.3 || !numberp($dc_R)) { # correlation statistics
&antsInfo("WARNING: low dc correlation (r = %.1f) between UL and DL data",$dc_R);
GMT_pstext('-Gwhite -Wred');
} elsif ($dc_R < 0.5) { GMT_pstext('-Gblack -Wyellow'); }
else { GMT_pstext('-Gblack -Wgreen'); }
printf(GMT "%f %f 12 0 0 BL %.1f\n",-0.07,0.9*$ymax,$dc_R);
}
GMT_pstext('-G255/127/80');
printf(GMT "%f %f 12 0 0 BL dc\n",-0.095,0.9*$ymax);
printf(GMT "%f %f 12 0 0 BL [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
0.02,0.9*$ymax,100*$dc_esig,100*$dc_rsig) if ($dual_head);
if ($dual_head) {
if ($uc_R < 0.3 || !numberp($uc_R)) {
&antsInfo("WARNING: low uc correlation (r = %.1f) between UL and DL data",$uc_R);
GMT_pstext('-Gwhite -Wred');
} elsif ($uc_R < 0.5) { GMT_pstext('-Gblack -Wyellow'); }
else { GMT_pstext('-Gblack -Wgreen'); }
printf(GMT "%f %f 12 0 0 BL %.1f\n",-0.07,0.95*$ymax,$uc_R);
}
GMT_pstext('-G46/139/87');
printf(GMT "%f %f 12 0 0 BL uc\n",-0.095,0.95*$ymax);
printf(GMT "%f %f 12 0 0 BL [%.1f/%.1f cm/s @~s@~\@-e/r\@-]\n",
0.02,0.95*$ymax,100*$uc_esig,100*$uc_rsig) if ($dual_head);
GMT_setR($R);
GMT_psxy('-W4,255/127/80'); # median profiles
for (my($bi)=0; $bi<=$#dcw; $bi++) {
printf(GMT "%f %f\n",(($dcns[$bi]>=$opt_k)?$dcwm[$bi]:nan),($bi+0.5)*$opt_o);
}
GMT_psxy('-W4,46/139/87');
for (my($bi)=0; $bi<=$#ucw; $bi++) {
printf(GMT "%f %f\n",(($ucns[$bi]>=$opt_k)?$ucwm[$bi]:nan),($bi+0.5)*$opt_o);
}
GMT_psxy('-Sc0.1 -G255/127/80'); # m.a.d. profiles
for (my($bi)=0; $bi<=$#dcw; $bi++) {
printf(GMT "%f %f\n",(($dcns[$bi]>=$opt_k)?$dcwmad[$bi]:nan),($bi+0.5)*$opt_o);
}
GMT_psxy('-Sc0.1 -G46/139/87');
for (my($bi)=0; $bi<=$#ucw; $bi++) {
printf(GMT "%f %f\n",(($ucns[$bi]>=$opt_k)?$ucwmad[$bi]:nan),($bi+0.5)*$opt_o);
}
GMT_setR($R2);
GMT_psxy('-Mn -W1,255/127/80');
for (my($bi)=0; $bi<=$#dcw; $bi++) { # number of samples
printf(GMT "%f %f\n",$dcns[$bi],($bi+0.5)*$opt_o);
}
GMT_psxy('-Mn -W1,46/139/87');
for (my($bi)=0; $bi<=$#dcw; $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 [$P{run_label}]\n"); }
else { printf(GMT "0.01 -0.06 14 0 0 TL %03d\n [$P{run_label}]",$id); }
GMT_pstext();
print(GMT '0.62 0.98 12 0 0 MR m.a.d.');
GMT_end();
if ($dual_head) { # correlation plot
my($mwm) = 0.05; # max(|w|) for axes
for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
$mwm = abs($DL_dc_median[$bi]) if abs($DL_dc_median[$bi]) > $mwm;
$mwm = abs($UL_dc_median[$bi]) if abs($UL_dc_median[$bi]) > $mwm;
}
for (my($bi)=0; $bi<@DL_uc_median; $bi++) {
next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
$mwm = abs($DL_uc_median[$bi]) if abs($DL_uc_median[$bi]) > $mwm;
$mwm = abs($UL_uc_median[$bi]) if abs($UL_uc_median[$bi]) > $mwm;
}
$mwm = int(100*$mwm+0.9999) / 100;
$R = "-R-$mwm/$mwm/-$mwm/$mwm";
GMT_begin(sprintf($corrPF,$id),"-JX$plotsize/$plotsize",$R,'-P -X6 -Y4');
GMT_psxy('-Ggrey80 -L');
printf(GMT "%g %g\n%g %g\n%g %g\n%g %g\n%g %g\n%g %g\n",
-$mwm, -$mwm+0.01,
-$mwm, -$mwm,
-$mwm+0.01, -$mwm,
$mwm, $mwm-0.01,
$mwm, $mwm,
$mwm-0.01, $mwm);
GMT_psxy('-W4,grey50');
print(GMT "-$mwm -$mwm\n$mwm $mwm\n");
GMT_psxy('-Sc0.12c -G255/127/80 -W1,blue,-');
for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
my($depth) = ($bi+0.5)*$opt_o;
last if ($depth > $opt_s);
print(GMT "$DL_dc_median[$bi] $UL_dc_median[$bi]\n");
}
GMT_psxy('-Sc0.12c -G255/127/80');
for (my($bi)=0; $bi<@DL_dc_median; $bi++) {
next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
my($depth) = ($bi+0.5)*$opt_o;
next unless ($depth > $opt_s);
print(GMT "$DL_dc_median[$bi] $UL_dc_median[$bi]\n");
}
GMT_psxy('-Sc0.12c -G46/139/87 -W1,blue,-');
for (my($bi)=0; $bi<@DL_uc_median; $bi++) {
next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
my($depth) = ($bi+0.5)*$opt_o;
last if ($depth > $opt_s);
print(GMT "$DL_uc_median[$bi] $UL_uc_median[$bi]\n");
}
GMT_psxy('-Sc0.12c -G46/139/87');
for (my($bi)=0; $bi<@DL_uc_median; $bi++) {
next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
my($depth) = ($bi+0.5)*$opt_o;
next unless ($depth > $opt_s);
print(GMT "$DL_uc_median[$bi] $UL_uc_median[$bi]\n");
}
GMT_pstext('-Gblue -N');
if (defined($outfile)) { printf(GMT "%f %f 14 0 0 TL $outfile [$P{run_label}]\n",-$mwm,1.1*$mwm); }
else { printf(GMT "%f %f 14 0 0 TL %03d\n [$P{run_label}]",$id,-$mwm,1.1*$wmw); }
GMT_psbasemap('-Bf0.01a0.05:"DL Vertical Velocity [m/s]":/f0.01a0.05:"UL Vertical Velocity [m/s]":WeSn');
GMT_end();
} # if dual_head
}
&antsExit(0);