#!/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: Fri Jul 23 14:03:05 2021
#                    (c) 2015 A.M. Thurnherr
#                    uE-Info: 740 38 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
#	Mar 16, 2016: - adapted to gmt5
#   Mar 31, 2016: - changed version %PARAM
#	Apr 14, 2016: - added profile id to warning messages
#	May 24, 2016: - improved plot
#	May 26, 2016: - added automatic directory creation on -d
#	Nov 28, 2017: - removed wcorr plot
#	Dec  9, 2017: - added $antsSuppressCommonOptions = 1;
#	Oct 31, 2018: - improved label (no longer explained/residual stddev)
#	Nov  1, 2018: - made layout consistent for dual- and single-head profiles
#	Jul  1, 2021: - made %PARAMs more standard
#				  - added -z to remove biases
#				  - added %?c_w_diff.rms to output
#				  - BUG: dc_corr returned strange rms
#	Jul  7, 2021: - reversed logic of -z (enables bias correction by default)
#				  - BUG: plot had label in wrong location for single-head profiles
#	Jul  9, 2021: - added window correlation stats
#				  - updated %PARAM names
#	Jul 13, 2021: - BUG: dc_sig, dc_rms confusion
#	Jul 23, 2021: - added summary info to plot
#				  - added seabed to plot
#				  - added annotations to plot
# HISTORY END


($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);

