#!/usr/bin/perl
#======================================================================
#                    L A D C P _ W _ R E G R I D 
#                    doc: Fri Apr 24 17:15:59 2015
#                    dlm: Mon Jul 27 17:20:23 2015
#                    (c) 2015 A.M. Thurnherr
#                    uE-Info: 182 21 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 %outgrid_*
#	Jul 27, 2015: - B-t)rack <wprof>

($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('b:i:k:o:p:t:',1,
	'[profile -i)d <id>]',
	'[-o)utput bin <resolution[10m]>] [-k) require <min> samples]',
#	'[valid LADCP -b)ins <bin[2],bin[*][,bin[2],bin[*]]>',
	'[-p)lot <[%03d_wprof.eps]> [B-t)rack <wprof>]]',
	'<.samp file> [.samp file]');

$dual_head = (@ARGV==1);

$opt_o_override = defined($opt_o);
$opt_o = $P{outgrid_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);

if (defined($opt_b)) {
	croak("$0: -b not implemented yet\n");
	($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b);
	croak("$0: cannot decode -b $opt_b\n")
		unless (numberp($LADCP_firstBin) &&
				($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));
}

#----------------------------------------------------------------------
# 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");
}

croak("$0: -t only makes sense when a plot is produced\n")
	if $opt_t && !defined($opt_p);

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

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

	if (defined($opt_t)) {
		open(BT,$opt_t) || croak("$opt_t: $!\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_t: 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 = '';												# 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 %outgrid_dz in input files")
			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}) && $P{outgrid_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);
		}
	}
	next if ($ants_[0][$bF]<$P{outgrid_firstbin}-1 || $ants_[0][$bF]>$P{outgrid_lastbin}-1);
	$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('outgrid_dz',$opt_o,'outgrid_minsamp',$opt_k);
	&antsAddParams('outgrid_firstbin',$P{outgrid_firstbin},'outgrid_lastbin',$P{outgrid_lastbin});
	&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);