#!/usr/bin/perl

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

#                    L I S T B I N S 

#                    doc: Fri Aug 25 15:57:05 2006

#                    (c) 2006 A.M. Thurnherr

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

# Split data file into per-bin time series.

# HISTORY:

#	Aug 25, 2006: - created from [listEns]

#	Aug 26, 2006: - added -M)agdecl

#				  - changed -b to -f

#	Aug 27, 2006: - added %bin

#	Aug 28, 2006: - BUG: options were confused

#	Jan  4, 2007: - improved usage message for -a

#				  - added %mag_decl

#				  - BUG: roundoff error in %pct_good_vels

#	Sep 19, 2007: - adapted to new [RDI_BB_Read.pl]

#	Feb  7, 2008: - added sound-speed correction

#				  - enabled 3-beam solutions

#	Feb  8, 2008: - added -d)iscard <beam>

#				  - added -b)eam coordinate output

#	Feb 12, 2008: - modified 3-beam output

#				  - added -p)ct_good <min>

#	Feb 13, 2008: - various improvements

#	Feb 19, 2008: - BUG: division by zero

#					BUG: min() did not work with 1st elt undef

#	Feb 21, 2008: - BUG: had forgotten to undo debugging changes

#				  - removed missing magdecl warning on -b

#	Feb 22, 2008: - moved ssCorr() to [RDI_Utils.pl]

#				  - BUG: %d complete ensembles was written to STDOUT

#				  - BUG: %N_ensembles was total number of ensembles, not

#						 only the ones used

#				  - BUG: 0 errvel was output for 3-beam solutions =>

#						 wrong time-average statistics

#	May 19, 2009: - added -w to calculate vertical velocities from

#					two beam pairs separately

#	May 21, 2009: - added horizontal beampair velocities on -w

#				  - -P)itchRoll <bias/bias>

#	May 22, 2009: - added -B) <bias/bias/bias/bias>

#	May 23, 2009: - adapted to changed beampair-velocity fun name

#	Aug 22, 2010: - added -R

#	Apr 29, 2013: - cosmetics

#				  - added warning on missing -S

#	Aug  7, 2013: - BUG: -w did not respect -d

#   Mar  4, 2014: - added support for missing PITCH/ROLL/HEADING

#	Apr 22, 2014: - BUG: extraneous semicolon

#	Nov 24, 2014: - enabled -w always

#	Mar 22, 2015: - replaced -f by -o (allowing for pipes)

#	Jun 16, 2015: - BUG: velocity bias code did not respect bad velocities

#	Jan  5, 2016: - adapted to [ANTS_tools_lib.pl]

#				  - adapted to calculation of w12, w34 from earth-coordinate data

#				  - several other changes to the code that should not affect the results

#	Jan  6, 2016: - -b removed (always output beamvels)

#	Jan 31, 2016: - started debugging the obviously wrong Earth2Beam() transformation

#	Feb 29, 2016: - continued debugging; removed debugging code

#   Mar 17, 2016: - adapted to new Getopt library

#	Jun  9, 2016: - minor improvements

#				  - BUG: velBeamToEarth() has new interface

#	Jan 27, 2017: - BUG: earth-coordinate label was wrong

#				  - changed semantics to include all %-good fields of earth-coordinate data

#					and also report total pcg (3+4 beam) in min_pcg

#				  - BUG: 3-beam %ages were incorrect: 1) they were based on goodvels instead of the

#						 entire ensemble range; 2) pings-per-ensemble were not considered

#				  - BUG: output layout was all messed up for non-valid velocities

# General Notes:

#	- everything (e.g. beams) is numbered from 1

#	- no support for BT data

# Post-Processing:

#	- e.g. in a Makefile,

#		BFILT = | list -Lconv ensemble elapsed dn12='str2dec_time(\$$date,\$$time,2012)' heading pitch roll temp u v w err_vel w12 w34 corr1 corr2 corr3 corr4 amp1 amp2 amp3 amp4 > `basename $@ .b1`.b%d

#	- used as follows

#		listBins -o "${BFILT}"

#	- this writes only the interesting fields and calculates a 2012-referenced decimal day field

# Soundspeed Correction:

#	- applied as described in the RDI coord-trans manual

#	- sound-speed variation over range is ignored (valid for small gradients)

#	=> - same simple correction for all velocity components

#	   - simple correction for cell depths

# Min %-good (min_pcg):

#	- nan for records w/o valid velocities

