#!/usr/bin/perl

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

#                    L I S T W 

#                    doc: Wed Mar 24 06:45:09 2004

#                    dlm: Mon Apr  2 18:32:03 2018

#                    (c) 2004 A.M. Thurnherr

#                    uE-Info: 25 61 NIL 0 0 72 0 2 4 NIL ofnI

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

# dump vertical velocities

# NB: currently broken

# HISTORY:

#	Mar 24, 2004: - created from [listens] and [mkprofile]

#	Mar 27, 2004: - added elapsed field

#				  - floatized time

#	Apr  3, 2004: - cosmetics

#	Nov  8, 2005: - UNIXTIME => UNIX_TIME

#	Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)

#	Jul 30, 2009: - NaN => nan

#   Nov 25, 2013: - checkEnsemble() expunged

#	Mar 17, 2016: - removed warning

#				  - updated ancient library names

#	Apr  2, 2018: - BUG: velBeamToInstrument() used old usage

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

require "$1/RDI_PD0_IO.pl";

require "$1/RDI_Coords.pl";

require "$1/RDI_Utils.pl";

use Getopt::Std;

$USAGE = "$0 @ARGV";

die("Usage: $0 " .

		"[-A)nts] " .

		"[-F)ilter <script>] " .

		"[bin -r)ange  <bin|0,bin|*>] " .

		"[-e)rr-vel <max|0.1>] [-c)orrelation <min|70>] " .

		"[-S)alin <val|35>] [-t)emp <bias>] " .

		"[output -f)ields <field[,...]> " .

		"<RDI file>\n")

	unless (&getopts("Ac:e:F:f:r:S:t:") && @ARGV == 1);

$opt_e = 0.1 unless defined($opt_e);				# defaults

$opt_c = 70	 unless defined($opt_c);

$opt_S = 35  unless defined($opt_S);

print(STDERR "WARNING: Using uncalibrated ADCP temperature!\n"),$opt_t = 0

	 unless defined($opt_t);

require $opt_F if defined($opt_F);					# load filter

if ($opt_f) {										# additional fields

	@f = split(',',$opt_f);

	foreach $f (@f) {

		my($f) = $f;								# copy it

		$f =~ s{\[.*$}{};							# remove indices

		$addFields .= " {$f}";

	}

}

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

print(STDERR "Reading $ARGV[0]...");				# read data

readData($ARGV[0],\%dta);

print(STDERR "done\n");

if (defined($opt_r)) {								# bin range

	($minb,$maxb) = split(',',$opt_r);

	die("$0: can't decode -r $opt_r\n") unless defined($maxb);

} else {

	$minb = 0;

	$maxb = $dta{N_BINS} - 1;

}

die("$ARGV[0]: not enough bins for choice of -r\n")	# enough bins?

	unless ($dta{N_BINS} >= $maxb);

if ($dta{BEAM_COORDINATES}) {						# coords used

	$beamCoords = 1;

} elsif (!$dta{EARTH_COORDINATES}) {

	die("$ARGV[0]: only beam and earth coordinates implemented so far\n");

}

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

# Reference-Layer w (from [mkprofile])

#	- also sets W field when valid

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

sub ref_lr_w($)

{

	my($ens) = @_;

	my($i,$n,$w);

	for ($i=$minb; $i<=$maxb; $i++) {

		next if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $opt_c ||

				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $opt_c ||

				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $opt_c ||

				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $opt_c);

		if ($beamCoords) {

			next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100 ||

					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100 ||

					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100 ||

					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);

			@v = velInstrumentToEarth(\%dta,$ens,

					velBeamToInstrument(\%dta,$ens,

						@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));

		} else {

			next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||

					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||

					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||

					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);

			@v = @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]};

		}

		next if (!defined($v[3]) || abs($v[3]) > $opt_e);

		if (defined($v[2])) {							# valid w

			$dta{ENSEMBLE}[$ens]->{W}[$i] = $v[2];

			$w += $v[2]; $n++;

		}

	}

	return $n ? $w/$n : undef;

}

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

print(STDERR "Generating profile by integrating w...");

for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {

	filterEnsemble(\%dta,$e)								# filter ensemble 

		if (defined($opt_F) &&

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

	$dta{ENSEMBLE}[$e]->{REFW} = ref_lr_w($e);

	next unless defined($dta{ENSEMBLE}[$e]->{REFW});

	unless (defined($firstgood)) {							# init profile

		$firstgood = $lastgood = $e;			

		$depth = 0;

	}

	my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} -			# time step since

			  $dta{ENSEMBLE}[$lastgood]->{UNIX_TIME};		# ... last good ens

	if ($dt > 120) {

		printf(STDERR "\nWARNING: %d-s gap too long, profile restarted at ensemble $e...",$dt);

		$firstgood = $lastgood = $e;			

		$dt = $depth = 0;

	}

	$depth += $dta{ENSEMBLE}[$lastgood]->{REFW} * $dt		# integrate depth

		if ($dt > 0);

	$dta{ENSEMBLE}[$e]->{DEPTH} = $depth;

	$atbottom = $e, $maxdepth = $depth if ($depth > $maxdepth);	

	my($ss) = soundSpeed($opt_S,							# sound-speed corr

				  	     $dta{ENSEMBLE}[$e]->{TEMPERATURE}-$opt_t,

				   	     $depth);

	$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} =

		$ss / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};

	$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND} = $ss;

	$dta{ENSEMBLE}[$e]->{REFW} *=

		$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION};

	$lastgood = $e;

}

printf(STDERR "done (max depth = %.1fm, depth at end of cast = %.1fm)\n",

			$maxdepth,$depth);

filterEnsembleStats() if defined($opt_F);

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

print(STDERR "Writing output...");

if ($opt_A) {

	print("#ANTS# [] $USAGE\n");

	print("#ANTS#FIELDS# {ens} {unix_time} {time} {bin} {depth} {dz} {w} {ref_w} {dw} $addFields\n");

	printf("#ANTS#PARAMS# maxdepth{$max_depth} bottom_time{%d}\n",

	    $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME} - 

		 	$dta{ENSEMBLE}[$firstgood]->{UNIX_TIME});

} else {

	print("# ens-no time elapsed bin-no depth dz w ref-w dw $addFields\n");

	print("#----------------------------------------------------------------------\n");

}

for ($e=$firstgood; $e<=$lastgood; $e++) {

	for ($i=$minb; $i<=$maxb; $i++) {						# dump valid

		next unless defined($dta{ENSEMBLE}[$e]->{W}[$i]);

		printf("%d %f %f %d %.1f %.1f %g %g %g ",

			$e,$dta{ENSEMBLE}[$e]->{UNIX_TIME},

			$dta{ENSEMBLE}[$e]->{UNIX_TIME} -

				$dta{ENSEMBLE}[$firstgood]->{UNIX_TIME},$i,

			$dta{ENSEMBLE}[$e]->{DEPTH} +

				$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *

					($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),

			$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *

				($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),

			$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *

				$dta{ENSEMBLE}[$e]->{W}[$i],

			$dta{ENSEMBLE}[$e]->{REFW},

			$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *

				$dta{ENSEMBLE}[$e]->{W}[$i] - $dta{ENSEMBLE}[$e]->{REFW},

		);

		sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }

		if (@f) {

			foreach $f (@f) {

				my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});

				$fn = $f unless defined($fn);

				p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));

			}

		}

		print("\n");

	}

}

print(STDERR "done\n");

exit(0);