#!/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

#                    (c) 2015 A.M. Thurnherr

#======================================================================

$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

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

croak("$0: no profile_id in first file => -i required\n")

	unless defined($id);

}

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

}

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

}

	$opt_o_override = 1;

	&antsCardOpt($opt_o);

} else {

	$opt_o = &antsRequireParam('outgrid_dz');

}

	$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 = -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($P{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($P{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($P{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});

		}

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

		}

	}