#	- min(%-good) of the beams used for the velocity solution

#	- min_pcg does not have to decrease monotonically with distance,

#	  at least when 3-beam solutions are allowed and when -p is used to

#	  edit the data

#	- non-monotonic min_pcg is particularly obvious with the DYNAMUCK BM_ADCP

#	  data, where one of the beams performed much worse than the others

# %-good of earth-coordinate data:

#	pcg1	3-beam (correlation threshold)

#	pcg2	bad (errvel threshold)

#	pcg3	bad (>1 bad beam)

#	pcg4	4-beam

use Getopt::Std;

$ADCP_tools_minVersion = 1.4;

($ADCP_TOOLS) = ($0 =~ m{(.*/)[^/]+});

require "$ADCP_TOOLS/ADCP_tools_lib.pl";

die("Usage: $0 [-r)ange <first_ens,last_ens>] [-R)enumber ensembles] " .

			  "[-o)utput <redirection[>bin%d.raw]>] " .

			  "[output -a)ll ens (not just those with good vels)] " .

			  "[-M)agnetic <declination>] " .

			  "[-S)oundspeed correction <salin|*,temp|*,depth|*> " .

			  "[-P)itch/Roll <bias/bias>] [-B)eamvel <bias/bias/bias/bias>] " .

		 	  "[require -4)-beam solutions] [-d)iscard <beam#>] " .

		 	  "[-p)ct-good <min>] " .

			  "<RDI file>\n")

	unless (&getopts("4aB:d:M:o:p:r:P:RS:") && @ARGV == 1);

($P{pitch_bias},$P{roll_bias}) = split('[,/]',$opt_P);

($P{velbias_b1},$P{velbias_b2},$P{velbias_b3},$P{velbias_b4}) = split('[,/]',$opt_B);

die("$0: -4 and -d are mutually exclusive\n")

	if ($opt_4 && defined($opt_d));

$opt_p = 0 unless defined($opt_p);

$RDI_Coords::minValidVels = 4 if ($opt_4);			# no 3-beam solutions

print(STDERR "WARNING: magnetic declination not set!\n")

	unless defined($opt_M);

$opt_o = '>bin%d.raw' unless defined($opt_o);

$ifn = $ARGV[0];

($first_ens,$last_ens) = split(',',$opt_r)

	if defined($opt_r);

if (defined($opt_S)) {

	($SS_salin,$SS_temp,$SS_depth) = split(',',$opt_S);

	$variable_ssCorr = ($SS_salin eq '*' || $SS_temp eq '*' || $SS_depth eq '*');

} else {

	print(STDERR "WARNING: no soundspeed correction applied!\n");

}

#----------------------------------------------------------------------

sub min(@)								# return minimum

{

	my($min) = 99e99;

	for (my($i)=0; $i<=$#_; $i++) {

		$min = $_[$i] if defined($_[$i]) && ($_[$i] < $min);

	}

	return ($min == 99e99) ? nan : $min;

}

sub dumpBin($$$)						# write time series of single bin

