#!/usr/bin/perl
#======================================================================
#                    L I S T B I N S 
#                    doc: Fri Aug 25 15:57:05 2006
#                    dlm: Thu Feb 21 15:05:40 2008
#                    (c) 2006 A.M. Thurnherr
#                    uE-Info: 32 57 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================

# 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

# General Notes:
#	- everything (e.g. beams) is numbered from 1
#	- no support for BT data

# 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

require "getopts.pl";
$0 =~ m{(.*/)[^/]+};
require "$1RDI_BB_Read.pl";
require "$1RDI_Coords.pl";
require "$1RDI_Utils.pl";

die("Usage: $0 [-r)ange <first_ens,last_ens>] " .
			  "[output -f)ile <pat[bin%d.raw]>] " .
			  "[-a)ll ens (not just those with good vels)] " .
			  "[-M)agnetic <declination>] " .
			  "[-S)oundspeed correction <salin|*,temp|*,depth|*> " .
		 	  "[require -4)-beam solutions] [-d)iscard <beam#>] " .
		 	  "[-%)good <min>] " .
		 	  "[output -b)eam coordinates] " .
			  "<RDI file>\n")
	unless (&Getopts("4abd:f:M:p:r:S:") && @ARGV == 1);

die("$0: -4 and -d are mutually exclusive\n")
	if ($opt_4 && defined($opt_d));

die("$0: -p and -b are mutually exclusive\n")
	if ($opt_b && defined($opt_p));

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

$opt_f = 'bin%d.raw' unless defined($opt_f);
$ifn = $ARGV[0];

($first_ens,$last_ens) = split(',',$opt_r)
	if defined($opt_r);

($salin,$temp,$depth) = split(',',$opt_S)
	if defined($opt_S);
$variable_ssCorr = ($salin eq '*' || $temp eq '*' || $depth eq '*');

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

sub ssCorr($$$$)						# sound velocity correction factor
{
	my($e,$S,$T,$D) = @_;
	$S = $dta{ENSEMBLE}[$e]->{SALINITY} if ($S eq '*');
	$T = $dta{ENSEMBLE}[$e]->{TEMPERATURE} if ($T eq '*');
	if ($D eq '*') {
		print(STDERR "WARNING: soundspeed correction using instrument pressure instead of depth!\n")
			unless ($warned);
		$warned = 1;
		$D = $dta{ENSEMBLE}[$e]->{PRESSURE};
	}
	return soundSpeed($S,$T,$D) / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
}
	
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($file) = sprintf($opt_f,$b+1);

	open(P,">$file") || die("$file: $!\n");
	print(P "#ANTS#PARAMS# ");
	foreach my $k (keys(%P)) {
		print(P "$k\{$P{$k}\} ");
	}
	my($pct3b) = ($good_vels[$b] > 0) ? 100*$three_beam[$b]/$good_vels[$b] : nan;
	printf(STDERR "%.0f%%/%.0f%% ",100*$good_vels[$b]/($le-$fe+1),$pct3b);
	printf(P "pct_good_vels{%.0f} ",100*$good_vels[$b]/($le-$fe+1));
	printf(P "pct_3_beam{%.0f} ",$pct3b);
	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($fe,$salin,$temp,$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} " .
			($opt_b ? "{v1} {v2} {v3} {v4} " : "{u} {v} {w} {err_vel} ") .
			"{corr1} {corr2} {corr3} {corr4} " .
			"{amp1} {amp2} {amp3} {amp4} " .
			($opt_b ? "{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($e,$salin,$temp,$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);
		print(P "$dta{ENSEMBLE}[$e]->{HEADING} ");
		print(P "$dta{ENSEMBLE}[$e]->{PITCH} ");
		print(P "$dta{ENSEMBLE}[$e]->{ROLL} ");
		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 ",$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);
			printf(P "%g ",$dta{ENSEMBLE}[$e]->{VELOCITY}[$b][3] * $ssCorr);
		} else {
			print(P "nan nan nan nan ");
		}
		print(P "@{$dta{ENSEMBLE}[$e]->{CORRELATION}[$b]} ");
		print(P "@{$dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$b]} ");
		if ($opt_b) {
			print(P "@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]} ");
		} else {
			printf(P "%d ",$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b]);
			printf(P "%s ",min(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]}));
		}
		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("%d complete ensembles...\n",scalar(@{$dta{ENSEMBLE}}));
$dta{HEADING_BIAS} = -$opt_M;						# magnetic declination

if ($dta{BEAM_COORDINATES}) {						# coords
	$beamCoords = 1;
} else {
	die("$0: -b requires input in beam coordinates\n")
		if ($opt_b);
	die("$ifn: only beam and earth coordinates implemented so far\n")
		if (!$dta{EARTH_COORDINATES});
}

$P{N_ensembles} = @{$dta{ENSEMBLE}};

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
	next if (defined($first_ens) &&
			 $dta{ENSEMBLE}[$e]->{NUMBER} < $first_ens);
	$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 (defined($opt_d)) {
			undef($dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b][$opt_d-1]);
			undef($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][$opt_d-1]);
		}
		for (my($i)=0; $i<4; $i++) {
			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]} = $beamCoords
			? velInstrumentToEarth(\%dta,$e,
				  velBeamToInstrument(\%dta,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
			  )
			: velApplyHdgBias(\%dta,$e,@{$dta{ENSEMBLE}[$e]->{VELOCITY}[$b]})
				unless ($opt_b);
		unless (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$b][0])) {
			undef(@{$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[$b]});
			next;
		}
		$dta{ENSEMBLE}[$e]->{THREE_BEAM}[$b] = $RDI_Coords::threeBeamFlag;
		$three_beam[$b] += $RDI_Coords::threeBeamFlag;
		$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;
}

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