#!/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: Mon Sep 19 18:33:30 2022

#                    (c) 2015 A.M. Thurnherr

#                    uE-Info: 92 33 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

#	Aug  7, 2021: - BUG: some wsamp params were wrong

#	Oct 12, 2021: - reduced window thickness for correlations from 500m to 320m

#	Oct 13, 2021: - add filter to include correlation results only if there are

#				    less than 20% gaps

#	May 10, 2022: - added -f

#	May 17, 2022: - changed semantics to take output file name from input file names

#	May 18, 2022: - BUG: new semantics did not work (oops)

#	Sep 19, 2022: - BUG: plot did not work any more (unclear since when) for dual-head

#						 profiles

# 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:c:d:f: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]>]',

	'[-c)orrelation <window size[320m]>]',

	'[ouptput -d)ir <name> -f)ile <fmt[%03d.wprof]>]',

	'[-p)lot <[%03d_wprof.eps]> [-b)t <wprof>]]',

	'<DL.wsamp file> [UL.wsamp file] (or only <UL.wsamp file>)');

($basename) = ($ARGV =~ m{([^/]*)\.[^\.]*$});						# determine output file name

$opt_f = '%03d.wprof'													# output file name

	unless defined($opt_f);

$dual_head = (@ARGV==1);												# single or dual head

if ($dual_head) {

	my($basename2) = ($ARGV[0] =~ m{([^/]*)\.[^\.]*$});

	undef($basename) unless ($basename2 eq $basename);

}

$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

&antsCardOpt(\$opc_c,320);												# window thickness for correlation estimates

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) {

	if (defined($basename)) {

		$outfile = defined($opt_d) ? "$opt_d/$basename.wprof" 	   : "$basename.wprof";

		$opt_p 	 = defined($opt_d) ? "$opt_d/${basename}_wprof.ps" : "${basename}_wprof.ps"

			unless defined($opt_p);

	} else {

		$opt_p = defined($opt_d) ? "$opt_d/%03d_wprof.ps" : '%03d_wprof.ps'

			unless defined($opt_p);

		$outfile = defined($opt_d) ? sprintf("%s/$opt_f",$opt_d,$id)

								   : sprintf($opt_f,$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;

	return (nan,nan,nan,nan,nan) unless ($li > $fi);				# for shallow profiles this fails

	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 > 0.8*($li-$fi+1));

	die("$n,$li,$fi [$n > 0.8*($li-$fi+1)]\n") unless ($n > 0);

	$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('depth.min');

	$ymin = round($ymin-25,50);

	$ymax = antsParam('water_depth');

	$ymax = antsRequireParam('depth.max') 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'))) {