{

	my($b,$fe,$le) = @_;

	my($out) = sprintf($opt_o,$b+1);

	open(P,"$out") || die("$out: $!\n");

	print(P "#ANTS#PARAMS# ");

	foreach my $k (keys(%P)) {

		print(P "$k\{$P{$k}\} ");

	}

	my($pct3b) = 100*$three_beam[$b]/(($le-$fe+1) * $dta{PINGS_PER_ENSEMBLE});

	printf(STDERR "%02d:%.0f%%/%.0f%% ",$b+1,100*$good_vels[$b]/($le-$fe+1),$pct3b);

	printf(P "pct_3_beam{%.0f} ",$pct3b);

	printf(P "pct_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1));

	printf(P "bin{%d}",$b+1);

	printf(P " soundspeed_correction{%s}",defined($opt_S) ? $opt_S : 'NONE!');

	printf(P " dz{%g}",$dz[$b] *

				(defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$fe],$SS_salin,$SS_temp,$SS_depth) : 1)

	) unless ($variable_ssCorr);

	print( P "\n");

	print(P "#ANTS#FIELDS# " .

			"{ensemble} {date} {time} {elapsed} " .

			"{heading} {pitch} {roll} " .

			"{sig_heading} {sig_pitch} {sig_roll} " .

			"{xmit_current} {xmit_voltage} " .

			"{temp} " .

			"{bv1} {bv2} {bv3} {bv4} {u} {v} {w} {err_vel} " .

			"{v12} {w12} {v34} {w34} " .

			"{corr1} {corr2} {corr3} {corr4} " .

			"{amp1} {amp2} {amp3} {amp4} " .

			"{pcg1} {pcg2} {pcg3} {pcg4} {3_beam} {min_pcg}"

	);

	print(P " {dz}") if ($variable_ssCorr);

	print(P "\n");

	my($t0) = $dta{ENSEMBLE}[$fe]->{UNIX_TIME};

	for (my($e)=$fe; $e<=$le; $e++) {

		next unless ($opt_a || $dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]);

		my($ssCorr) = defined($opt_S) ? ssCorr($dta{ENSEMBLE}[$e],$SS_salin,$SS_temp,$SS_depth) : 1;

		print(P "$dta{ENSEMBLE}[$e]->{NUMBER} ");

		print(P "$dta{ENSEMBLE}[$e]->{DATE} ");

		print(P "$dta{ENSEMBLE}[$e]->{TIME} ");

		printf(P "%d ",$dta{ENSEMBLE}[$e]->{UNIX_TIME}-$t0);		# elapsed time

		print(P defined($dta{ENSEMBLE}[$e]->{HEADING}) ? "$dta{ENSEMBLE}[$e]->{HEADING} " : 'nan ');

		print(P defined($dta{ENSEMBLE}[$e]->{PITCH}) ? "$dta{ENSEMBLE}[$e]->{PITCH} " : 'nan ');

		print(P defined($dta{ENSEMBLE}[$e]->{ROLL}) ? "$dta{ENSEMBLE}[$e]->{ROLL} " : 'nan ');

		print(P "$dta{ENSEMBLE}[$e]->{HEADING_STDDEV} ");

		print(P "$dta{ENSEMBLE}[$e]->{PITCH_STDDEV} ");

		print(P "$dta{ENSEMBLE}[$e]->{ROLL_STDDEV} ");

		print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_CURRENT} ");

		print(P "$dta{ENSEMBLE}[$e]->{ADC_XMIT_VOLTAGE} ");

		print(P "$dta{ENSEMBLE}[$e]->{TEMPERATURE} ");

		if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {

			printf(P "%g %g %g %g ",@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});

			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] * $ssCorr);

			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] * $ssCorr);

			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] * $ssCorr);

			if ($dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]) {

				print(P "nan ");

			} else {

				printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr);

			}

			printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]) ? "%g " : "nan ",

						$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][0]);

			printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]) ? "%g " : "nan ",

						$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][1]);

			printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]) ? "%g " : "nan ",

						$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][2]);

			printf(P defined($dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]) ? "%g " : "nan ",

						$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b][3]);

		} else {

			print(P "nan nan nan nan ");

			print(P "nan nan nan nan ");

			print(P "nan nan nan nan ");

		}

		print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");

		print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");

		if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {

			print(P "@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]} ");

		} else {

			print(P "nan nan nan nan ");

		}

		if ($beamCoords) {

			printf(P "%d ",$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]);

			printf(P "%s ",min(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}));

		} else {

			if ($dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b]) {

				print(P "$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0] ");

			} else {

				print(P "nan ");

			}

			printf(P "%g ",$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]+$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][3]);

		}

		printf(P "%g ",$dz[$b]*$ssCorr) if ($variable_ssCorr);

		print(P "\n");

	}

	close(P);

}

#----------------------------------------------------------------------

# MAIN

#----------------------------------------------------------------------

$P{RDI_file} = $ifn;

$P{mag_decl} = $opt_M if defined($opt_M);

readData($ifn,\%dta);

