#!/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: Tue Oct 13 10:51:31 2015
#                    (c) 2015 A.M. Thurnherr
#                    uE-Info: 50 0 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]

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

&antsUsage('b:i:k:o:p:',1,
	'[profile -i)d <id>]',
	'[-o)utput bin <resolution>] [-k) require <min> samples]',
	'[-p)lot <[%03d_wprof.eps]> [-b)t <wprof>]]',
	'<.wsamp file> [.wsamp file]');

$dual_head = (@ARGV==1);

if (defined($opt_o)) {
	$opt_o_override = 1;
	&antsCardOpt($opt_o);
} else {
	$opt_o = &antsRequireParam('outgrid_dz');
}

if (defined($opt_k)) {
	$opt_k_override = 1;
	&antsCardOpt($opt_k);
} else {
	$opt_k = &antsRequireParam('outgrid_minsamp');
	$opt_k *= 2 if ($dual_head);
}
	
#----------------------------------------------------------------------
# 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: -b only makes sense when a plot is 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)
#	good_bins(first_good,last_good)						good_bins(2,'*')
#	
#----------------------------------------------------------------------

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 good_bins($$)
{
	push(@gbFirst,shift); push(@gbLast,shift);
	$gbFirst[$#gbFirst] = 1    if ($gbFirst[$#gbFirst] eq '*');
	$gbLast[$#gbLast]   = 9e99 if ($gbLast[$#gbLast] eq '*');
}	
	
	

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

#----------------------------------------------------------------------
# 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 numberp($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,'-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 = '';												# 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}) && !$opt_k_override &&
					(( $dual_head && $P{outgrid_minsamp}*2!=$opt_k) ||
					 (!$dual_head && $P{outgrid_minsamp}!=$opt_k));
		
		$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 = $id; $RUN = antsRequireParam('run_label');
		undef(@rngMin); undef(@rngMax); undef(@bins);
		unless ($return = do "./EditParams") {					# man perlfunc
			croak("./EditParams: $@\n") if ($@);
#	        croak("do ./EditParams: $!\n") unless defined($return);		# happens when EditParams does not exist
#	        croak("couldn't run ./EditParams\n") unless $return;		# happens when EditParams does not exist
		}

		if ($opt_p && @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 ($opt_p && @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)?$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]}));
	&antsOut(@{$out[$bi]});				 
}

if ($opt_p) {																		# complete plot
	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();
}

&antsExit(0);