$antsSuppressCommonOptions = 1;
&antsUsage('b:d:i:k:l:o:p:v:w:z',1,
	'[profile -i)d <id>]',
	'[disable -z)eroing of <w> (disable bias correction)]',
	'[-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 exists
croak("$0: no profile_id in first file => -i required\n")
	unless defined($id);
	
if (defined($opt_d)) {													# select output directory
	unless (-d $opt_d) {
	    unless ($opt_d =~ m{/}) {
			print(STDERR "Warning: Creating output sub-directory ./$opt_d\n");
		    mkdir($opt_d);
		}
		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" : '%03d_wprof.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;
}

#----------------------------------------------------------------------
# Correlation Statistics
#
#	return values:
#		   R	correlation coefficient w_DL,w_UL
#	     var	variance of avg(w_DL,w_UL)
#	  DL_var	variance of w_DL
#	  UL_var	variance of w_DL
#	rms_diff	rms of w_DL-w_UL
#----------------------------------------------------------------------

sub dc_corr($$)
{
	my($fi,$li) = @_;
	my($n) = 0;
	my($ax,$ay) = (0,0);
	my($ssq_diff) = 0;

	for (my($bi)=$fi; $bi<=$li; $bi++) {
		next unless numberp($DL_dc_median[$bi]) && numberp($UL_dc_median[$bi]);
		$n++;
		$ax += $DL_dc_median[$bi];
		$ay += $UL_dc_median[$bi];
		$ssq_diff += ($DL_dc_median[$bi] - $UL_dc_median[$bi])**2;
	}
	return (nan,nan,nan,nan,nan) unless ($n > 2);
	$ax /= $n;
	$ay /= $n;

	my($syy,$sxy,$sxx) = (0,0,0);
	for (my($bi)=$fi; $bi<=$li; $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($var) = ($sxx + $syy) / (2*$n);			# variance of avg(w_DL,w_UL)
	my($var_DL) = $sxx/$n;
	my($var_UL) = $syy/$n;
	my($rms_diff) = sqrt($ssq_diff/$n);
	return ($R,$var,$var_DL,$var_UL,$rms_diff);
}

sub uc_corr($$)
{
	my($fi,$li) = @_;
	my($n) = 0;
	my($ax,$ay) = (0,0);
	my($ssq_diff) = 0;

	for (my($bi)=$fi; $bi<=$li; $bi++) {
		next unless numberp($DL_uc_median[$bi]) && numberp($UL_uc_median[$bi]);
		$n++;
		$ax += $DL_uc_median[$bi];
		$ay += $UL_uc_median[$bi];
		$ssq_diff += ($DL_uc_median[$bi] - $UL_uc_median[$bi])**2;
	}
	return (nan,nan,nan,nan,nan) unless ($n > 2);
	$ax /= $n;
	$ay /= $n;

	my($syy,$sxy,$sxx) = (0,0,0);
	for (my($bi)=$fi; $bi<=$li; $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($var) = ($sxx + $syy) / (2*$n);			# variance of avg(w_DL,w_UL)
	my($var_DL) = $sxx/$n;
	my($var_UL) = $syy/$n;
	my($rms_diff) = sqrt($ssq_diff/$n);
	return ($R,$var,$var_DL,$var_UL,$rms_diff);
}

#----------------------------------------------------------------------
# Main Program
#----------------------------------------------------------------------

$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');
$instrument_type 	= &antsRequireParam('ADCP_type');
$xducer_frequency	= &antsRequireParam('ADCP_frequency');
$orientation		= &antsRequireParam('ADCP_orientation');

$dc_var				= &antsRequireParam('dc_w.var');				# for dual-head LADCPs, variables will be 
$uc_var				= &antsRequireParam('uc_w.var');				# overwritten by [du]c_corr()

($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($sumPF)>0);
}

my($R,$R2);
if (defined($opt_p)) {												# begin summary plot
	$xmin = -0.1; $x2min = -700;
	$xmax = 0.35; $x2max =	500;
	$ymin = antsParam('min_depth');
	$ymin = round($ymin-25,50);
	$ymax = antsParam('water_depth');
	$ymax = antsRequireParam('max_depth') unless numberp($ymax);
	$ymax = round($ymax+25,50);
	$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('-W0.5');
	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('-W0.5');
	print(GMT ">\n50 $ymin\n50 $ymax\n");
	print(GMT ">\n250 $ymin\n250 $ymax\n");
	print(GMT ">\n450 $ymin\n450 $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('-W1.5');
		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: #$id: inconsistent %%outgrid_dz in profile #$id")	# consistency checks
			if (defined($P{outgrid_dz}) && $P{outgrid_dz}!=$opt_o &&!$opt_o_override);
		&antsInfo("WARNING: inconsistent %%outgrid_minsamp in profile #$id")
			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(&antsRequireParam('ADCP_bin_length'))) {
			unless ($warned) {
				&antsInfo("WARNING: inconsistent ADCP sampling parameters in profile #$id --- using conservative values");
				$warned = 1;
			}
			$bin_length = min($bin_length,$P{ADCP_bin_length});
		}
		unless (round($pulse_length) == round(&antsRequireParam('ADCP_pulse_length'))) {
			unless ($warned) {
				&antsInfo("WARNING: inconsistent ADCP sampling parameters in profile #$id --- using conservative values");
				$warned = 1;
			}
			$pulse_length = min($pulse_length,$P{ADCP_pulse_length});
		}
		unless (round($blanking_dist) == round(&antsRequireParam('ADCP_blanking_distance'))) {
			unless ($warned) {
				&antsInfo("WARNING: inconsistent ADCP sampling parameters in profile #$id --- using conservative values");
				$warned = 1;
			}
			$blanking_dist = min($blanking_dist,$P{ADCP_blanking_distance});
		}

		$instrument_type2 	= &antsRequireParam('ADCP_type');					# for summary info
		$xducer_frequency2	= &antsRequireParam('ADCP_frequency');
		$orientation2		= &antsRequireParam('ADCP_orientation');

		$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('-W1,coral,-');
			for (my($bi)=0; $bi<=$#dcw1; $bi++) {
				printf(GMT "%f %f\n",$DL_dc_median[$bi],($bi+0.5)*$opt_o);
			}
			GMT_psxy('-W1,SeaGreen,-');
			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] - (!$opt_z ? $P{'dc_w.mu'} : 0));	# 	vertical velocity
		push(@{$dcw1[$bi]},$ants_[0][$dcwF] - (!$opt_z ? $P{'dc_w.mu'} : 0))	# 	single-instrument w
			if ($dual_head);					
		push(@{$dce[$bi]}, $ants_[0][$eF]);	 									# 	elapsed time
	} else {								 									# upcast
		push(@{$ucw[$bi]}, $ants_[0][$ucwF] - (!$opt_z ? $P{'uc_w.mu'} : 0));
		push(@{$ucw1[$bi]},$ants_[0][$ucwF] - (!$opt_z ? $P{'uc_w.mu'} : 0))
			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_var,$dc_var_DL,$dc_var_UL,$dc_rms_wdiff) = &dc_corr(int($opt_s/$opt_o),$#dcw1);		# correlation statistics
	($uc_R,$uc_var,$uc_var_DL,$uc_var_UL,$uc_rms_wdiff) = &uc_corr(int($opt_s/$opt_o),$#ucw1);
	&antsAddParams('dc_w.R',$dc_R,'uc_w.R',$uc_R,
				   'DL_dc_w.var',$dc_var_DL,'UL_dc_w.var',$dc_var_UL,
				   'DL_uc_w.var',$uc_var_DL,'UL_uc_w.var',$uc_var_UL,
				   'dc_w.var',$dc_var,'uc_w.var',$uc_var,
				   'dc_wdiff.rms',$dc_rms_wdiff,'uc_wdiff.rms',$uc_rms_wdiff);

	my($last_depth,$last_bi,$dc_sumsq_res,$dc_n,$uc_sumsq_res,$uc_n);			# window correlation
	for (my($bi)=0; $bi<=$#dcw1; $bi++) {
		($dc_R[$bi],$dc_var[$bi],$dummy,$dummy,$dc_rms_wdiff[$bi]) =
			&dc_corr(max(0,$bi-int(250/$opt_o)),min($#dcw1,$bi+int(250/$opt_o)));
		($uc_R[$bi],$uc_var[$bi],$dummy,$dummy,$uc_rms_wdiff[$bi]) =
			&uc_corr(max(0,$bi-int(250/$opt_o)),min($#dcw1,$bi+int(250/$opt_o)));
	}
		
#		my($depth) = ($bi+0.5) * $opt_o;
#		unless (defined($last_bi)) {
#			$last_bi = $bi;
#			$last_depth = $depth;
#			next;
#		}
#		if ($depth >= $last_depth+500 || $bi == $#dcw1) {
#			($dc_R[$bi],$dc_rms[$bi],$dc_rms_wdiff[$bi]) = &dc_corr($last_bi,$bi);
#			($uc_R[$bi],$uc_rms[$bi],$uc_rms_wdiff[$bi]) = &uc_corr($last_bi,$bi);
#			undef($last_bi);
#			next;
#		}
#	}

	if (defined($opt_p)) {														# plot 2nd-instrument profiles
		GMT_psxy('-W1,coral,.');
		for (my($bi)=0; $bi<=$#dcw1; $bi++) {
			printf(GMT "%f %f\n",$UL_dc_median[$bi],($bi+0.5)*$opt_o);
		}
		GMT_psxy('-W1,SeaGreen,.');
		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
#	- same output file layout for single- and dual-head systems
#----------------------------------------------------------------------

@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
				  'dc_w.R','dc_w.var','dc_wdiff.rms',
				  'uc_w.R','uc_w.var','uc_wdiff.rms');

if ($dual_head) {																# dual-head output
	&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('depth.min',round($min_depth),'depth.max',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);
}

#&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]}));
	if ($dual_head) {
		push(@{$out[$bi]},$dc_diff[$bi],$uc_diff[$bi],
						  $dc_R[$bi],$dc_var[$bi],$dc_rms_wdiff[$bi],
						  $uc_R[$bi],$uc_var[$bi],$uc_rms_wdiff[$bi]);
	} else {
		push(@{$out[$bi]},nan,nan,nan,nan,nan,nan);
	}
	&antsOut(@{$out[$bi]});				 
}

if (defined($opt_p)) {																# complete summary plot
	GMT_setR($R);

	if ($P{water_depth} > 0) {															# SEABED
		GMT_psxy('-G204/153/102');
		print(GMT "$xmin $ymax\n0.07 $ymax\n0.07 $P{water_depth}\n $xmin $P{water_depth}\n");
	}

	GMT_psxy('-W1.5,coral');														# 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('-W1.5,SeaGreen');
	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 -Gcoral');														# 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 -GSeaGreen');
	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('-W0.7,coral');
	for (my($bi)=0; $bi<=$#dcw; $bi++) {											# number of samples
		if ($dcns[$bi]) { printf(GMT "%f %f\n",$dcns[$bi],($bi+0.5)*$opt_o); }
		else 			{ print(GMT "nan nan\n"); }
	}
	GMT_psxy('-W0.7,SeaGreen');
	for (my($bi)=0; $bi<=$#dcw; $bi++) {
		if ($ucns[$bi]) { printf(GMT "%f %f\n",$ucns[$bi],($bi+0.5)*$opt_o); }
		else 			{ print(GMT "nan nan\n"); }
	}

	GMT_psbasemap('-Bf10a1000-950:"                                                  # of Samples":N');
	GMT_psbasemap('-Ba2000-1550N'); GMT_psbasemap('-Ba1000-750N');

	$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');

	if ($dual_head) {
		GMT_psxy('-W1,100/100/255');												# surface layer limit
			print(GMT "-0.1 $opt_s\n0.07 $opt_s\n");
	}
	
	GMT_unitcoords();
	my(@y) = (1.018,1.052,1.076,1.109);
	
	GMT_pstext('-F+f9,Helvetica,CornFlowerBlue+jTL -N');							# summary information
	if ($dual_head) {
		printf(GMT "0.64 $y[0] Dual-Head (%d / %d kHz)\n",
					round($xducer_frequency,50),round($xducer_frequency2,50));
	} else {
		printf(GMT "0.64 $y[0] %d kHz $instrument_type $orientation\n",
					round($xducer_frequency,50));
	}
	print( GMT "0.64 $y[1] rms <w>\n		0.77 $y[1] :\n");
	if ($dual_head) {
		printf(GMT "0.64 $y[2] rms @~D@~w\n		0.77 %f :\n",$y[2]+0.007);
		printf(GMT "0.64 %f correl. (r)\n		0.77 $y[3] :\n",$y[3]-0.005);
	}

	if ($dual_head) {
		if ($dc_rms_wdiff > sqrt($dc_var)) {
			GMT_pstext('-F+f9,Helvetica,coral+jTR -N -Gyellow');
		} else {
			GMT_pstext('-F+f9,Helvetica,coral+jTR -N');
		}
		if (numberp($dc_R)) {
			&antsInfo("WARNING: low dc correlation (r = %.1f) between UL and DL data in profile #$id",$dc_R)
				if ($dc_R < 0.3);
			printf(GMT "0.88 %f %.1fmm/s\n",$y[1]-0.005,round(sqrt($dc_var)*1000,.1));
			printf(GMT "0.88 %f %.1fmm/s\n",$y[2]+0.001,round($dc_rms_wdiff*1000,.1));
            printf(GMT "0.88 %f %.1f",$y[3]-0.005,$dc_R);
        } else {
			&antsInfo("WARNING: no overlap between UL and DL dc data below the surface layer in profile #$id");
        }
	} else {
		GMT_pstext('-F+f9,Helvetica,coral+jTR -N');
		printf(GMT "0.88 %f %.1fmm/s\n",$y[1]-0.005,round(sqrt($dc_var)*1000,.1));
	}

	if ($dual_head) {
		if ($uc_rms_wdiff > sqrt($uc_var)) {
			GMT_pstext('-F+f9,Helvetica,SeaGreen+jTR -N -Gyellow');
		} else {
			GMT_pstext('-F+f9,Helvetica,SeaGreen+jTR -N');
		}
		if (numberp($uc_R)) {
			&antsInfo("WARNING: low uc correlation (r = %.1f) between UL and DL data in profile #$id",$uc_R)
				if ($uc_R < 0.3);
			printf(GMT "0.99 %f %.1fmm/s\n",$y[1]-0.005,round(sqrt($uc_var)*1000,.1));
			printf(GMT "0.99 %f %.1fmm/s\n",$y[2]+0.001,round($uc_rms_wdiff*1000,.1));
            printf(GMT "0.99 %f %.1f",$y[3]-0.005,$uc_R);
        } else {
			&antsInfo("WARNING: no overlap between UL and DL uc data below the surface layer in profile #$id");
        }
	} else {
		GMT_pstext('-F+f9,Helvetica,SeaGreen+jTR -N');
		printf(GMT "0.99 %f %.1fmm/s\n",$y[1]-0.005,round(sqrt($uc_var)*1000,.1));
	}

	GMT_pstext('-F+f14,Helvetica,blue+jTL -N');											# annotations
	if (defined($outfile)) { print(GMT "0.01 -0.06 $outfile [$P{run_label}]\n"); }
	else { 					printf(GMT "0.01 -0.06 %03d\n [$P{run_label}]",$id); }
	GMT_pstext('-F+f12,Helvetica+jMR');
		print(GMT '0.62 0.98 m.abs.dev.');
	GMT_pstext('-F+f9,Helvetica,orange+jBR -N -Gwhite');
		print(GMT "0.99 0.99 V$VERSION\n");
	GMT_pstext('-F+f12,Helvetica,coral+jTL -Gwhite');
		print(GMT "0.02 0.02 downcast\n");
	GMT_pstext('-F+f12,Helvetica,SeaGreen+jTL -Gwhite');
		print(GMT "0.24 0.02 upcast\n");
	GMT_pstext('-F+f12,Helvetica+jBL -Gwhite');
		print(GMT "0.02 0.98 b.track\n");
        
	GMT_end();

	if ($dual_head && length($corrPF)>0) {												# 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('-W2,grey50');
			print(GMT "-$mwm -$mwm\n$mwm $mwm\n");
		GMT_psxy('-Sc0.12c -Gcoral -W0.3,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 -Gcoral');
			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 -GSeaGreen -W0.3,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 -GSeaGreen');
			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('-F+f14,Helvetica,blue+jTL -N');
			if (defined($outfile)) { printf(GMT "%f %f $outfile [$P{run_label}]\n",-$mwm,1.1*$mwm); }
		    else { 					 printf(GMT "%f %f %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 && length(corrPF) > 0 

}

&antsExit(0);