printf(STDERR "%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}));

$dta{HEADING_BIAS} = -$opt_M;						# magnetic declination

if ($dta{BEAM_COORDINATES}) {						# coords

	$beamCoords = 1;

} else {

	die("$ifn: only beam and earth coordinates supported\n")

		if (!$dta{EARTH_COORDINATES});

	die("$ifn: -p requires beam-coordinate data\n")

		if ($opt_p > 0);

	die("$ifn: -d requires beam-coordinate data\n")

		if defined($opt_d);

	die("$ifn: -B requires beam-coordinate data\n")

		if defined($opt_B);

}

for (my($b)=0; $b<$dta{N_BINS}; $b++) {				# calc dz

	$dz[$b] = $dta{DISTANCE_TO_BIN1_CENTER} + $b*$dta{BIN_LENGTH};

}

$lastGoodBin = 0;

for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {				# check/transform velocities

	$dta{ENSEMBLE}[$e]->{NUMBER} = $e+1 if ($opt_R);	# renumber ensembles

	next if (defined($first_ens) &&

			 $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);

	$dta{ENSEMBLE}[$e]->{PITCH} -= $P{pitch_bias}

		if defined($dta{ENSEMBLE}[$e]->{PITCH});

	$dta{ENSEMBLE}[$e]->{ROLL}  -= $P{roll_bias}

		if defined($dta{ENSEMBLE}[$e]->{ROLL});

	$P{first_ens} = $dta{ENSEMBLE}[$e]->{NUMBER},$fe = $e

		unless defined($P{first_ens});

	last if (defined($last_ens) &&

			 $dta{ENSEMBLE}[$e]->{NUMBER} > $last_ens);

	$P{last_ens} = $dta{ENSEMBLE}[$e]->{NUMBER};

	$le = $e;

	die("3-beams used in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")

		if ($dta{ENSEMBLE}[$e]->{N_BEAMS_USED} < 4);

	die("BIT error in ensemble $dta{ENSEMBLE}[$e]->{NUMBER}\n")

		if defined($dta{ENSEMBLE}[$e]->{BUILT_IN_TEST_ERROR});

	die("Low gain in ensemble #$dta{ENSEMBLE}[$e]->{NUMBER}\n")

        if ($dta{ENSEMBLE}[$e]->{LOW_GAIN});

	for (my($b)=0; $b<$dta{N_BINS}; $b++) {

		if ($beamCoords) {

			for (my($i)=0; $i<4; $i++) {									# percent-good editing (-p)

				if ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i] < $opt_p) {

					undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$i]);

					undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$i]);

				}

	        }

			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0] -= $P{velbias_b1}			# beam-velocity biases (-B)

				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0]);

			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1] -= $P{velbias_b2}

				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][1]);

			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2] -= $P{velbias_b3}

				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][2]);

			$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] -= $P{velbias_b4}

				if defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3]);

			if (defined($opt_d)) {											# discard data from given beam (-d)

				undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);

				undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);

			}

			@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} =					# save beam velocities

				@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]};

			@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =				# calculate w12, w34

				velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});

			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = 						# calculate earth velocities

				velBeamToEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});

			$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag;

			$three_beam[$b] += $RDI_Coords::threeBeamFlag;

			unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {		# not a valid transformation

				undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});

				next;

			}

		} else { 															# Earth coordinates

			@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]} =					# calculate beam velocities

				velEarthToBeam(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});

			@{$dta{ENSEMBLE}[$e]->{BEAMPAIR_VELOCITY}[$b]} =				# calculate w12, w34

				velBeamToBPEarth(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{BEAM_VELOCITY}[$b]});

			@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]} = 						# correct for heading bias

				velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]});

			$three_beam[$b] += ($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][0]/100 * $dta{PINGS_PER_ENSEMBLE});

			unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {		# no valid velocity

				undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});

				next;

			}

		}

		$dta{ENSEMBLE}[$e]->{GOOD_VEL}[$b] = 1;

		$good_vels[$b]++; 

		$lastGoodBin = $b if ($b > $lastGoodBin);

		$firstGoodEns = $e unless defined($firstGoodEns);

		$lastGoodEns = $e;

    }

}

unless (defined($opt_r)) {

	$fe = $firstGoodEns;

	$le = $lastGoodEns;

}

$P{N_ensembles} = $le - $fe + 1;

$firstBin = 0;

$lastBin = $lastGoodBin;

print( STDERR "Start      : $dta{ENSEMBLE}[$fe]->{DATE} $dta{ENSEMBLE}[$fe]->{TIME}\n");

print( STDERR "End        : $dta{ENSEMBLE}[$le]->{DATE} $dta{ENSEMBLE}[$le]->{TIME}\n");

printf(STDERR "Bins       : %d-%d\n",$firstBin+1,$lastBin+1);

if ($beamCoords) {

	printf(STDERR "3-Beam     : %d %d %d %d\n",$RDI_Coords::threeBeam_1,

											   $RDI_Coords::threeBeam_2,

											   $RDI_Coords::threeBeam_3,

											   $RDI_Coords::threeBeam_4);

}

print(STDERR "Good/3-Beam: ");

for ($b=$firstBin; $b<=$lastBin; $b++) {				# generate output

	dumpBin($b,$fe,$le);

}

print(STDERR "\n");

exit(0);