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